Solving Overdetermined Systems with the QR Decomposition
Join the DZone community and get the full member experience.Join For Free
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) # 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.
Published at DZone with permission of Giuseppe Vettigli, DZone MVB. See the original article here.
Opinions expressed by DZone contributors are their own.