QR_SOLVE is a C++ library which computes a least squares solution of a linear system of the form A*x=b.
There are many possible cases that can arise with the matrix A. Formally, we distinguish the cases M < N, M = N, and M > N, and we expect trouble whenever M is not equal to N. Trouble may also arise when M = N but the matrix is singular. Even if the matrix is, mathematically speaking, non-singular, it may be so close to singularity that an accurate solution is difficult to achieve.
When M > N, we are placing more conditions than we have degrees of freedom, so we suppose that such a linear system cannot be solved. However, it is possible that the extra conditions are illusory, being constructed from linear combinations of a fundamental set of N conditions. Thus, a system that we typically call "overdetermined" can have a solution in the ordinary sense, that satisfies all the conditions, as long as the right hand side is "consistent". Another way of saying this is that the system is solvable if the right hand side lies in the column space of A...although that simply says that it is a linear combination of the columns of A, which just says A*x=b.
If A does not have full column rank, however, then even if the right hand side lies in the column space of A, there will be more than one linear combination of columns that produce b. Thus, the equations are consistent, the system is solvable, but not uniquely so.
If M < N, then we are placing fewer conditions than we have degrees of freedom. As long as the right hand side lies in the column space of A, we can guarantee that there will be multiple solutions.
Thus, the question of a "solution" to the problem A*x=b is complicated enough that it seems to defy a common algorithmic approach. Nonetheless, there are some sensible, robust procedures for producing an answer that corresponds to the classical solution, or solves the overdetermined problem when the right hand side is consistent. This is the linear least squares solution, which finds a vector x which minimizes the Euclidean norm of the residual: ||r|| = ||A*x-b||. This solution is produced by computing the QR factorization of the matrix A
When there are multiple solutions to the problem, the QR approach used here produces a solution. A more satisfactory approach, using the pseudoinverse, will produce a solution x which satisfies the additional constraint that it has the minimum norm of the entire family of solutions. That pseudoinverse approach is not implemented in this library. The singular value decomposition (SVD) can also produce this minimal solution.
For comparison, a solver that applies the normal equations is included. This approach requires M >= N, and that A have full column rank. It constructs and solves the NxN system A'*A*x=A'*b. This system has a condition number that is the square of the original system, so it also suffers from a significant loss in accuracy.
We also include an SVD solver, which uses the pseudoinverse approach. First compute A = U * S * V', where U and V are orthogonal, and S is MxN diagonal, then to solve A*x=b write x = V * S^ * U' * b, where S^ is the matrix formed by transposing S and then replacing each nonzero diagonal element s by 1/s. (However, very small elements should probably be zeroed instead of inverted.) This procedure will also produce a vector x which minimizes the Euclidean norm. However, it has one feature that the QR solver does not: in cases where the solution x is not unique (because A does not have full column rank) the SVD solver returns the unique vector x of minimum Euclidean norm.
QR_SOLVE needs the R8LIB library. The test program also needs the TEST_LS library.
The computer code and data files described and made available on this web page are distributed under the GNU LGPL license.
QR_SOLVE is available in a C version and a C++ version and a FORTRAN77 version and a FORTRAN90 version.
LINPACK_D, a C++ library which solves linear systems using double precision real arithmetic;
LLSQ, a C++ library which solves the simple linear least squares problem of finding the formula of a straight line y=a*x+b which minimizes the root-mean-square error to a set of N data points.
MINPACK, a C++ library which solves systems of nonlinear equations, or the least squares minimization of the residual of a set of linear or nonlinear equations.
R8LIB, a C++ library which contains many utility routines using double precision real (R8) arithmetic.
TEST_LS, a C++ library which implements linear least squares test problems of the form A*x=b.
You can go up one level to the C++ source codes.