Solving Linear Systems¶
In [1]:
import numpy as np
import scipy.linalg as la
Solve \(Ax = b\) where \(A\) is full rank square matrix¶
Solve \(Ax = b\) for all the \(b\) vectors
In [2]:
np.random.seed(123)
In [3]:
A = np.random.randint(0, 10, (3,3))
A
Out[3]:
array([[2, 2, 6],
[1, 3, 9],
[6, 1, 0]])
Check if A is invertible.
In [4]:
la.det(A)
Out[4]:
-12.0
In [5]:
b = np.random.randint(0, 10, (3, 10))
b
Out[5]:
array([[1, 9, 0, 0, 9, 3, 4, 0, 0, 4],
[1, 7, 3, 2, 4, 7, 2, 4, 8, 0],
[7, 9, 3, 4, 6, 1, 5, 6, 2, 1]])
Direct matrix inversion¶
Not recommended. Less numerically stable and slower when there are multiple \(b\) to solve for.
In [6]:
%%time
x = np.empty((3, 10))
for i in range(10):
x[:, i] = la.inv(A) @ b[:, i]
CPU times: user 1.39 ms, sys: 0 ns, total: 1.39 ms
Wall time: 1.12 ms
In [7]:
x
Out[7]:
array([[ 0.25 , 3.25 , -1.5 , -1. ,
4.75 , -1.25 , 2. , -2. ,
-4. , 3. ],
[ 5.5 , -10.5 , 12. , 10. ,
-22.5 , 8.5 , -7. , 18. ,
26. , -17. ],
[ -1.75 , 3.91666667, -3.5 , -3. ,
7.41666667, -1.91666667, 2.33333333, -5.33333333,
-7.33333333, 5.33333333]])
Using solve¶
Recommended method.
In [24]:
%%time
la.solve(A, b)
CPU times: user 741 µs, sys: 1.13 ms, total: 1.88 ms
Wall time: 779 µs
Out[24]:
array([[ 0.25 , 3.25 , -1.5 , -1. ,
4.75 , -1.25 , 2. , -2. ,
-4. , 3. ],
[ 5.5 , -10.5 , 12. , 10. ,
-22.5 , 8.5 , -7. , 18. ,
26. , -17. ],
[ -1.75 , 3.91666667, -3.5 , -3. ,
7.41666667, -1.91666667, 2.33333333, -5.33333333,
-7.33333333, 5.33333333]])
Using PLU decomposition¶
The permuation matrix \(P\) indicates the row switches. Solving requires
- forward substitution \(Ly = P^Tb\)
- backward substitution \(Ux = y\)
In [21]:
%%time
P, L, U = la.lu(A)
y = la.solve_triangular(L, P.T@b, lower=True)
x = la.solve_triangular(U, y, lower=False)
CPU times: user 842 µs, sys: 1.3 ms, total: 2.15 ms
Wall time: 795 µs
In [22]:
x
Out[22]:
array([[ 0.25 , 3.25 , -1.5 , -1. ,
4.75 , -1.25 , 2. , -2. ,
-4. , 3. ],
[ 5.5 , -10.5 , 12. , 10. ,
-22.5 , 8.5 , -7. , 18. ,
26. , -17. ],
[ -1.75 , 3.91666667, -3.5 , -3. ,
7.41666667, -1.91666667, 2.33333333, -5.33333333,
-7.33333333, 5.33333333]])
Solving least squares equations¶
We solve \(X\beta = y\) where the dimensions of \(X\) are such that \(m > n\). Here we interpreet the solution as finding the best coefficinets for fittting the eequation
\(y = \beta_0 + \beta_1 x_1 + \beta_2 x_2\)
Note that we are looking for the projection of \(y\) onto the column space (image) of the column vectors of \(X\).
In [26]:
X = np.c_[np.ones((10, 1)), np.random.randint(0, 10, (10, 2))]
In [27]:
X
Out[27]:
array([[1., 8., 3.],
[1., 5., 0.],
[1., 2., 6.],
[1., 2., 4.],
[1., 4., 6.],
[1., 3., 0.],
[1., 6., 4.],
[1., 7., 6.],
[1., 7., 1.],
[1., 5., 7.]])
In [28]:
y = np.random.randint(0, 10, (10, 1))
In [29]:
y
Out[29]:
array([[9],
[2],
[4],
[8],
[1],
[2],
[1],
[1],
[3],
[5]])
Using lstsq method¶
Recommended.
In [31]:
%%time
β, res, rank, s = la.lstsq(X, y)
CPU times: user 643 µs, sys: 0 ns, total: 643 µs
Wall time: 623 µs
In [32]:
β
Out[32]:
array([[ 3.51230193],
[-0.02659447],
[ 0.05892189]])
Using the normal equations¶
In [35]:
la.solve(X.T@X, X.T@y)
Out[35]:
array([[ 3.51230193],
[-0.02659447],
[ 0.05892189]])
Using optimization¶
In [56]:
def res(β, X, y):
β = β.reshape((-1,1))
return (y - X@β).T @ (y - X@β)
In [57]:
β0 = np.zeros(3)
In [58]:
import scipy.optimize as opt
In [59]:
opt.fmin(res, β0, args=(X, y))
Optimization terminated successfully.
Current function value: 76.138865
Iterations: 169
Function evaluations: 304
Out[59]:
array([ 3.51234729, -0.02660293, 0.05892217])
In [ ]: