qr_solve, a Python 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 test program 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.
r8lib, a Python code which contains many utility routines using double precision real (R8) arithmetic.
test_lls, a Python code which implements linear least squares test problems of the form A*x=b.