In [106]:
import numpy as np
import matplotlib.pyplot as plt
import scipy.linalg as sla
import numpy.linalg as LA
import pandas as pd
from time import clock
%matplotlib inline
In [46]:
def example(m, n):
A = np.vander(np.linspace(-1., 1., m), n)
b = A.dot(np.ones(n))
return A, b
def solve_lsq(A, b, method='qr'):
if method=='normal':
try:
U = sla.cholesky(A.T.dot(A))
except:
ATA = A.T.dot(A)
ATA[np.abs(ATA) < 10 ** -14] = 0.
U = sla.cholesky(ATA)
y = sla.solve_triangular(U.T, A.T.dot(b), lower=True)
return sla.solve_triangular(U, y)
if method=='qr':
Q, R = sla.qr(A, mode='economic')
return sla.solve_triangular(R, Q.T.dot(b)[:R.shape[0]])
if method=='svd':
U, s, V = sla.svd(A, full_matrices=False)
return V.dot(np.diag(1 / s).dot(U.T.dot(b)))
def test_lsq():
A = np.random.random((4, 4))
b = np.random.random(4)
x = solve_lsq(A, b, 'normal')
print(A.T.dot(b) - A.T.dot(A).dot(x))
In [50]:
A, b = example(100, 27)
x = solve_lsq(A, b, method='normal')
np.isclose(A.dot(x), b)
Out[50]:
In [119]:
m = 100
times = np.zeros(m)
for mode in ['qr', 'svd']:
for n in range(1, m):
A, b = example(m, n)
t_start = clock()
x = solve_lsq(A, b, method=mode)
times[n] = clock() - t_start
plt.plot(times, label=mode)
times = np.zeros(m)
for n in range(1, 26):
A, b = example(m, n)
t_start = clock()
x = solve_lsq(A, b, method='normal')
times[n] = clock() - t_start
plt.plot(times[:26])
Out[119]:
In [109]:
for n in range(1, 26):
A, b = example(m, n)
Out[109]:
In [79]:
m = 10
n = 9
A, b = example(m, n)
t_start = process_time()
x = solve_lsq(A, b, method='svd')
np.allclose(A.dot(x), b)
Out[79]:
In [81]:
U, s, V = sla.svd(A, full_matrices=False)
V.dot(np.diag(1 / s).dot(U.T.dot(b)))
Out[81]:
In [82]:
b
Out[82]:
In [26]:
1e-15
Out[26]:
In [99]:
A, b = example(m, 10)
U,s,V = LA.svd(A)
In [84]:
U.shape, s.shape, V.shape
Out[84]:
In [85]:
N = np.diag(1 /s)
W = np.zeros((10, 90))
v = np.hstack((W, N))
In [90]:
dot(,U.T)
In [100]:
c = np.dot(U.T,b) # c = U^t*b
w = LA.solve(np.diag(s),c) # w = V^t*c
xSVD = np.dot(V.T,w)
In [89]:
A.shape, x.shape
Out[89]:
In [93]:
x
Out[93]:
In [ ]: