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.
This page has no comments yet. Be the first one!