Question

In general, the linear algebra functions can be found in scipy.linalg. You can also get access to BLAS and LAPACK function via scipy.linagl.blas and scipy.linalg.lapack.

import scipy.linalg as la

A = np.array([[1,2],[3,4]])
b = np.array([1,4])
print(A)
print(b)

[[1 2]
 [3 4]]
[1 4]

# Matrix operations
import numpy as np
import scipy.linalg as la
from functools import reduce

A = np.array([[1,2],[3,4]])
print(np.dot(A, A))
print(A)
print(la.inv(A))
print(A.T)

[[ 7 10]
 [15 22]]
[[1 2]
 [3 4]]
[[-2.   1. ]
 [ 1.5 -0.5]]
[[1 3]
 [2 4]]

x = la.solve(A, b) # do not use x = dot(inv(A), b) as it is inefficient and numerically unstable
print(x)
print(np.dot(A, x) - b)

[ 2.  -0.5]
[ 0.  0.]

Matrix decompositions

A = np.floor(npr.normal(100, 15, (6, 10)))
print(A)

[[  94.   82.  125.  108.  105.   88.   99.   82.   97.  112.]
 [  83.  124.   67.  103.   73.  111.  125.   81.  122.   62.]
 [  93.   84.  107.  107.   80.   85.   96.   89.   85.  102.]
 [ 116.  116.   64.   98.   82.   98.  121.   70.  122.   98.]
 [ 118.  108.  103.  102.   68.   98.   88.   78.  103.   95.]
 [ 112.  115.   74.   80.  106.  104.  114.  105.   80.   99.]]

P, L, U = la.lu(A)
print(np.dot(P.T, A))
print
print(np.dot(L, U))

[[ 118.  108.  103.  102.   68.   98.   88.   78.  103.   95.]
 [  83.  124.   67.  103.   73.  111.  125.   81.  122.   62.]
 [  94.   82.  125.  108.  105.   88.   99.   82.   97.  112.]
 [ 116.  116.   64.   98.   82.   98.  121.   70.  122.   98.]
 [ 112.  115.   74.   80.  106.  104.  114.  105.   80.   99.]
 [  93.   84.  107.  107.   80.   85.   96.   89.   85.  102.]]

[[ 118.  108.  103.  102.   68.   98.   88.   78.  103.   95.]
 [  83.  124.   67.  103.   73.  111.  125.   81.  122.   62.]
 [  94.   82.  125.  108.  105.   88.   99.   82.   97.  112.]
 [ 116.  116.   64.   98.   82.   98.  121.   70.  122.   98.]
 [ 112.  115.   74.   80.  106.  104.  114.  105.   80.   99.]
 [  93.   84.  107.  107.   80.   85.   96.   89.   85.  102.]]

Q, R = la.qr(A)
print(A)
print
print(np.dot(Q, R))

[[  94.   82.  125.  108.  105.   88.   99.   82.   97.  112.]
 [  83.  124.   67.  103.   73.  111.  125.   81.  122.   62.]
 [  93.   84.  107.  107.   80.   85.   96.   89.   85.  102.]
 [ 116.  116.   64.   98.   82.   98.  121.   70.  122.   98.]
 [ 118.  108.  103.  102.   68.   98.   88.   78.  103.   95.]
 [ 112.  115.   74.   80.  106.  104.  114.  105.   80.   99.]]

[[  94.   82.  125.  108.  105.   88.   99.   82.   97.  112.]
 [  83.  124.   67.  103.   73.  111.  125.   81.  122.   62.]
 [  93.   84.  107.  107.   80.   85.   96.   89.   85.  102.]
 [ 116.  116.   64.   98.   82.   98.  121.   70.  122.   98.]
 [ 118.  108.  103.  102.   68.   98.   88.   78.  103.   95.]
 [ 112.  115.   74.   80.  106.  104.  114.  105.   80.   99.]]

U, s, V = la.svd(A)
m, n = A.shape
S =  np.zeros((m, n))
for i, _s in enumerate(s):
    S[i,i] = _s
