The LA_LEAST_SQUARES function is used to solve the linear least-squares problem:
Minimizex ||Ax - b||2
where A is a (possibly rank-deficient) n-column by m-row array, b is an m-element input vector, and x is the n-element solution vector. There are three possible cases:
- If m ≥ n and the rank of A is n, then the system is overdetermined and a unique solution may be found, known as the least-squares solution.
- If m < n and the rank of A is m, then the system is under determined and an infinite number of solutions satisfy Ax - b = 0. In this case, the solution is found which minimizes ||x||2, known as the minimum norm solution.
- If A is rank deficient, such that the rank of A is less than MIN(m, n), then the solution is found which minimizes both ||Ax - b||2 and ||x||2, known as the minimum-norm least-squares solution.
The LA_LEAST_SQUARES function may also be used to solve for multiple systems of least squares, with each column of b representing a different set of equations. In this case, the result is a k-by-n array where each of the k columns represents the solution vector for that set of equations.
LA_ LEAST_SQUARES is based on the following LAPACK routines:
sgels, sgelsy, sgelss, sgelsd
dgels, dgelsy, dgelss, dgelsd
cgels, cgelsy, cgelss, cgelsd
zgels, zgelsy, zgelss, zgelsd
Given the following under determined system of equations:
2t + 5u + 3v + 4w = 3
7t + u + 3v + 5w = 1
4t + 3u + 6v + 2w = 6
The following program can be used to find the solution:
; Define the coefficient array:
a = [[2, 5, 3, 4], $
[7, 1, 3, 5], $
[4, 3, 6, 2]]
; Define the right-hand side vector b:
b = [3, 1, 6]
; Find and print the minimum norm solution of a:
x = LA_LEAST_SQUARES(a, b)
PRINT, 'LA_LEAST_SQUARES solution:', x
-0.0376844 0.350628 0.986164 -0.409066
The result is an n-element vector or k-by-n array.
The n-by-m array used in the least-squares system.
An m-element input vector containing the right-hand side of the linear least-squares system, or a k-by-m array, where each of the k columns represents a different least-squares system.
Set this keyword to use double-precision for computations and to return a double-precision (real or complex) result. Set DOUBLE = 0 to use single-precision for computations and to return a single-precision (real or complex) result. The default is /DOUBLE if A is double precision, otherwise the default is DOUBLE = 0.
Set this keyword to indicate which computation method to use. Possible values are:
- METHOD = 0 (the default): Assume that array A has full rank equal to min(m, n). If m ≥ n, find the least-squares solution to the overdetermined system. If m < n, find the minimum norm solution to the under determined system. Both cases use QR or LQ factorization of A.
- METHOD = 1: Assume that array A may be rank deficient; use a complete orthogonal factorization of A to find the minimum norm least-squares solution.
- METHOD = 2: Assume that array A may be rank deficient; use singular value decomposition (SVD) to find the minimum norm least-squares solution.
- METHOD = 3: Assume that array A may be rank deficient; use SVD with a divide-and-conquer algorithm to find the minimum norm least-squares solution. The divide-and-conquer method is faster than regular SVD, but may require more memory.
Set this keyword to a named variable in which to return the effective rank of A. If METHOD = 0 or the array is full rank, then RANK will have the value MIN(m, n).
Set this keyword to the reciprocal condition number used as a cutoff value in determining the effective rank of A. Arrays with condition numbers larger than 1/RCONDITION are assumed to be rank deficient. If RCONDITION is set to zero or omitted, then array A is assumed to be of full rank. This keyword is ignored for METHOD = 0.
If m > n and the rank of A is n (the system is overdetermined), then set this keyword to a named variable in which to return the residual sum-of-squares for Result. If B is an m-element vector then RESIDUAL will be a scalar; if B is a k-by-m array then RESIDUAL will be a k-element vector containing the residual sum-of-squares for each system of equations. If m ≤ n or A is rank deficient (rank < n) then the values in RESIDUAL will be zero.
Set this keyword to a named variable that will contain the status of the computation. Possible values are:
- STATUS = 0: The computation was successful.
- STATUS > 0: For METHOD=2 or METHOD=3, this indicates that the SVD algorithm failed to converge, and STATUS off-diagonal elements of an intermediate bidiagonal form did not converge to zero. For METHOD=0 or METHOD=1 the STATUS will always be zero.
Resources and References
For details see Anderson et al., LAPACK Users' Guide, 3rd ed., SIAM, 1999.