Basic Linear Algebra with numpy

Linear Algebra

  • Scipy and numpy have powerful linear algebra functionality.
  • Vector operations
    • basic operations (multiplication by scalar, addition, etc.)
    • dot product
    • cross product
  • Matrix operations
    • basic operations (multiplication by scalar, mult. by vector, mult. by matrix, addition, etc.)
    • inverse
    • transpose
    • determinant
    • trace

Linear Algebra Operations

Linear Albebra Operations

  • contained in scipy.linalg or numpy.linalg module
  • Solving linear systems: A x = b with A as a matrix and x, b as vectors.
  • Finding eigenvalues, eigenvectors.
  • Singular value decomposition (SVD).
  • Various matrix factorizations (LU, Cholesky, etc.)

Making matrices from 2-D arrays

  • use the mat or matrix functions.

  • Can take an array, list of lists, or string to construct:

    >>> A = mat([[3, 1, 4], [1, 5, 9], [2, 6, 5]])
    >>> A = mat('3 1 4; 1 5 9; 2 6 5') # equivalent
    >>> print A
    [[3 1 4]
     [1 5 9]
     [2 6 5]]
    >>> print type(A), A.dtype, A.shape
    <class 'numpy.core.defmatrix.matrix'> int32 (3, 3)
    

Special Matrix Properties

Objects of the numpy.core.defmatrix.matrix class have several additional built-in properties.

Property Gives
.T transpose
.H Hermitian transpose (transpose with complex conjugate)
.I matrix inverse
.A matrix as a basic ndarray
.A1 matrix as a one-dimensional ndarray

Examples

>>> A
matrix([[3, 1, 4],
        [1, 5, 9],
        [2, 6, 5]])
>>> A.T # transpose
matrix([[3, 1, 2],
        [1, 5, 6],
        [4, 9, 5]])
>>> A.H # Hermitian transpose (same for real matrix)
matrix([[3, 1, 2],
        [1, 5, 6],
        [4, 9, 5]])
>>> A.I # matrix inverse
matrix([[ 0.32222222, -0.21111111,  0.12222222],
        [-0.14444444, -0.07777778,  0.25555556],
        [ 0.04444444,  0.17777778, -0.15555556]])
>>> A.A # as an array
array([[3, 1, 4],
       [1, 5, 9],
       [2, 6, 5]])
>>> A.A1 # as a 1-dimensional array
array([3, 1, 4, 1, 5, 9, 2, 6, 5])

Basic Matrix Operations

>>> A
matrix([[3, 1, 4],
        [1, 5, 9],
        [2, 6, 5]])

Matrix addition:

>>> A + eye(3) # eye(3) = 3x3 identity
matrix([[4, 1, 4],
        [1, 6, 9],
        [2, 6, 6]])
>>> set_printoptions(suppress=True) # print very small numbers as 0.

Matrix multiplication:

>>> A*A.I
matrix([[ 1.,  0., -0.],
        [ 0.,  1., -0.],
        [-0.,  0.,  1.]])

Matrix powers:

>>> A**3 # A*A*A
matrix([[ 168,  424,  565],
        [ 346,  990, 1294],
        [ 302,  854, 1081]])

Trancendental functions (from the Taylor series):

>>> set_printoptions(precision=3)
>>> linalg.expm(A) # matrix exponential
array([[  33159.108,   93457.367,  120428.208],
       [  75981.418,  214248.526,  276060.62 ],
       [  64583.091,  182098.641,  234637.133]])
>>> linalg.sinm(A) # matrix sine
array([[ 0.473, -0.203,  0.251],
       [-0.13 ,  0.092,  0.543],
       [ 0.107,  0.37 ,  0.006]])

Matrix-vector multiplication

  • Use a 1-column or 1-row matrix as a vector:

    >>> b = mat('1; 2; 3')
    
    >>> b # column vector
    matrix([[1],
            [2],
            [3]])
    
    >>> A*b # matrix multiply by vector
    matrix([[17],
            [38],
            [29]])
    
    >>> b.T * A # gives row vector
    matrix([[17, 41, 55]])
    
    >>> (b.T * A).A1 # turn back to 1-D array
    array([17, 41, 55])
    

Solving a set of linear equations

Problem:
solve A x = b for x.
Solution #1 (only recommended for small matrices):

A-1A x = A-1 b

I x = A-1 b

x = A-1 b

Example:

>>> set_printoptions(precision = 6)
>>> A
matrix([[3, 1, 4],
        [1, 5, 9],
        [2, 6, 5]])
>>> b # RHS
matrix([[1],
        [2],
        [3]])
>>> x = A.I * b
>>> x # solution
matrix([[ 0.266667],
        [ 0.466667],
        [-0.066667]])
>>> A*x - b # solution verification
matrix([[ 0.],
        [-0.],
        [ 0.]])

Better way to solve linear equations

  • Use linalg.solve(A, b) instead of multiplication by inverse, especially for large matrices.

  • Faster and more robust than using inverse.

  • Example:

    >>> linalg.solve(A, b)
    array([[ 0.266667],
           [ 0.466667],
           [-0.066667]])