Linear algebra with NumPy

Linear algebra is an important subdivision of mathematics. We can use linear algebra, for instance, to perform linear regression. The numpy.linalg subpackage holds linear algebra routines. With this subpackage, you can invert matrices, compute eigenvalues, solve linear equations, and find determinants among other matters. Matrices in NumPy are represented by a subclass of ndarray.

Inverting matrices with NumPy

The inverse of a square and invertible matrix A in linear algebra is the matrix A-1, which when multiplied with the original matrix is equal to the identity matrix I. This can be written down as the following mathematical equation:

A A-1 = I

The inv() function in the numpy.linalg subpackage can do this for us. Let's invert an example matrix. To invert matrices, follow the ensuing steps:

1. Create the example matrix.

   We will create the demonstration matrix with the mat() function:

   A = np.mat("2 4 6;4 2 6;10 -4 18")
   print "A\n", A

   The A matrix is printed as follows:

   A
   [[ 2 4 6]
    [ 4 2 6]
    [10 -4 18]]]
   

2. Invert the matrix.

   Now, we can view the inv() subroutine in action:

   inverse = np.linalg.inv(A)
   print "inverse of A\n", inverse

   The inverse matrix is displayed as follows:

   inverse of A
   [[-0.41666667 0.66666667 -0.08333333]
    [ 0.08333333 0.16666667 -0.08333333]
    [ 0.25 -0.33333333 0.08333333]]
   

   Tip

   If the matrix is singular, or not square, a LinAlgError is raised. If you wish, you can check the solution manually. This is left as a drill for you. The pinv() NumPy function performs a pseudo inversion, which can be applied to any matrix, including matrices that are not square.

3. Check by multiplication.

   Let's check what we get when we multiply the original matrix with the result of the inv() function:

   print "Check\n", A * inverse

   The result is the identity matrix, as expected (ignoring small differences):

   Check
   [[ 1.00000000e+00 0.00000000e+00 -5.55111512e-17]
    [ -2.22044605e-16 1.00000000e+00 -5.55111512e-17]
    [ -8.88178420e-16 8.88178420e-16 1.00000000e+00]]
   

By subtracting the 3 x 3 identity matrix from the previous result, we get the errors of the inversion process:

print "Error\n", A * inverse - np.eye(3)

The errors should be negligible in general, but in some cases small errors could be propagated with undesirable side effects:

[[ -1.11022302e-16 0.00000000e+00 -5.55111512e-17]
 [ -2.22044605e-16 4.44089210e-16 -5.55111512e-17]
 [ -8.88178420e-16 8.88178420e-16 -1.11022302e-16]]

In such cases, higher precision data types might help or switch to a superior algorithm. We computed the inverse of a matrix with the inv() routine of the numpy.linalg subpackage. We made certain, with matrix multiplication, whether this is indeed the inverse matrix (see inversion.py in this book's code bundle):

import numpy as np

A = np.mat("2 4 6;4 2 6;10 -4 18")
print "A\n", A

inverse = np.linalg.inv(A)
print "inverse of A\n", inverse

print "Check\n", A * inverse
print "Error\n", A * inverse - np.eye(3)

Solving linear systems with NumPy

A matrix transforms a vector into another vector in a linear fashion. This operation numerically corresponds to a system of linear equations. The solve() subroutine of numpy.linalg solves systems of linear equations of the form Ax = b; here, A is a matrix, b can be a one-dimensional or two-dimensional array, and x is an unknown quantity. We will witness the dot() subroutine in action. This function computes the dot product of two floating-point numbers' arrays.

Let's solve an instance of a linear system. To solve a linear system, follow the ensuing steps:

1. Create the matrix A and array b.

   The following code will create A and b:

   A = np.mat("1 -2 1;0 2 -8;-4 5 9")
   print "A\n", A
   b = np.array([0, 8, -9])
   print "b\n", b

   The matrix A and array (vector) b are defined as follows:

   
2. Call the solve() function.

   Solve this linear system with the solve() function:

   x = np.linalg.solve(A, b)
   print "Solution", x

   The solution of the linear system is as follows:

   Solution [ 29. 16. 3.]
   

3. Check with the dot() function.

   Check whether the solution is correct with the dot() function:

   print "Check\n", np.dot(A , x)

   The result is as expected:

   Check
   [[ 0. 8. -9.]]
   

We solved a linear system by applying the solve() function from the linalg subpackage of NumPy and checking the result with the dot() function (please refer to solution.py in this book's code bundle):

import numpy as np

A = np.mat("1 -2 1;0 2 -8;-4 5 9")
print "A\n", A

b = np.array([0, 8, -9])
print "b\n", b

x = np.linalg.solve(A, b)
print "Solution", x

print "Check\n", np.dot(A , x)