Saturday, March 10, 2012

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.56961487]]
lstqs solution
[[ 0.08704059]
 [ 0.56961487]]
As we can see, the solutions are the same.


  1. Thanks! This is helping me get some image processing homework done.

  2. Succinct. Thank you.

  3. If speed is important to you, you can use the triangular nature of the `R` matrix and solve it much faster using solve_triangular (in scipy.linalg):

    import scipy
    x_qr2 = scipy.linalg.solve_triangular(R, Qb, check_finite=False)

    This is 5.5x (or more) faster than numpy.linalg.solve, for 160 by 160 matrixes at least.

  4. Is there a way to get nonnegative solution for underdertermined system for Ax=b

  5. Is there a way to get nonnegative solution for underdertermined system for Ax=b

    1. Of course, but it is a different kind of optimization problem.