Stability of least squares algorithms

Stability of least squares algorithms#

We consider the problem of fitting a polynomial to a given function on \(m\) uniformly spaced points.

\[ p(t) = \exp(\sin(4t), \qquad t \in [0,1] \]

\(p\) is of degree \(n-1\).

\[ p(t_i) = b_i = \exp(\sin(4t_i), \qquad 0 \le i \le m-1 \]

This leads to over-determined system

\[\begin{split} A x = b \qquad \textrm{where} \qquad t = \begin{bmatrix} t_0 \\ t_1 \\ \vdots \\ t_{m-1} \end{bmatrix} \in R^m \qquad A = [t^0, t^1, \ldots, t^{n-1}] \in R^{m \times n} \end{split}\]

The elements of \(x\) are the coefficients of the monomials \(t^i\) in \(p(t)\).

from numpy import linspace,empty,exp,sin,arcsin,column_stack
from numpy.linalg import lstsq,norm,qr,svd
from scipy.linalg import solve,solve_triangular
from qr import mgs,house

Generate matrix \(A\) and right hand side \(b\)

m, n = 100, 15
t = linspace(0,1,m)
A = empty((m,n))
for i in range(n):
    A[:,i] = t**i
b = exp(sin(4*t))
b = b/2006.787453080206

Because of the normalization, the last element of \(x\) which is coefficient of \(t^{n-1}\) should be equal to 1.

Solve the problem using Python functions and examine the condition number \(\kappa\), the angle \(\theta\) between \(b\) and range(A) and \(\eta\).

x,_,_,s = lstsq(A,b,rcond=None); y = A @ x
kappa = s.max() / s.min()
print('kappa = %10.4e' % kappa)
theta = arcsin(norm(b-y)/norm(b))
print('theta = %10.4e' % theta)
eta   = s.max() * norm(x) / norm(y)
print('eta   = %10.4e' % eta)

kappa = 2.2718e+10
theta = 3.7461e-06
eta   = 2.1036e+05

Using reduced QR, the least squares solution is obtained from

\[ A = Q R, \qquad R x = Q^* b, \qquad R = \textrm{upper triangular} \]

In the augmented algorithm, form the augmented matrix \([A \ b]\) and find its QR decomposition

\[ [A \ b] = Q R \]

Then

\[ Q^* b = \textrm{last column of R} \]

Using numpy.linalg.qr#

# numpy.linalg.qr
Q,R = qr(A)
x = solve_triangular(R, Q.T @ b)
print('\nnumpy.linalg.qr')
print('x_15, x_15-1 =',x[-1], x[-1]-1)

numpy.linalg.qr
x_15, x_15-1 = 1.0000000926361714 9.263617140042868e-08

# numpy.linalg.qr augmented
Q2,R2 = qr(column_stack((A,b)))
Qb = R2[0:n,n]   # Last column of R2
R  = R2[0:n,0:n]
x = solve_triangular(R,Qb)
print('\nnumpy.linalg.qr augmented')
print('x_15, x_15-1 =',x[-1], x[-1]-1)

numpy.linalg.qr augmented
x_15, x_15-1 = 1.0000000926395132 9.26395131717328e-08

Using our own Householder triangularization#

# Householder triangularization
Q,R = house(A,mode='reduced')
x = solve_triangular(R, Q.T @ b)
print('\nHouseholder triangularization')
print('x_15, x_15-1 =',x[-1], x[-1]-1)

Householder triangularization
x_15, x_15-1 = 0.9999997622620326 -2.3773796742343478e-07

numpy.linalg.qr is doing some column pivoting so it has an order of magnitude smaller error.

Using modified Gram-Schmidt#

# Modified GS
Q, R = mgs(A)
x = solve_triangular(R, Q.T@b)
print('\nModified GS')
print('x_15, x_15-1 =',x[-1], x[-1]-1)

Modified GS
x_15, x_15-1 = 0.943853356170208 -0.05614664382979195

# Modified GS, augmented
Q2,R2 = mgs(column_stack((A,b)))
Qb = R2[0:n,n]   # Last column of R2
R  = R2[0:n,0:n]
x = solve_triangular(R,Qb)
print('\nModified GS augmented')
print('x_15, x_15-1 =',x[-1], x[-1]-1)

Modified GS augmented
x_15, x_15-1 = 1.0000000397100135 3.971001349967196e-08

The augmented mgs has similar error to numpy.linalg.qr which is based on householder.

Using normal equations#

# Normal equation
x = solve(A.T@A, A.T@b, assume_a='sym')
print('\nNormal equations')
print('x_15, x_15-1 =',x[-1], x[-1]-1)

Normal equations
x_15, x_15-1 = -0.8970336811736241 -1.8970336811736241

/var/folders/tc/17nf7v3j53jgg4db9vxl40g00000gn/T/ipykernel_91921/1749421503.py:2: LinAlgWarning: Ill-conditioned matrix (rcond=3.88679e-19): result may not be accurate.
  x = solve(A.T@A, A.T@b, assume_a='sym')

Using SVD#

Using reduced SVD

\[ A = U \Sigma V^*, \qquad x = V \Sigma^{-1} U^* b \]
# SVD
U,S,VT = svd(A,full_matrices=False) # reduced SVD
x = VT.T @ ((U.T@b)/S)
print('\nSVD')
print('x_15,x_15-1=',x[-1], x[-1]-1)

SVD
x_15,x_15-1= 1.000000092629737 9.262973699186716e-08

Exercise#

Instead of using monomials as a basis for polynomials, use Legendre polynomials and solve the least squares problem using the different methods.

Contents