In [1]:

    

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 process_time
%matplotlib inline


    

In [74]:

    

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':
        U = sla.cholesky(A.T.dot(A))
        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 [42]:

    

A, b = example(4, 4)


    

In [80]:

    

m = 100
times = np.zeros(m)
for n in range(m, m):
    A, b = example(m, n)
    print(n)
    t_start = process_time()
    x = solve_lsq(A, b, method='normal')
    times[n] = process_time() - t_start


    

In [111]:

    

plt.plot(times)


    

    
Out[111]:





[<matplotlib.lines.Line2D at 0x7f28cb6ac9e8>]






    

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]:





False

In [81]:

    

U, s, V = sla.svd(A, full_matrices=False)
V.dot(np.diag(1 / s).dot(U.T.dot(b)))


    

    
Out[81]:





array([ 0.87888801,  0.73738623,  1.2775084 , -0.94054901,  1.30438628,
        1.13090557,  0.91060691,  1.13696161,  0.25475408])

In [82]:

    

b


    

    
Out[82]:





array([ 1.        ,  0.62108984,  0.64609802,  0.7500381 ,  0.9       ,
        1.125     ,  1.49992379,  2.23865695,  4.03128129,  9.        ])

In [ ]: