Documentation |
Least squares solution of linear equations
This functionality does not run in MATLAB.
numeric::leastSquares(A, B, <mode>, <method>, options)
numeric::leastSquares(A, B) computes a matrix X that solves the linear matrix equation A X = B in the least squares sense: the columns X_{j} of X minimize where the B_{j} are the columns of B.
For a given vector B, a vector X minimizes if and only if X is a solution of the "normal equations" A^{H} A X = A^{H} B, where A^{H} is the Hermitian transpose of the m×n matrix A. The solution is unique if rank(A) = n.
numeric::leastSquares allows to solve several least squares problems simultaneously by combining several ‘right hand sides' B_{j} columnwise to a matrix B.
If no return type is specified via the option ReturnType = d, the domain type of the return data depends on the type of the input matrix A:
The special solution X as well as the kernel of an array A are returned as arrays.
The special solution and the kernel of an hfarray of domain type DOM_HFARRAY are returned as hfarrays.
For a dense matrix A of type Dom::DenseMatrix(), both the special solution X as well as the kernel are returned as matrices of type Dom::DenseMatrix() over the ring of MuPAD^{®} expressions.
For all other matrices of category Cat::Matrix, both the special solution X as well as the kernel are returned as matrices of type Dom::Matrix() over the ring of MuPAD expressions. This includes input matrices A of type Dom::Matrix(...), Dom::SquareMatrix(...), Dom::MatrixGroup(...) etc.
Without Symbolic, the input data are converted to floating-point numbers. The matrix A must not contain non-convertible parameters, unless Symbolic is used. If such objects are found, numeric::leastSquares automatically switches to its symbolic mode, issuing a warning. This warning may be suppressed via NoWarning.
Symbolic parameters in B are accepted without warning. However, HardwareFloats cannot be used if there are any symbolic parameters in A or B.
If A^{H} A has a non-trivial kernel, the least squares solution X is not unique. The return value X is a special solution of the equation A^{H} A X = A^{H} B. With the SVD method, X is the special solution with columns of minimal Euclidean length.
Note: The result computed with HardwareFloats may differ from the solution computed with SoftwareFloats or Symbolic! In particular, this is the case for systems with a non-trivial kernel. Further, the The results computed with QRD and SVD may differ. |
The kernel is computed only in the symbolic mode (option Symbolic). All floating-point methods return the value NIL for the kernel.
With Symbolic, the n×d matrix KernelBasis is the most general solution of A^{H} A X = 0. Its columns span the d-dimensional kernel of A^{H} A.
If the kernel is 0-dimensional, the return value of KernelBasis is the integer 0. If KernelBasis is returned as an array, the dimension d of the kernel is d = op(KernelBasis, [0, 3, 2]]). If KernelBasis is returned as a matrix of type Dom::Matrix() or Dom::DenseMatrix(), the dimension d of the kernel is d = KernelBasis::dom::matdim(KernelBasis)[2].
Note: Without the option Symbolic, the implemented algorithms take care of numerical stabilization. With Symbolic, exact data are assumed. The least squares solutions is computed via numeric::matlinsolve( A^{H} A, A^{H} B , Symbolic). The symbolic strategy tries do maximize speed and does not take care of numerical stabilization! Do not use Symbolic for systems involving floating-point entries! In particluar, due to round-off, it may happen that no solution of A^{H} A X = A^{H} B is found. In such a case, [FAIL, NIL, NIL] is returned. Cf. Example 4. |
All entries of A and B must be arithmetical expressions.
Note: Apart from matrices of type Dom::Matrix(...), Cat::Matrix objects A from matrix domains such as Dom::DenseMatrix(...) or Dom::SquareMatrix(...) are internally converted to arrays over expressions via expr(A). Note that the option Symbolic should be used if the entries cannot be converted to numerical expressions. The same holds true for matrices B passed as Cat::Matrix objects. |
The function is sensitive to the environment variable DIGITS, which determines the numerical working precision.
We consider a matrix A of rank 1:
A := array(1..3, 1..2, [[1, 2], [1, 2], [1, 2]]): B := [3, 4, 5]:
The normal equations have a 1-parameter set of of solutions:
[X, KernelBasis, Res] := numeric::leastSquares(A, B, Symbolic)
The numerical method QRD produces a special solution:
[X, KernelBasis, Res] := numeric::leastSquares(A, B, QRD)
The numerical method SVD produces a solution X of minimal norm:
[X, KernelBasis, Res] := numeric::leastSquares(A, B, SVD)
delete A, B, X, KernelBasis, Res:
We consider an ill-conditioned least squares problem. By construction, the following overdetermined system has an exact solution X = [1, 2, …, n]:
m := 10: n := 8: A := array(1..m, 1..n, [[1/(i + j + 100) $ j=1..n] $ i=1..m]): B := array(1..m, [_plus(A[i,j]*j $ j=1..n) $ i=1..m]): numeric::leastSquares(A, B, Symbolic)
The coefficient matrix A is rather ill-conditioned:
singvals := numeric::singularvalues(A): conditionOfA := max(op(singvals))/min(op(singvals))
Consequently, round-off has a drastic effect in a numerical approximation. The methods yield results of different quality:
numeric::leastSquares(A, B, QRD)
numeric::leastSquares(A, B, SVD)
delete m, n, A, B, singvals, conditionOfA:
This example involves a symbolic parameter c in the matrix A. The option Symbolic must be used:
A:= matrix([[c, 2], [1/3, 2/3], [1/7, 2/7]]): B:= [1, 2, 3]: numeric::leastSquares(A, B, Symbolic)
normal(%)
delete A, B:
Floating point entries may cause problems in conjunction with the option Symbolic, because the computation is not stabilized numerically in the symbolic node. The following matrix A has rank 2:
A := matrix([[1, 30], [10.0^(-15), 31*10.0^(-15)]]):
However, due to round-off, the ‘normal matrix' A^{H} A has rank 1. No solution is found with Symbolic:
A::dom::transpose(A) * A
numeric::leastSquares(A, [31, 32*10^5], Symbolic)
No such problem arises in the numerical schemes. Note, however, that the numerical methods yield different results in this extremely ill-conditioned problem:
numeric::leastSquares(A, [31, 32*10^5], QRD)
numeric::leastSquares(A, [31, 32*10^5], SVD)
delete A:
A |
An m×n matrix of domain type DOM_ARRAY or of category Cat::Matrix |
B |
An m×p matrix of domain type DOM_ARRAY or of category Cat::Matrix. Column vectors B may also be represented by a 1-dimensional array(1..m, [B_{1}, B_{2}, …] ) or by a list [B_{1}, B_{2}, …]. |
mode |
One of the flags Hard, HardwareFloats, Soft, SoftwareFloats, or Symbolic |
method |
One of the flags QRD, SVD |
Hard, HardwareFloats, Soft, SoftwareFloats |
With Hard (or HardwareFloats), computations are done using fast hardware float arithmetic from within a MuPAD session. Hard and HardwareFloats are equivalent. With this option, the input data are converted to hardware floats and processed by compiled C code. The result is reconverted to MuPAD floats and returned to the MuPAD session. With Soft (or SoftwareFloats) computations are dome using software float arithmetic provided by the MuPAD kernel. Soft and SoftwareFloats are equivalent. SoftwareFloats is used by default if the current value of DIGITS is larger than 15 and the input matrix A is not of domain type DOM_HFARRAY. Compared to the SoftwareFloats used by the MuPAD kernel, the computation with HardwareFloats may be many times faster. Note, however, that the precision of hardware arithmetic is limited to about 15 digits. Further, the size of floating-point numbers may not be larger than approximately 10^{308} and not smaller than approximately 10^{- 308}. If no HardwareFloats or SoftwareFloats are requested explicitly, the following strategy is used: If the current value of DIGITS is smaller than 16 or if the matrix A is a hardware float array of domain type DOM_HFARRAY, then hardware arithmetic is tried. If this is successful, the result is returned. If the result cannot be computed with hardware floats, software arithmetic by the MuPAD kernel is tried. If the current value of DIGITS is larger than 15 and the input matrix A is not of domain type DOM_HFARRAY, or if one of the options Soft, SoftwareFloats or Symbolic is specified, MuPAD computes the result with its software arithmetic without trying to use hardware floats first. There may be several reasons for hardware arithmetic to fail:
If neither HardwareFloats nor SoftwareFloats is specified, the user is not informed whether hardware floats or software floats are used. If HardwareFloats are specified but fail due to one of the reasons above, a warning is issued that the (much slower) software floating-point arithmetic of the MuPAD kernel is used. Note that HardwareFloats can only be used if all input data can be converted to floating-point numbers. The trailing digits in floating-point results computed with HardwareFloats and SoftwareFloats may differ. | |
Symbolic |
Prevents the conversion of the input data to floating-point numbers. Exact arithmetic is used. This option overrides HardwareFloats and SoftwareFloats. This option must be used, if the matrix A contains symbolic parameters that cannot be converted to floating-point numbers. The normal equations A^{H} A X = A^{H} B are passed to numeric::matlinsolve with the option Symbolic. If the least squares problem does not have a unique solution, a special solution X is returned together with the kernel of A^{H} A. Cf. Example 1.
| |
QRD |
Use a QR decomposition. All entries of A must be convertible to floating-point values. This is the default method. The matrix A must not contain symbolic parameters that cannot be converted to floating point numbers. If such objects are found, then numeric::leastSquares automatically switches to its Symbolic mode, issuing a warning. The computation proceeds with exact arithmetic, using the input data without floating-point conversions. The warning may be suppressed by the option NoWarning. Symbolic parameters in B are accepted without warning. They are processed by the floating-point algorithm. Numerical expressions such as etc. are accepted and converted to floats. If the least squares problem does not have a unique solution, only a special solution is returned. The kernel is not computed: it is returned as NIL. The method QRD provides a numerically stable way of solving the normal equations A^{H} A X = A^{H} B by a QR decomposition. In extremely ill-conditioned situations, it may be worthwhile to consider the slower, yet more stable method SVD. The conditioning is given by the ratio of the largest singular value of A divided by the smallest singular value of A. If this value is large, the problem is ill-conditioned. Cf. Example 2. | |
SVD |
Use a singular value decomposition. All entries of A must be convertible to floating-point values. The matrix A must not contain symbolic parameters that cannot be converted to floating point numbers. If such objects are found, then numeric::leastSquares automatically switches to its symbolic mode, issuing a warning. The computation proceeds with exact arithmetic, using the input data without floating point conversions. The warning may be suppressed by the option NoWarning. Symbolic parameters in B are accepted without warning. They are processed by the floating-point algorithm. Numerical expressions such as etc. are accepted and converted to floats. If the least squares problem does not have a unique solution, the columns X_{j} of the solution X have a minimal Euclidean length . The kernel is not computed: it is returned as NIL. A singular value decomposition A = U D V^{H} is used to solve the normal equations in the form D^{2} V^{H} X = D U^{H} B. For small or zero singular values d_{j} in D = diag(d_{1}, d_{2}, …), the corresponding components of V^{H} x are set to zero. Usually, the numerical method SVD is slower than the default method QRD. However, in ill-conditioned situations, it is numerically more stable. The conditioning is given by the ratio of the largest singular value of A divided by the smallest singular value of A. If this value is large, the problem is ill-conditioned. | |
NoWarning |
Suppresses warnings If symbolic coefficients are found in A, numeric::leastSquares automatically switches to the Symbolic mode with a warning. With this option, this warning is suppressed. | |
ReturnType |
Option, specified as ReturnType = d Return the (special) solution and the kernel as matrices of domain type d. The following return types d are available: DOM_ARRAY, or DOM_HFARRAY, Dom::Matrix(), or Dom::DenseMatrix(). |
A list [X, KernelBasis, Residues] is returned.
The (special) least squares solution X is an n×p matrix.
With Symbolic, KernelBasis is an n×d matrix (d is the dimension of the kernel of A^{H} A). Its columns span the kernel of A^{H} A. If the kernel is trivial, KernelBasis is the integer 0.
Without Symbolic, the kernel is not computed. The value NIL is returned for the KernelBasis.
The list of arithmetical expressions Residues consists of the minimized least squares deviations corresponding to the columns of X and B.