Solving overdetermined systems with the QR decomposition

A system of linear equations is considered overdetermined if there are more equations than unknowns. In practice, we have a system Ax=b where A is a m by n matrix and b is a m dimensional vector b but m is greater than n. In this case, the vector b cannot be expressed as a linear combination of the columns of A. Hence, we can't find x so that satisfies the problem Ax=b (except in specific cases) but it is possible to determine x so that Ax is as close to b as possible. So we wish to find x which minimizes the following error

Considering the QR decomposition of A we have that Ax=b becomes

multiplying by Q^T we obtain

and since Q^T is orthogonal (this means that Q^T*Q=I) we have

Now, this is a well defined system, R is an upper triangular matrix and Q^T*b is a vector. More precisely b is the orthogonal projection of b onto the range of A. And,

The function linalg.lstsq() provided by numpy returns the least-squares solution to a linear system equation and is able to solve overdetermined systems. Let's compare the solutions of linalg.lstsq() with the ones computed using the QR decomposition:

from numpy import *

# generating a random overdetermined system
A = random.rand(5,3)
b = random.rand(5,1) 

x_lstsq = linalg.lstsq(A,b)[0] # computing the numpy solution

Q,R = linalg.qr(A) # qr decomposition of A
Qb = dot(Q.T,b) # computing Q^T*b (project b onto the range of A)
x_qr = linalg.solve(R,Qb) # solving R*x = Q^T*b

# comparing the solutions
print 'qr solution'
print x_qr
print 'lstqs solution'
print x_lstsq

This is the output of the script above:

qr solution
[[ 0.08704059]
 [-0.10106932]
 [ 0.56961487]]
lstqs solution
[[ 0.08704059]
 [-0.10106932]
 [ 0.56961487]]

As we can see, the solutions are the same.

Finite differences with Toeplitz matrix

A Toeplitz matrix is a band matrix in which each descending diagonal from left to right is constant. In this post we will see how to approximate the derivative of a function f(x) as matrix-vector products between a Toeplitz matrix and a vector of equally spaced values of f. Let's see how to generate the matrices we need using the function toeplitz(...) provided by numpy:

from numpy import *
from scipy.linalg import toeplitz
import pylab

def forward(size):
 """ returns a toeplitz matrix
   for forward differences
 """
 r = zeros(size)
 c = zeros(size)
 r[0] = -1
 r[size-1] = 1
 c[1] = 1
 return toeplitz(r,c)

def backward(size):
 """ returns a toeplitz matrix
   for backward differences
 """
 r = zeros(size)
 c = zeros(size)
 r[0] = 1
 r[size-1] = -1
 c[1] = -1
 return toeplitz(r,c).T

def central(size):
 """ returns a toeplitz matrix
   for central differences
 """
 r = zeros(size)
 c = zeros(size)
 r[1] = .5
 r[size-1] = -.5
 c[1] = -.5
 c[size-1] = .5
 return toeplitz(r,c).T

# testing the functions printing some 4-by-4 matrices
print 'Forward matrix'
print forward(4)
print 'Backward matrix'
print backward(4)
print 'Central matrix'
print central(4)

The result of the test above is as follows:

Forward matrix
[[-1.  1.  0.  0.]
 [ 0. -1.  1.  0.]
 [ 0.  0. -1.  1.]
 [ 1.  0.  0. -1.]]
Backward matrix
[[ 1.  0.  0. -1.]
 [-1.  1.  0.  0.]
 [ 0. -1.  1.  0.]
 [ 0.  0. -1.  1.]]
Central matrix
[[ 0.   0.5  0.  -0.5]
 [-0.5  0.   0.5  0. ]
 [ 0.  -0.5  0.   0.5]
 [ 0.5  0.  -0.5  0. ]]

We can observe that the matrix-vector product between those matrices and the vector of equally spaced values of f(x) implements, respectively, the following equations:

Forward difference,

Backward difference,

And central difference,

where h is the step size between the samples. Those equations are called Finite Differences and they give us an approximate derivative of f. So, let's approximate some derivatives!

x = linspace(0,10,15)
y = cos(x) # recall, the derivative of cos(x) is sin(x)
# we need the step h to compute f'(x) 
# because the product gives h*f'(x)
h = x[1]-x[2]
# generating the matrices
Tf = forward(15)/h 
Tb = backward(15)/h
Tc = central(15)/h

pylab.subplot(211)
# approximation and plotting
pylab.plot(x,dot(Tf,y),'g',x,dot(Tb,y),'r',x,dot(Tc,y),'m')
pylab.plot(x,sin(x),'b--',linewidth=3)
pylab.axis([0,10,-1,1])

# the same experiment with more samples (h is smaller)
x = linspace(0,10,50)
y = cos(x)
h = x[1]-x[2]
Tf = forward(50)/h
Tb = backward(50)/h
Tc = central(50)/h

pylab.subplot(212)
pylab.plot(x,dot(Tf,y),'g',x,dot(Tb,y),'r',x,dot(Tc,y),'m')
pylab.plot(x,sin(x),'b--',linewidth=3)
pylab.axis([0,10,-1,1])
pylab.legend(['Forward', 'Backward', 'Central', 'True f prime'],loc=4)
pylab.show()

The resulting plot would appear as follows:

As the theory suggests, the approximation is better when h is smaller and the central differences are more accurate (note that, they have a higher order of accuracy respect to the backward and forward ones).

Convolution with numpy

A convolution is a way to combine two sequences, x and w, to get a third sequence, y, that is a filtered version of x. The convolution of the sample x_t is computed as follows:

It is the mean of the weighted summation over a window of length k and w_t are the weights. Usually, the sequence w is generated using a window function. Numpy has a number of window functions already implemented: bartlett, blackman, hamming, hanning and kaiser. So, let's plot some Kaiser windows varying the parameter beta:

import numpy
import pylab

beta = [2,4,16,32]

pylab.figure()
for b in beta:
 w = numpy.kaiser(101,b) 
 pylab.plot(range(len(w)),w,label="beta = "+str(b))
pylab.xlabel('n')
pylab.ylabel('W_K')
pylab.legend()
pylab.show()

The graph would appear as follows:

And now, we can use the function convolve(...) to compute the convolution between a vector x and one of the Kaiser window we have seen above:

def smooth(x,beta):
 """ kaiser window smoothing """
 window_len=11
 # extending the data at beginning and at the end
 # to apply the window at the borders
 s = numpy.r_[x[window_len-1:0:-1],x,x[-1:-window_len:-1]]
 w = numpy.kaiser(window_len,beta)
 y = numpy.convolve(w/w.sum(),s,mode='valid')
 return y[5:len(y)-5]

Let's test it on a random sequence:

# random data generation
y = numpy.random.random(100)*100 
for i in range(100):
 y[i]=y[i]+i**((150-i)/80.0) # modifies the trend

# smoothing the data
pylab.figure(1)
pylab.plot(y,'-k',label="original signal",alpha=.3)
for b in beta:
 yy = smooth(y,b) 
 pylab.plot(yy,label="filtered (beta = "+str(b)+")")
pylab.legend()
pylab.show()

The program would have an output similar to the following:

As we can see, the original sequence have been smoothed by the windows.