print(reduce(np.dot, [U, S, V]))

[[  94.   82.  125.  108.  105.   88.   99.   82.   97.  112.]
 [  83.  124.   67.  103.   73.  111.  125.   81.  122.   62.]
 [  93.   84.  107.  107.   80.   85.   96.   89.   85.  102.]
 [ 116.  116.   64.   98.   82.   98.  121.   70.  122.   98.]
 [ 118.  108.  103.  102.   68.   98.   88.   78.  103.   95.]
 [ 112.  115.   74.   80.  106.  104.  114.  105.   80.   99.]]

B = np.cov(A)
print(B)

[[ 187.7333 -182.4667   94.9333 -105.4444    1.2    -137.2   ]
 [-182.4667  609.6556  -83.3111  371.0556   90.8778   70.5667]
 [  94.9333  -83.3111   97.2889  -48.8889   45.0222  -79.8   ]
 [-105.4444  371.0556  -48.8889  438.5     145.5     109.0556]
 [   1.2      90.8778   45.0222  145.5     215.4333  -39.7667]
 [-137.2      70.5667  -79.8     109.0556  -39.7667  234.1   ]]

u, V = la.eig(B)
print(np.dot(B, V))
print
print(np.real(np.dot(V, np.diag(u))))

[[-280.8911  157.1032   12.1003  -60.7161    8.8142   -1.5134]
 [ 739.1179   34.4268    3.8974    4.3778   14.9092 -122.8749]
 [-134.1449  128.3162  -11.0569   -6.6382   37.3675   13.4467]
 [ 598.7992   77.4348   -5.3372  -52.7843  -14.996    94.553 ]
 [ 170.8339  193.7335    5.8732   67.6135    1.1042   90.1451]
 [ 199.7105 -218.1547    6.1467   -5.6295   26.3372  101.0444]]

[[-280.8911  157.1032   12.1003  -60.7161    8.8142   -1.5134]
 [ 739.1179   34.4268    3.8974    4.3778   14.9092 -122.8749]
 [-134.1449  128.3162  -11.0569   -6.6382   37.3675   13.4467]
 [ 598.7992   77.4348   -5.3372  -52.7843  -14.996    94.553 ]
 [ 170.8339  193.7335    5.8732   67.6135    1.1042   90.1451]
 [ 199.7105 -218.1547    6.1467   -5.6295   26.3372  101.0444]]

C = la.cholesky(B)
print(np.dot(C.T, C))
print
print(B)

[[ 187.7333 -182.4667   94.9333 -105.4444    1.2    -137.2   ]
 [-182.4667  609.6556  -83.3111  371.0556   90.8778   70.5667]
 [  94.9333  -83.3111   97.2889  -48.8889   45.0222  -79.8   ]
 [-105.4444  371.0556  -48.8889  438.5     145.5     109.0556]
 [   1.2      90.8778   45.0222  145.5     215.4333  -39.7667]
 [-137.2      70.5667  -79.8     109.0556  -39.7667  234.1   ]]

[[ 187.7333 -182.4667   94.9333 -105.4444    1.2    -137.2   ]
 [-182.4667  609.6556  -83.3111  371.0556   90.8778   70.5667]
 [  94.9333  -83.3111   97.2889  -48.8889   45.0222  -79.8   ]
 [-105.4444  371.0556  -48.8889  438.5     145.5     109.0556]
 [   1.2      90.8778   45.0222  145.5     215.4333  -39.7667]
 [-137.2      70.5667  -79.8     109.0556  -39.7667  234.1   ]]

Finding the covariance matrix

np.random.seed(123)
x = np.random.multivariate_normal([10,10], np.array([[3,1],[1,5]]), 10)
# create a zero mean array
u = x - x.mean(0)
cov = np.dot(u.T, u)/(10-1)
print cov, '\n'
print np.cov(x.T)

[[ 5.1286  3.0701]
 [ 3.0701  9.0755]]

