qr_solve, a FORTRAN90 code which computes a linear least squares (LLS) solution of a system 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 code. 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 code needs the R8LIB code. The test code also needs the TEST_LLS code.
The computer code and data files described and made available on this web page are distributed under the MIT license
qr_solve is available in a C version and a C++ version and a FORTRAN90 version and a MATLAB version and a Python version.
BVLS, a FORTRAN90 code which applies linear least squares (LLS) 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.
LAWSON, a FORTRAN90 code which contains routines for solving least squares problems and singular value decompositions (SVD), by Charles Lawson, Richard Hanson.
LINPACK_D, a FORTRAN90 code which solves linear systems using double precision real arithmetic;
LLSQ, a FORTRAN90 code 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.
NMS, a FORTRAN90 code 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 FORTRAN90 code which contains many utility routines using double precision real (R8) arithmetic.
test_lls, a FORTRAN90 code which implements linear least squares (LLS) test problems which are given A and b, and seek a vector x which minimizes the L2 norm of the error ||Ax-b||.