Linear Algebra

Original Tutorial may be found here.

In [1]:
import numpy as np

This Tutorial will use also the Linear Algebra (linalg) sub-module (sitting inside numpy).

Refer to the for full numpy.linalg utilities.

Simple Array Operations

In [2]:
a = np.array([[1.0, 2.0],
              [3.0, 4.0]])
In [3]:
a
Out[3]:
array([[ 1.,  2.],
       [ 3.,  4.]])
In [4]:
a.transpose()
Out[4]:
array([[ 1.,  3.],
       [ 2.,  4.]])
In [5]:
a
Out[5]:
array([[ 1.,  2.],
       [ 3.,  4.]])

Compute the inverse. This MathIsFun.com article provides some excellent examples in explaining inverse of a matrix.

In [6]:
a_inv = np.linalg.inv(a)
In [7]:
a_inv
Out[7]:
array([[-2. ,  1. ],
       [ 1.5, -0.5]])

We know that, multiplying (dot multiplication) a matrix and its inverse gives us an identify matrix (see the MathIsFun.com website!).

In [8]:
np.dot(a, a_inv)
Out[8]:
array([[  1.00000000e+00,   1.11022302e-16],
       [  0.00000000e+00,   1.00000000e+00]])

See? :)

In [9]:
u = np.eye(2) # unit 2x2 matrix; "eye" represents "I"  (Identify Matrix)
In [10]:
u
Out[10]:
array([[ 1.,  0.],
       [ 0.,  1.]])
In [11]:
j = np.array([[0.0, -1.0],
              [1.0, 0.0]])
In [12]:
j
Out[12]:
array([[ 0., -1.],
       [ 1.,  0.]])
In [13]:
np.dot(j, j) # matrix product
Out[13]:
array([[-1.,  0.],
       [ 0., -1.]])

Use the trace() function to return the sum along the diagonals of the array. See this NumPy Doc.

In [14]:
np.trace(u) # trace
Out[14]:
2.0

Say we have three matrices:

ax = y

in which a and y are known.

We want to solve x.

In [15]:
a = np.array([[1.0, 2.0],
              [3.0, 4.0]])
In [16]:
y = np.array([[5.], [7.]])
In [17]:
a
Out[17]:
array([[ 1.,  2.],
       [ 3.,  4.]])
In [18]:
y
Out[18]:
array([[ 5.],
       [ 7.]])
In [19]:
x = np.linalg.solve(a, y)
In [20]:
x
Out[20]:
array([[-3.],
       [ 4.]])

To prove that this is correct, this must be true:

ax = y
In [21]:
s = (np.dot(a,x) == y)
In [22]:
s
Out[22]:
array([[ True],
       [ True]], dtype=bool)
In [23]:
s.all()
Out[23]:
True

Or just visualize it:

In [24]:
np.dot(a,x)
Out[24]:
array([[ 5.],
       [ 7.]])
In [25]:
y
Out[25]:
array([[ 5.],
       [ 7.]])

It certain looks like:

ax = y

Now, eigenvalue...

In [26]:
j
Out[26]:
array([[ 0., -1.],
       [ 1.,  0.]])
In [27]:
np.linalg.eig(j)
Out[27]:
(array([ 0.+1.j,  0.-1.j]),
 array([[ 0.70710678+0.j        ,  0.70710678-0.j        ],
        [ 0.00000000-0.70710678j,  0.00000000+0.70710678j]]))

See this article to learn more about eigenvalue and eigenvector.

Matrix Class

Here is a short intro to the Matrix class.

In [28]:
A = np.matrix('1.0 2.0; 3.0 4.0')
In [29]:
A
Out[29]:
matrix([[ 1.,  2.],
        [ 3.,  4.]])
In [30]:
type(A)
Out[30]:
numpy.matrixlib.defmatrix.matrix
In [31]:
A
Out[31]:
matrix([[ 1.,  2.],
        [ 3.,  4.]])
In [32]:
A.T  # Transpose
Out[32]:
matrix([[ 1.,  3.],
        [ 2.,  4.]])
In [33]:
X = np.matrix('5.0 7.0')
In [34]:
X
Out[34]:
matrix([[ 5.,  7.]])
In [35]:
Y = X.T
In [36]:
Y
Out[36]:
matrix([[ 5.],
        [ 7.]])
In [37]:
A*Y   # matrix multiplication
Out[37]:
matrix([[ 19.],
        [ 43.]])
In [38]:
A.I   # inverse
Out[38]:
matrix([[-2. ,  1. ],
        [ 1.5, -0.5]])
In [39]:
np.linalg.solve(A,Y) # solving linear equation
Out[39]:
matrix([[-3.],
        [ 4.]])

Indexing: Comparing Matrices and 2D Arrays

