Thursday, May 5, 2011

QR decomposition with numpy

We will see how to compute the QR decomposition of a matrix A and how to use Q and R to solve the linear equation system Ax=b using the from described here.
from numpy import *

A = floor(random.rand(4,4)*20-10) # random matrix

Q,R = linalg.qr(A) # qr decomposition of A

b = floor(random.rand(4,1)*20-10)
# solve Ax = b using the standard numpy function
x = linalg.solve(A,b)

# solve Ax = b using Q and R
y = dot(Q.T,b)
xQR = linalg.solve(R,y) 

print "\nSolution compared"
print x.T,'Ax=b'
print xQR.T,'Rx=y'
And now we can compare the solutions obtained
Solution compared
[[ 0.69207502  0.05565638  0.68965517 -0.2183908 ]] Ax=b
[[ 0.69207502  0.05565638  0.68965517 -0.2183908 ]] Rx=y

2 comments:

  1. When using QR decomposition in Numpy, the first basis vector that it chooses can sometimes affect the numerical accuracy of the solution. See this post for an example where the L1-norm of the difference between the QR decomp solution and the "exact" solution was not zero:

    http://bpchesney.org/?p=657

    ReplyDelete