QR decomposition
Calculating and using QR decomposition
Intro
QR decomposition (factorization) is decomposition of a matrix into orthogonal (Q) and upper triangular (R) matrices. QR factorization is used in solving linear least square problems and finding eigenvalues. This post shows how QR decomposition is computed and how to use it to solve practical problems.
Gram-Schmidt Process
QR decomposition has following formula:
A = QR, where:
- A is original matrix we want to decompose
- Q is orthogonal matrix
- R is upper triangular matrix
Main goal is rather simple, decompose matrix into matrices Q and R. To find a orthogonal matrix Q, we could used Gram-Schmidt process. This process takes input matrix and makes columns orthogonal (perpendicular to each other). In the following plot we turn vectors a and b to orthogonal vectors v1 and v2.
To make orthogonal vectors, we could use following information. Vector a is first vector and we can match it with v1. v2 must be orthogonal to vector v1, so we should project vector b to v2 to make b orthogonal to a (v1). Projection of vector b to a (v1) is given by:
Information we have is in the following plot:
Now we need to project b to v2 to get orthogonality between a and b. For this we could use e, which is our projection error and could be calculated using following formula:
One thing to note is that e has the same length and the same direction as v2. We can could conclude:
This is basically how Gram-Schmidt process works. For more than two dimensions (assume we have vectors a, b and c) we could calculate process in a following way (source):
What is happening in more dimensions is that we are subtracting from vector directions to previous components (v1, v2 and etc. if we have more dimensions). Now we know how to get matrix Q.
To compute matrix R we can use following changes in QR formula:
A = QR,
QᵗA = QᵗQR (because Q is orthogonal Qᵗ=inverse(Q))
QᵗA = R, thus R = QᵗA
In python QR decomposition using Gram-Schmidt process could be implemented as following:
Example usage and output:
A = np.array([[60, 91, 26], [60, 3, 75], [45, 90, 31]], dtype='float')
Q, R = qr_factorization(A)
Q
#[[ 0.62469505, 0.71080736, -0.63928383],
# [ 0.62469505, 0.02343321, 0.75368929],
# [ 0.46852129, 0.70299629, -0.15254061]]
R
#[[ 96.04686356, 0. , 77.61835966],
# [ 0. , 128.02343535, 0. ],
# [ 0. , 0. , 35.17655816]]np.allclose(A, Q@R)
#True
As usual this is not the only way to calculate QR decomposition. Next let’s look how to use Householder reflections for QR decomposition.
Householder reflections
Gram-Schmidt might be numerically unstable, meaning that small changes input might lead to relatively large changes in output (source). More stable way is to use Householder reflections. Householder projects vector through a “mirror”. We have vector x which we want to make reflect to vector Qx. For reflection we’ll use orthogonal matrix Q.
From the previous we know that x projection to u is:
From the Householder reflection plot we can see that if we subtract twice the component parallel to u from the x, we’ll get Qx (source):
Parallel component happens to be our x projection to u, so we can write:
Because of the associativity of matrix multiplication:
We can write:
Which becomes:
And finally we have:
This is a way how to get matrix Q. Full householder process has following steps:
- select first column (a) of a initial matrix A we want to decompose:
[[ 12, -51, 4],
[ 6, 167, -68],
[ -4, 24, -41]]
- calculate vector u which is perpendicular to mirror, using following formula. Sign(a1) means sign of first element in a and used for avoiding numerical cancellation errors. e1 stands for standard basis vector ([1, 0, 0, …], meaning all except first entry will become zero):
- obtain Householder matrix H for column by using previous formula:
- obtain Q by calculating dot product Q = QH
- update matrix A with householder matrix: A = HA
- take next column from matrix A from the ith row and follow previous steps. Example selection of second column of updated matrix A:
[[-14., -21., 14. ]
[ -0., 173.923, -65.692]
[ 0., 19.385, -42.538]]
- Finally if process has finished, we’ll return Q and updated version A as Q and R matrices.
In python this process is implemented in the following way:
If we use this with function we obtain:
A = np.array([[60, 91, 26], [60, 3, 75], [45, 90, 31]], dtype='float')
Q, R = qr(a)
Q.round(3)
#[[-0.625, 0.355, -0.696],
# [-0.625, -0.762, 0.172],
# [-0.469, 0.542, 0.698]]
R.round(3)
#[[ -96.047, -100.888, -77.618],
# [ -0. , 78.813, -31.083],
# [ 0. , 0. , 16.469]]Q_np, R_np=np.linalg.qr(A)
np.allclose(Q, Q_np)
#True
np.allclose(R, R_np)
#True
Our implementation gives similar results as numpy.
Now we can have look how QR decomposition could be used in practice.
Solving linear equations
We have matrices A, x and B. We know A and B but would like to find x. We can use some algebra (source):
- Ax = B
- QRx = B
- Rx = QᵗB (because Qᵗ = inverse(Q) as Q is orthogonal)
Solving last equation becomes trivial using back substitution (see more information from here), because R is upper triangular matrix. Simple example:
A = np.array([
[2., 1., 1.],
[1., 3., 2.],
[1., 0., 0]])
B = np.array([4., 5., 6.])
#decompose
Q, R = np.linalg.qr(A)
#Q.TB
y = Q.T @ B
#solve for x, code for back_substitution is from here https://gist.github.com/RRisto/2cddc23d829877248915d66ecbe09a39
back_substitution(R, y)
#[ 6., 15., -23.]#using np solver
np.linalg.solve(R, y)
#[ 6., 15., -23.]
QR decomposition makes solving linear equations simple. We can use this decomposition for finding eigenvalues and eigenvectors.
Finding eigenvalues and eigenvectors
Finding eigenvalues and vectors is not very difficult. We can use following procedure (source):
- create matrix X from matrix A and decompose using QR decomposition:
Q, R = np.linalg.qr(X)
- calculate product of Q and previous products (pQ is initialized as diagonal matrix):
pQ = pQ @ Q
- calculate new X using R and Q:
X = R @ Q
- repeat the process n times
Eigenvalues will be diagonal elements of X and eigenvectors are in matrix pQ.
Simple Python implementation for previous algorithm:
Using this function gives us:
A = np.array([[2., 1., 1.],[1., 3., 2.],[1., 0., 0]])
find_eig_qr(A)
#eigenvalues
#[ 3.91222918, 1.28646207, -0.19869124]
#eigenvectors
#[[ 0.51233387, -0.72297297, -0.46349121],
# [ 0.84874276, 0.50856627, 0.14490026],
# [ 0.13095702, -0.46762212, 0.87417379]]#compare with numpy implementation
np.linalg.eig(A)
#eigenvalues
#[ 3.91222918, 1.28646207, -0.19869124]
#eigenvectors
#[[-0.51233387, -0.51118241, -0.17058945],
# [-0.84874276, 0.76210326, -0.48349198],
# [-0.13095702, -0.39735522, 0.85856551]]
We can see that eigenvectors differ from implementation provided by numpy. This a very simple implementation but produces pretty good results.
Conclusion
QR decomposition could be calculated using Gram-Schmidt process and Householder projections, last being more numerically stable. There are more variants to those algorithms, here I presented simplest ways to use them. QR decomposition could be used to solving linear equations and finding eigenvalues and eigenvectors.
References
Finding matrix eigenvectors using QR decomposition, stats.stackexchange.com
Gram-Schmidt Process, ML Wiki
Linear equations with Python: the QR decomposition, b2bvoice.com
QR decomposition, rosettacode.com
The Householder transformation in numerical linear algebra, John Kerl
Why is Householder computationally more stable than modified Gram-Schmidt? math.stackexchange.com