LA_LEAST_SQUARES

The LA_LEAST_SQUARES function is used to solve the linear least-squares problem:

Minimize_x||Ax - b||₂

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||₂, 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||₂ and ||x||₂, 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:

Output Type	LAPACK Routines
Float	sgels, sgelsy, sgelss, sgelsd
Double	dgels, dgelsy, dgelss, dgelsd
Complex	cgels, cgelsy, cgelss, cgelsd
Double complex	zgels, zgelsy, zgelss, zgelsd

Examples

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

IDL prints:

LA_LEAST_SQUARES solution:

-0.0376844     0.350628     0.986164    -0.409066

Syntax

Result = LA_LEAST_SQUARES( A, B [, /DOUBLE] [, METHOD=value] [, RANK=variable] [, RCONDITION=value] [, RESIDUAL=variable] [, STATUS=variable] )

Return Value

The result is an n-element vector or k-by-n array.

Arguments

A

The n-by-m array used in the least-squares system.

B

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.

Keywords

DOUBLE

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.

METHOD

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.

RANK

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).

RCONDITION

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.

RESIDUAL

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.

STATUS

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.

Version History

5.6	Introduced

Resources and References

For details see Anderson et al., LAPACK Users' Guide, 3rd ed., SIAM, 1999.