[[ 5.1286  3.0701]
 [ 3.0701  9.0755]]

Least squares solution

Suppose we want to solve a system of noisy linear equations

\[\begin{split}y_1 = b_0 x_1 + b_1 \\ y_2 = b_0 x_2 + b_1 \\ y_3 = b_0 x_2 + b_1 \\ y_4 = b_0 x_4 + b_1 \\\end{split}\]

Since the system is noisy (implies full rank) and overdetermined, we cannot find an exact solution. Instead, we will look for the least squares solution. First we can rewrrite in matrix notation \(Y = AB\), treating \(b_1\) as the coefficient of \(x^0 = 1\):

\[\begin{split}\left( \begin{array}{c} y_1 \\ y_2 \\ y_3 \\ y_4 \end{array} \right) = \left( \begin{array}{cc} x_1 & 1 \\ x_2 & 1 \\ x_3 & 1 \\ x_4 & 1 \end{array} \right) \left( \begin{array}{cc} b_0 & b_1 \end{array} \right)\end{split}\]

The solution of this (i.e. the \(B\) matrix) is solved by multipling the psudoinverse of \(A\) (the Vandermonde matrix) with \(Y\)

\[(A^\text{T}A)^{-1}A^\text{T} Y\]

Note that higher order polynomials have the same structure and can be solved in the same way

\[\begin{split}\left( \begin{array}{c} y_1 \\ y_2 \\ y_3 \\ y_4 \end{array} \right) = \left( \begin{array}{ccc} x_1^2 & x_1 & 1 \\ x_2^2 & x_2 & 1 \\ x_3^2 & x_3 & 1 \\ x_4^2 & x_4 & 1 \end{array} \right) \left( \begin{array}{ccc} b_0 & b_1 & b_2 \end{array} \right)\end{split}\]

# Set up a system of 11 linear equations
x = np.linspace(1,2,11)
y = 6*x - 2 + npr.normal(0, 0.3, len(x))

# Form the VanderMonde matrix
A = np.vstack([x, np.ones(len(x))]).T

# The linear algebra librayr has a lstsq() function
# that will do the above calculaitons for us

b, resids, rank, sv = la.lstsq(A, y)

# Check against pseudoinverse and the normal equation
print("lstsq solution".ljust(30), b)
print("pseudoinverse solution".ljust(30), np.dot(la.pinv(A), y))
print("normal euqation solution".ljust(30), np.dot(np.dot(la.inv(np.dot(A.T, A)), A.T), y))

# Now plot the solution
xi = np.linspace(1,2,11)
yi = b[0]*xi + b[1]

plt.plot(x, y, 'o')
plt.plot(xi, yi, 'r-');

('lstsq solution                ', array([ 5.5899, -1.4177]))
('pseudoinverse solution        ', array([ 5.5899, -1.4177]))
('normal euqation solution      ', array([ 5.5899, -1.4177]))

# As advertised, this works for finding coeefficeints of a polynomial too

x = np.linspace(0,2,11)
y = 6*x*x + .5*x + 2 + npr.normal(0, 0.6, len(x))
plt.plot(x, y, 'o')
A = np.vstack([x*x, x, np.ones(len(x))]).T
b = la.lstsq(A, y)[0]

xi = np.linspace(0,2,11)
yi = b[0]*xi*xi + b[1]*xi + b[2]
plt.plot(xi, yi, 'r-');

# It is important to understand what is going on,
# but we don't have to work so hard to fit a polynomial

b = np.random.randint(0, 10, 6)
x = np.linspace(0, 1, 25)
y = np.poly1d(b)(x)
y += np.random.normal(0, 5, y.shape)

p = np.poly1d(np.polyfit(x, y, len(b)-1))
plt.plot(x, y, 'bo')
plt.plot(x, p(x), 'r-')
list(zip(b, p.coeffs))

[(6, -250.9964),
 (7, 819.7606),
 (1, -909.5724),
 (5, 449.7862),
 (7, -91.2660),
 (9, 15.5274)]