qr_solve, a FORTRAN77 code 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.
The test program also needs the TEST_LLS 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 Fortran90 version and a MATLAB version and an Octave version and a Python version.
BVLS, a FORTRAN77 library which applies least squares methods to solve a linear system for which lower and upper constraints may have been placed on every variable, by Charles Lawson and Richard Hanson.
DQED, a FORTRAN77 library which solves constrained least squares problems, by Richard Hanson and Fred Krogh.
GEQP3, a FORTRAN77 library which contains the portion of the LAPACK library that carries out the QR factorization, with column pivoting, of an M by N rectangular matrix, with N <= M.
LAPACK_EXAMPLES, a FORTRAN77 program which demonstrates the use of the LAPACK linear algebra library.
LAWSON, a FORTRAN77 library which contains routines for solving least squares problems and singular value decompositions (SVD), by Charles Lawson, Richard Hanson.
LINPACK_D, a FORTRAN77 library which solves linear systems using double precision real arithmetic;
LLSQ, a FORTRAN77 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 FORTRAN77 library which solves systems of nonlinear equations, or the least squares minimization of the residual of a set of linear or nonlinear equations.
NMS, a FORTRAN77 library which includes a wide variety of numerical software, including solvers for linear systems of equations, interpolation of data, numerical quadrature, linear least squares data fitting, the solution of nonlinear equations, ordinary differential equations, optimization and nonlinear least squares, simulation and random numbers, trigonometric approximation and Fast Fourier Transforms (FFT).
R8LIB, a FORTRAN77 library which contains many utility routines using double precision real (R8) arithmetic.
TEST_LS, a FORTRAN77 library which implements linear least squares test problems of the form A*x=b.