matrices in python



In [9]:

    
import numpy as np
import scipy as sp
import scipy.linalg as linalg

ndarray vs matrix

ndarray is recommended

creating



In [16]:

    
A = sp.matrix([[1,2,3],[3,1,2],[4,5,7]])
a = np.array([[1,2,3],[3,1,2],[4,5,7]])
A, a









    Out[16]:





(matrix([[1, 2, 3],
         [3, 1, 2],
         [4, 5, 7]]), array([[1, 2, 3],
        [3, 1, 2],
        [4, 5, 7]]))

transpose



In [13]:

    
A.T, a.T









    Out[13]:





(matrix([[1, 3, 4],
         [2, 1, 5],
         [3, 2, 7]]), array([[1, 3, 4],
        [2, 1, 5],
        [3, 2, 7]]))

inverse



In [14]:

    
A.I, linalg.inv(a)









    Out[14]:





(matrix([[-0.75,  0.25,  0.25],
         [-3.25, -1.25,  1.75],
         [ 2.75,  0.75, -1.25]]), array([[-0.75,  0.25,  0.25],
        [-3.25, -1.25,  1.75],
        [ 2.75,  0.75, -1.25]]))

matrix multiplication



In [24]:

    
B = np.matrix([2,2,3])
b = np.array([2,2,3])
A*B.T, a.dot(b)









    Out[24]:





(matrix([[15],
         [14],
         [39]]), array([15, 14, 39]))

as vertical vector



In [40]:

    
B.T, b[:,np.newaxis]









    Out[40]:





(matrix([[2],
         [2],
         [3]]), array([[2],
        [2],
        [3]]))

elementwise multiplication



In [27]:

    
np.multiply(A,B), a*b









    Out[27]:





(matrix([[ 2,  4,  9],
         [ 6,  2,  6],
         [ 8, 10, 21]]), array([[ 2,  4,  9],
        [ 6,  2,  6],
        [ 8, 10, 21]]))



In [33]:

    
np.multiply(A,B.T), a*b[:,np.newaxis]









    Out[33]:





(matrix([[ 2,  4,  6],
         [ 6,  2,  4],
         [12, 15, 21]]), array([[ 2,  4,  6],
        [ 6,  2,  4],
        [12, 15, 21]]))

solving linear system

using inverse matrix is not recommended

slow



In [54]:

    
a = np.random.randint(low=-10,high=10,size=(100,100))
b = np.random.randint(low=-10,high=10,size=(100))
%timeit linalg.inv(a).dot(b)









    



1000 loops, best of 3: 698 µs per loop

fast



In [55]:

    
%timeit linalg.solve(a,b)









    



1000 loops, best of 3: 324 µs per loop



In [77]:

    
a = np.random.randint(low=0, high=10, size=(2,2))
b = np.random.randint(low=0, high=10, size=2)
a,b









    Out[77]:





(array([[8, 2],
        [0, 6]]), array([8, 8]))



In [78]:

    
linalg.solve(a,b)









    Out[78]:





array([ 0.66666667,  1.33333333])

eigenvalues and eigenvector

If there exists a nonzero vector $x$ such that

$Ax = \lambda x$

then $\lambda$ is called $A$'s egenvalues and $x$ is $A$'s eigenvector.

how to find $\lambda$ and $x$:

$Ax = \lambda Ix$ add identidy matrix

$Ax - \lambda Ix = 0$ subtract one side

$(A - \lambda I)x = 0$ factorize

This have a solution when det$(A - \lambda I) = 0$. In case of a 2x2 matrix this gives the equation

$\begin{vmatrix}

\end{vmatrix}



In [79]:

    
help(linalg.eigh)









    



Help on function eigh in module scipy.linalg.decomp:

eigh(a, b=None, lower=True, eigvals_only=False, overwrite_a=False, overwrite_b=False, turbo=True, eigvals=None, type=1, check_finite=True)
    Solve an ordinary or generalized eigenvalue problem for a complex
    Hermitian or real symmetric matrix.
    
    Find eigenvalues w and optionally eigenvectors v of matrix `a`, where
    `b` is positive definite::
    
                      a v[:,i] = w[i] b v[:,i]
        v[i,:].conj() a v[:,i] = w[i]
        v[i,:].conj() b v[:,i] = 1
    
    Parameters
    ----------
    a : (M, M) array_like
        A complex Hermitian or real symmetric matrix whose eigenvalues and
        eigenvectors will be computed.
    b : (M, M) array_like, optional
        A complex Hermitian or real symmetric definite positive matrix in.
        If omitted, identity matrix is assumed.
    lower : bool, optional
        Whether the pertinent array data is taken from the lower or upper
        triangle of `a`. (Default: lower)
    eigvals_only : bool, optional
        Whether to calculate only eigenvalues and no eigenvectors.
        (Default: both are calculated)
    turbo : bool, optional
        Use divide and conquer algorithm (faster but expensive in memory,
        only for generalized eigenvalue problem and if eigvals=None)
    eigvals : tuple (lo, hi), optional
        Indexes of the smallest and largest (in ascending order) eigenvalues
        and corresponding eigenvectors to be returned: 0 <= lo <= hi <= M-1.
        If omitted, all eigenvalues and eigenvectors are returned.
    type : int, optional
        Specifies the problem type to be solved:
    
           type = 1: a   v[:,i] = w[i] b v[:,i]
    
           type = 2: a b v[:,i] = w[i]   v[:,i]
    
           type = 3: b a v[:,i] = w[i]   v[:,i]
    overwrite_a : bool, optional
        Whether to overwrite data in `a` (may improve performance)
    overwrite_b : bool, optional
        Whether to overwrite data in `b` (may improve performance)
    check_finite : boolean, optional
        Whether to check that the input matrices contain only finite numbers.
        Disabling may give a performance gain, but may result in problems
        (crashes, non-termination) if the inputs do contain infinities or NaNs.
    
    Returns
    -------
    w : (N,) float ndarray
        The N (1<=N<=M) selected eigenvalues, in ascending order, each
        repeated according to its multiplicity.
    v : (M, N) complex ndarray
        (if eigvals_only == False)
    
        The normalized selected eigenvector corresponding to the
        eigenvalue w[i] is the column v[:,i].
    
        Normalization:
    
            type 1 and 3: v.conj() a      v  = w
    
            type 2: inv(v).conj() a  inv(v) = w
    
            type = 1 or 2: v.conj() b      v  = I
    
            type = 3: v.conj() inv(b) v  = I
    
    Raises
    ------
    LinAlgError :
        If eigenvalue computation does not converge,
        an error occurred, or b matrix is not definite positive. Note that
        if input matrices are not symmetric or hermitian, no error is reported
        but results will be wrong.
    
    See Also
    --------
    eig : eigenvalues and right eigenvectors for non-symmetric arrays



In [80]:

    
values, vector = linalg.eigh(a)
values.shape, vector.shape









    Out[80]:





((2,), (2, 2))



In [81]:

    
print(values)









    



[ 6.  8.]



In [82]:

    
print(vector)









    



[[ 0. -1.]
 [-1. -0.]]