NumPy Matrix and Linear Algebra

What's difference between numpy dot() and inner()?

Let's look into 2D array as an example:

>>> a=np.array([[1,2],[3,4]])
>>> b=np.array([[11,12],[13,14]])
>>> np.dot(a,b)
array([[37, 40],
       [85, 92]])
>>> np.inner(a,b)
array([[35, 41],
       [81, 95]])

With dot():

With inner():

The chapters on NumPy have been using arrays (NumPy Array Basics A and NumPy Array Basics B). However, for certain areas such as linear algebra, we may instead want to use matrix.

>>> import numpy as np
>>> A = np.matrix([[1.,2], [3,4], [5,6]])
>>> A
matrix([[ 1.,  2.],
        [ 3.,  4.],
        [ 5.,  6.]])

We may also take the Matlab style by giving a string rather than a list:

>>> B = np.matrix("1.,2; 3,4; 5,6")
>>> B
matrix([[ 1.,  2.],
        [ 3.,  4.],
        [ 5.,  6.]])

Vectors are handled as matrices with one row or one column:

>>> x = np.matrix("10., 20.")
>>> x
matrix([[ 10.,  20.]])
>>> x.T
matrix([[ 10.],
        [ 20.]])

Here is an example for matrix and vector multiplication:

>>> x = np.matrix("4.;5.")
>>> x
matrix([[ 4.],
        [ 5.]])

>>> A = np.matrix([[1.,2], [3,4], [5,6]])
>>> A
matrix([[ 1.,  2.],
        [ 3.,  4.],
        [ 5.,  6.]])

>>> A*x
matrix([[ 14.],
        [ 32.],
        [ 50.]])

For vectors, indexing requires two indices:

>>> print x[0,0], x[1,0]
4.0 5.0

Though np.matrix takes a real matrix form and look pleasing, usually, for most of the cases, arrays are good enough.

>>> import numpy as np
>>> A = np.ones((4,3))
>>> A
array([[ 1.,  1.,  1.],
       [ 1.,  1.,  1.],
       [ 1.,  1.,  1.],
       [ 1.,  1.,  1.]])
>>> np.rank(A)
2

Note that the rank of the array is not the rank of the matrix in linear algebra (dimension of the column space) but the number of subscripts it takes!

Scalars have rank 0:

>>> x = np.array(10)
>>> x
array(10)
>>> np.rank(x)
0

NumPy supports arrays of any dimension such as rank 3 (2x2x2):

>>> A = np.ones((2,2,2))
>>> A
array([[[ 1.,  1.],
        [ 1.,  1.]],

       [[ 1.,  1.],
        [ 1.,  1.]]])
>>> A[1,0,1]
1.0

>>> A = np.array([[1,2],[3,4]])
>>> A
array([[1, 2],
       [3, 4]])
>>> b = np.array([10, 20])
>>> b
array([10, 20])
>>> ans = np.dot(A,b)
>>> ans
array([ 50, 110])

Another example:

>>> A = np.array([[ 1, 2 ,3], [ 4, 5 ,6]])
>>> B = np.array([7,8,9])
>>> 
>>> A
array([[1, 2, 3],
       [4, 5, 6]])
>>> 
>>> B
array([7, 8, 9])
>>> 
>>> A.shape
(2, 3)
>>> B.shape
(3,)
>>> 
>>> A.T
array([[1, 4],
       [2, 5],
       [3, 6]])
>>> 
>>> B.T
array([7, 8, 9])
>>> 
>>> A.dot(B)
array([ 50, 122])
>>> 
>>> np.dot(A,B)
array([ 50, 122])

Now we want to solve Ax = b:

>>> import numpy as np
>>> from numpy.linalg import solve
>>> A = np.array([[1,2],[3,4]])
>>> A
array([[1, 2],
       [3, 4]])
>>> b = np.array([10, 20])
>>> b
>>> x = solve(A,b)
>>> x
array([ 0.,  5.])

>>> import numpy as np
>>> from numpy.linalg import eig
>>> A = np.array([[1,2],[3,4]])
>>> eig(A)
(array([-0.37228132,  5.37228132]), array([[-0.82456484, -0.41597356],
       [ 0.56576746, -0.90937671]]))

The eig returns two tuples: the first one is the eigen values and the second one is a matrix whose columns are the two eigen vectors.

We can unpack the tuples:

>>> eigen_val, eigen_vec = eig(A)
>>> eigen_val
array([-0.37228132,  5.37228132])
>>> eigen_vec
array([[-0.82456484, -0.41597356],
       [ 0.56576746, -0.90937671]])

We want to solve $\int_0^{3}x^4dx = \frac{243}{4}$:

>>> from scipy.integrate import quad
>>> def f(x):
...     return x**4
... 
>>> quad(f, 0., 3.)
(48.599999999999994, 5.39568389967826e-13)

The returned tuple indicates (ans, error estimate).

We can get the same answer if we use lambda instead:

>>> quad(lambda x: x**4, 0, 3)
(48.599999999999994, 5.39568389967826e-13)