Here I investigate whether QR factorization can have superior accuracy to LU factorization with partial pivoting. QR Factorization has many stability advantages over LU factorization, which suggests that a well-implemented QR factorization could produce a more accurate solution to a system of linear equations \( Ax=b \). I construct a numerical experiment based on a family of matrices which represents a wide range of troublesome numerical behavior, such as lack of symmetry, severe non-normality, and poor conditioning and for a random sample of matrices in this family I solve the same system of linear equations with QR factorization and with LU factorization with partial pivoting. I do this experiment in C++ using LAPACK.
QR factorization has favorable worst-case stability properties compared to LU factorization. The big reason is that it can be implemented with stable orthogonal transformations like Householder reflectors, whose use does not depend on finding a sufficiently large pivot. LU factorization, on the other hand, uses Gaussian elimination; partial pivoting works very well in practice but admits unfavorable worst cases. This does not automatically imply a smaller forward error for QR, which is what I test here.
QR factorization is not typically used to solve a square, full-rank system because LU requires fewer floating-point operations. Both can spend most of their arithmetic in level 3 BLAS, so it is still interesting to ask whether QR’s stronger worst-case stability properties produce more accurate solutions.
I construct a family of matrices and sample from it randomly. For each sample I solve the same system of linear equations with QR factorization and again with LU factorization. I then compute the relative error from the true solution. I also compute condition numbers of each matrix to give an impression of the relative difficulty of solving systems with them.
The family of matrices is generated with the following Python code
def make_matrix(nx,ny,coeffs):
m=nx*ny
A=np.zeros((m,m))
def id(ix,iy):
return ix + nx*iy
for ix in range(0,nx):
for iy in range(0,ny):
A[id(ix,iy),id(ix,iy)]=coeffs[0]
if ix>0:
A[id(ix-1,iy),id(ix,iy)]=coeffs[1]
if ix<nx-1:
A[id(ix+1,iy),id(ix,iy)]=coeffs[2]
if iy>0:
A[id(ix,iy-1),id(ix,iy)]=coeffs[3]
if iy<ny-1:
A[id(ix,iy+1),id(ix,iy)]=coeffs[4]
return A
the family of matrices defined by the function make_matrix is fully defined by the parameters nx,ny,coeffs.
Ordinarily I would prefer to do everything in SciPy, but at the time I wrote this experiment I did not find a convenient SciPy interface to multiply by the implicitly represented matrix \( Q ^ T \) from a QR factorization without forming \( Q \) explicitly. Since I am specifically interested in the errors resulting from a quality implementation of a QR-based solve I think it’s appropriate that application of \(Q\) be done not by forming it explicitly but rather directly from its representation using Householder reflectors - which requires less arithmetic and storage. As a result I did the experiment in C++ to have easier access to LAPACK, but for simplicity I outline the approach taken in python-like pseudocode
import numpy as np
nsamples=10000
nx=32
ny=32
m=nx*ny
soln=np.random.rand(m)
for i in range(0,nsamples):
#Generate random matrix with uniform coefficients between -4 and 4
coeffs=np.random.uniform(-4.0,4.0,5)
A=make_matrix(nx,ny,coeffs)
#Manufacture RHS which results in desired solution
b=A@soln
#Solve with QR factorization
xqr=qrsolve(A,b)
#Solve with LU factorization
xlu=lusolve(A,b)
#get relative error for QR solve
qr_err = np.amax( np.abs(xqr - soln)/np.abs(soln) )
#get relative error for LU solve
lu_err = np.amax( np.abs(xlu - soln)/np.abs(soln) )
#Compute other quantities of interest below (conditioning, loss of symmetry, etc)
#Omitted
As mentioned above I had to implement this with raw LAPACK. There are many routines for computing orthogonal factorizations and doing various other calculations with LAPACK. I list all the routines I used below so the interested reader may easily find them by internet search:
I have made the code available at this repository. Note: I have not fully documented everything there.
To generate these results I computed \(10000\) matrices and computed conditioning, loss of symmetry, non-normality, and finally the errors from a QR solve and a LU solve for the same system of linear equations.
The purpose of the other computations is to establish that this family of matrices encompasses a wide range of potentially problematic numerical behavior for which we might expect a more stable calculation like QR factorization to perform better than a less stable one like LU factorization.
For this calculation I computed the condition number of each matrix sample and plotted its logarithm as a histogram.
For this calculation I computed the loss of symmetry \( \frac{| A ^ T - A |}{|A + A^T |} \) and plotted this as a histogram.
For this calculation I computed the loss of normality \( \frac{| AA ^ T - A^T A |}{|AA^T + A^TA |} \) and plotted this as a histogram.
Finally I generated a random solution \( x \) and manufactured a right-hand side \( b \) such that \( Ax = b \). I then tried to reconstruct \( x\) using first QR factorization and then LU factorization. I plot their errors below on the same scales.
and then the same errors but instead expressed as a ratio
Despite favorable worst-case stability properties, QR factorization did not outperform LU factorization in forward error when solving the same linear system on this family of matrices. Perhaps there are other families of matrices for which the improved stability of the QR factorization could materially improve on the accuracy of the LU factorization, but it did not for this family despite having a fairly wide range of potentially problematic numerical properties.
This was just a curiosity of mine, and even if QR factorization did turn out to be more accurate, that alone would not be a very compelling reason to use it instead of LU factorization for solving a system of linear equations. The reason is that one can apply iterative refinement for a relatively small extra cost and often recover additional digits of accuracy, subject mainly to the conditioning of the problem.
The following books provide complete and total coverage of QR factorization, LU factorization, their mathematical properties and proper implementation. I referenced them liberally while writing this blog post and code.