How to find QR decomposition of a rectangular matrix in overdetermined linear system solution?

Question

While trying to find cell-centered gradients in finite volume method computation of incompressible fluid flow I get over-determined linear system. This is a well known "cell based least-square" gradient reconstruction.

I know in theory that by using least squares we search for a solution which minimizes error vector in 2-norm, and that we can solve the problem by calculating 'pseudo-inverse'

$x=(A^TA)^{-1}A^Tb$

but I read that it is not typically how the problem is solved. Rather $QR$ decomposition is performed on original system matrix, and eventually we get

$Rx=Q^Tb$

What I'm looking for is a fast algorithm for $QR$ decomposition of rectangular matrices that I can implement in Fortran, or any other suggestion related to particular problem that I mentioned.

System matrix is usually 6x3.

Thanks!

Answer 1

If you are looking for an efficient and numerically stable algorithm to solve the least squares problem you can use Householder reflectors. Golub and van Loan "Matrix Computations'' gives a very nice treatment on Householder reflectors (in fact one of the dedications in the book is for Alston Scott Householder). We'll outline the basic idea behind the linear least squares problem and how we can use Householder reflectors to solve it. Suppose we have an overdetermined matrix $A\in\mathbb{R}^{n\times k}$ with $n\geq k$ and full rank $$ A=\left(\!\!{\begin{array}{cccc} * & * & * & * \\ * & * & * & * \\ * & * & * & * \\ * & * & * & * \\ * & * & * & * \\ * & * & * & * \\ \end{array}}\!\!\right). $$ Now let us take the vector $v$ to be the $i$th column of the matrix $A$ running from the diagonal elements $i$ to $n$. Then we calculate the $i$th Householder reflector $H_i$ to zero out all the entries below $A_{ii}$. We define $\gamma = sign(e_1^Tv)\|v\|_2$, where $e_1^T = [1\;0\;\cdots\;0]$, $$ x = v+\gamma e_1\qquad and \qquad \alpha = -\frac{2}{\|x\|_2^2}. $$ The $sign$ function on $\gamma$ is chosen for numerical stability. Then the reflector $H_i$ is given by $$ H_i = I+\alpha xx^T. $$ Note that the reflector $H_i$ will be an orthogonal matrix which will zero every element in the vector $v$ except the first value, i.e., $$ H_iv = (I+\alpha xx^T)v = -\gamma e_1. $$ This can be done in a computationally efficient way because of the structure of the reflector. We do not need to form the orthogonal matrix for each $H_i$ instead we can compute its action as $$ H_iv = (I+\alpha xx^T)v = v+\alpha x(x^Tv) = -\gamma e_1, $$ which is a single dot product and then a scaled vector addition. So we have turned an $\mathcal{O}(n^3)$ operation into an $\mathcal{O}(n)$ operation. We can then create a sequence of $k$ reflectors to drive the rectangular matrix $A$ to an upper trapezoidal form. For example, $$ H_4H_3H_2H_1A = \left(\!\!{\begin{array}{cccc} * & * & * & * \\ 0 & * & * & * \\ 0 & 0 & * & * \\ 0 & 0 & 0 & * \\ 0 & 0 & 0 & * \\ 0 & 0 & 0 & 0 \\ \end{array}}\!\!\right) = \left(\!{\begin{array}{c} R \\ 0 \end{array}}\!\right). $$ The product of these $k$ Householder reflectors is an $n\times n$ orthogonal matrix which we denote $H=H_kH_{k-1}\cdots H_1$. We wish to find the minimizing solution $x_{min}$ to the overdetermined system $Ax=b$. We transform the system using the matrix $H$ into the system $HAx = Hb$ or $$ \left(\!{\begin{array}{c} R \\ 0 \end{array}}\!\right)x = \left(\!{\begin{array}{c} c \\ d \end{array}}\!\right). $$ Then the minimal solution with respect to the $2$-norm is $x_{min} = R^{-1}c$ which can be solved efficiently using a back substitution method.

There are more details on the geometry of the least squares problem (given in Golub and van Loan) and how the Householder reflectors project onto two spaces, the range of $A$ and the perp space of $A$. One can prove interesting results on the residual of linear least squares problems, it is related to the value $d$ above. But the original post seemed more an implementation question versus a theoretical linear algebra question so we'll omit such discussion.

Answer 2

If you're using Fortran and don't know much about linear algebra, I'd suggest you use LAPACK.

The relevant routine in there for solving an over-determined system of equations using a QR-decomposition is DGELS.

Since your matrix is quite small, you may get significantly more preformance with a custom QR-decomposition, but LAPACK is probably the more robust solution, and definitely the fastest in terms of development time.

Answer 3

Modified Gram-Schmidt is fastest for your problem in terms of implementation time and performance. To avoid loss of orthogonality you must append the matrix with the right hand side. See Golub, van Loan for details.