Solving Overdetermined Systems with the QR Decomposition
Join the DZone community and get the full member experience.
Join For FreeA 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 leastsquares 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:
This is the output of the script above:
As we can see, the solutions are the same.
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 leastsquares 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.
Decomposition (computer science)
Published at DZone with permission of Giuseppe Vettigli, DZone MVB. See the original article here.
Opinions expressed by DZone contributors are their own.
Trending

How to Optimize CPU Performance Through Isolation and System Tuning

HashMap Performance Improvements in Java 8

Which Is Better for IoT: Azure RTOS or FreeRTOS?

How To Scan and Validate Image Uploads in Java
Comments