numeric::leastSquares

Least squares solution of linear equations

Use only in the MuPAD Notebook Interface.

This functionality does not run in MATLAB.

Syntax

```numeric::leastSquares(`A`, `B`, <`mode`>, <`method`>, `options`)
```

Description

`numeric::leastSquares(A, B)` computes a matrix X that solves the linear matrix equation AX = B in the least squares sense: the columns Xj of `X` minimize where the Bj 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" AHAX = AHB, where AH 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' Bj 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 AHA has a non-trivial kernel, the least squares solution X is not unique. The return value `X` is a special solution of the equation AHAX = AHB. 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 AHAX = 0. Its columns span the d-dimensional kernel of AHA.

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( AH A, AH 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 AH A X = AH 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.

Environment Interactions

The function is sensitive to the environment variable `DIGITS`, which determines the numerical working precision.

Examples

Example 1

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:`

Example 2

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:`

Example 3

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:`

Example 4

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' AHA 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:`

Parameters

 `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, [B1, B2, …] )``` or by a list ```[B1, B2, …]```. `mode` One of the flags `Hard`, `HardwareFloats`, `Soft`, `SoftwareFloats`, or `Symbolic` `method` One of the flags `QRD`, `SVD`

Options

`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 10308 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:

• The current value of `DIGITS` is larger than 15.

• The data contains symbolic objects.

• The data contains numbers larger than 10308 or smaller than 10- 308 that cannot be represented by hardware floats.

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 AHAX = AHB 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 AHA. Cf. Example 1.

 Note:   This option should not be used for systems with floating-point coefficients! Numerical instabilities may occur in floating-point operations. Further, if the rank of A is not maximal, then `numeric::leastSquares` may fail to find a solution due to numerical round-off. In such a case, ```[FAIL, NIL, NIL]``` is returned. Cf. Example 4.

`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 AHAX = AHB 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 Xj of the solution X have a minimal Euclidean length .

The kernel is not computed: it is returned as `NIL`.

A singular value decomposition A = UDVH is used to solve the normal equations in the form D2VHX = DUHB. For small or zero singular values dj in D = diag(d1, d2, …), the corresponding components of VHx 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()`.

Return Values

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 AHA). Its columns span the kernel of AHA. 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`.