*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_lstsqThis 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.

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

ReplyDeleteI'm always happy to help you David :)

DeleteSuccinct. Thank you.

ReplyDeleteIf 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):

ReplyDeleteimport 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.

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

ReplyDeleteIs there a way to get nonnegative solution for underdertermined system for Ax=b

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

Delete