Note that there are some important differences between NumPy arrays and matrices. NumPy provides two fundamental objects: an N-dimensional array object and a universal function object. Other objects are built on top of these. In particular, matrices are 2-dimensional array objects that inherit from the NumPy array object. For both arrays and matrices, indices must consist of a proper combination of one or more of the following: integer scalars, ellipses, a list of integers or boolean values, a tuple of integers or boolean values, and a 1-dimensional array of integer or boolean values. A matrix can be used as an index for matrices, but commonly an array, list, or other form is needed to accomplish a given task.

As usual in Python, indexing is zero-based. Traditionally we represent a 2D array or matrix as a rectangular array of rows and columns, where movement along axis 0 is movement across rows, while movement along axis 1 is movement across columns.

Let's make an array and matrix to slice:

In [40]:
A = np.arange(12)
In [41]:
A
Out[41]:
array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11])
In [42]:
A.shape = (3,4)
In [43]:
A
Out[43]:
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11]])
In [44]:
M = np.mat(A.copy())  # np.mat Interpret the input as a matrix.
In [45]:
M
Out[45]:
matrix([[ 0,  1,  2,  3],
        [ 4,  5,  6,  7],
        [ 8,  9, 10, 11]])
In [46]:
print type(A)," ",type(M)
<type 'numpy.ndarray'>   <class 'numpy.matrixlib.defmatrix.matrix'>
In [47]:
print A
[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]]
In [48]:
print M
[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]]

Now, let's take some simple slices. Basic slicing uses slice objects or integers.

For example, the evaluation of A[:] and M[:] will appear familiar from Python indexing, however it is important to note that slicing NumPy arrays does not make a copy of the data; slicing provides a new view of the same data.

In [49]:
print A[:]; print A[:].shape
[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]]
(3L, 4L)
In [50]:
print M[:]; print M[:].shape
[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]]
(3L, 4L)

Now for something that differs from Python indexing: you may use commaseparated indices to index along multiple axes at the same time.

In [51]:
print A[:,1]; print A[:,1].shape
[1 5 9]
(3L,)
In [52]:
print M[:,1]; print M[:,1].shape
[[1]
 [5]
 [9]]
(3L, 1L)

Notice the difference in the last two results.

Use of a single colon for the 2D array produces a 1-dimensional array, while for a matrix it produces a 2-dimensional matrix. A slice of a matrix will always produce a matrix.

For example, a slice M[2,:] produces a matrix of shape (1,4).

In contrast, a slice of an array will always produce an array of the lowest possible dimension.

For example, if C were a 3-dimensional array, C[...,1] produces a 2D array while C[1,:,1] produces a 1-dimensional array.

From this point on, we will show results only for the array slice if the results for the corresponding matrix slice are identical.

Lets say that we wanted the 1st and 3rd column of an array. One way is to slice using a list:

In [53]:
A
Out[53]:
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11]])
In [54]:
A[:,[1,3]]
Out[54]:
array([[ 1,  3],
       [ 5,  7],
       [ 9, 11]])
In [55]:
A[:,].take([1,3],axis=1)
Out[55]:
array([[ 1,  3],
       [ 5,  7],
       [ 9, 11]])

If we wanted to skip the first row, we could use:

In [56]:
A[1:,].take([1,3],axis=1)
Out[56]:
array([[ 5,  7],
       [ 9, 11]])

Or we could simply use A[1:,[1,3]]. Yet another way to slice the above is to use a cross product:

In [57]:
 A[1:,[1,3]]
Out[57]:
array([[ 5,  7],
       [ 9, 11]])
In [58]:
A
Out[58]:
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11]])

Now let's do something a bit more complicated. Lets say that we want to retain all columns where the first row is greater than 1. One way is to create a boolean index:

In [59]:
A[0,:] > 1
Out[59]:
array([False, False,  True,  True], dtype=bool)
In [60]:
A[:, A[0,:] > 1]
Out[60]:
array([[ 2,  3],
       [ 6,  7],
       [10, 11]])

Just what we wanted! But indexing the matrix is not so convenient.

In [61]:
M
Out[61]:
matrix([[ 0,  1,  2,  3],
        [ 4,  5,  6,  7],
        [ 8,  9, 10, 11]])
In [62]:
M[0,:] > 1
Out[62]:
matrix([[False, False,  True,  True]], dtype=bool)
In [63]:
M[:, M.A[0,:]>1]   # Note the M.A
Out[63]:
matrix([[ 2,  3],
        [ 6,  7],
        [10, 11]])

If we wanted to conditionally slice the matrix in two directions, we must adjust our strategy slightly. Instead of:

In [64]:
A[A[:,0]>2,A[0,:]>1]
Out[64]:
array([ 6, 11])
In [65]:
M[M.A[:,0]>2,M.A[0,:]>1]
Out[65]:
matrix([[ 6, 11]])

we need to use the cross product ix_():

In [66]:
A[np.ix_(A[:,0]>2,A[0,:]>1)]
Out[66]:
array([[ 6,  7],
       [10, 11]])
In [67]:
M[np.ix_(M.A[:,0]>2,M.A[0,:]>1)]
Out[67]:
matrix([[ 6,  7],
        [10, 11]])

Conclusion

We have now completed lesson 7 - Linear Algebra.