# numeric::matlinsolve

Solve a linear matrix equation

### Use only in the MuPAD Notebook Interface.

This functionality does not run in MATLAB.

## Syntax

```numeric::matlinsolve(`A`, `B`, `options`)
```

## Description

`numeric::matlinsolve(A, B)` returns the matrix solution X of the matrix equation AX = B together with the kernel of the matrix A.

`numeric::matlinsolve` is a fast numerical linear solver for both sparse and dense systems. It is also a recommended solver for linear systems with exact or symbolic coefficients (use option `Symbolic`).

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 X as well as the kernel of an hfarray A are returned as hfarrays.

• For a dense matrixA of type `Dom::DenseMatrix()`, both the special solution X as well as the kernel of A 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 of A 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`, exact numerical input data such as `PI + sqrt(2)`, `sin(3)` etc. are converted to floating-point numbers. Floating point arithmetic is used. Its precision is given the environment variable `DIGITS`. If symbolic data are found that cannot be converted to floating-point numbers, `numeric::matlinsolve` automatically switches to its symbolic mode, issuing a warning. This warning may be suppressed via `NoWarning`.

With `Symbolic`, symbolic parameters in the coefficient matrix `A` as well as in the right hand side `B` are accepted and processed without a warning.

With `SofwareFloats`, the right hand side `B` may contain symbolic parameters that cannot be converted to floating-point numbers. All entries of the coefficient matrix `A`, however, must be convertible to floating-point numbers.

With `HardwareFloats`, neither A nor B must contain symbolic parameters that cannot be converted to floating-point numbers.

`X` is a special solution of the equation AX = B. If A has a non-trivial kernel, the solution X is not unique.

 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.

Cf. Example 9.

The n×d matrix `KernelBasis` is the most general solution of AX = 0. Its columns span the d-dimensional kernel of A.

Thus, the kernel of A may be computed via `numeric::matlinsolve(A, [0, ..., 0])[2]`.

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:   Due to round-off errors, some or all basis vectors in the kernel of `A` may be missed in the numerical modes.

The special solution `X` in conjunction with `KernelBasis` provides the general solution of AX = B. Solutions generated without the option `ShowAssumptions` are valid for arbitrary complex values of the symbolic parameters which may be present in `A` and `B`. If no such solution exists, then `[FAIL,NIL]` is returned. Solutions that are valid only for special values of the symbolic parameters may be obtained with `ShowAssumptions`. See Example 3, Example 4, and Example 7.

`numeric::matlinsolve` internally uses a sparse representation of the matrices. It is suitable for solving large sparse systems. See Example 5.

 Note:   `numeric::matlinsolve` does not react to any assumptions on symbolic parameters in `A,B` that are set via `assume`.
 Note:   Gaussian elimination with partial pivoting is used. Without the option `Symbolic`, the pivoting strategy takes care of numerical stabilization. With `Symbolic`, exact data are assumed. The symbolic pivoting strategy tries do maximize speed and does not take care of numerical stabilization! Do not use `Symbolic` for linear systems involving floating-point entries! See Example 6.
 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.Note that `linalg::matlinsolve` must be used, when the solution is to be computed over the component domain. See . Example 8.

## Environment Interactions

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

## Examples

### Example 1

We use equivalent input formats ```(B1, B2)``` to represent a vector with components [a, π]. First, this vector is defined as a 2-dimensional array:

```A := array(1..2, 1..3, [[1, 2, 3],[1, 1, 2]]): B1 := array(1..2, 1..1, [[a], [PI]]): numeric::matlinsolve(A, B1)```

Next, we use a 1-dimensional array and compute an exact solution:

```B2 := array(1..2, [a, PI]): numeric::matlinsolve(A, B2, Symbolic)```

Now, a list is used to specify the vector. No internal assumptions were used by `numeric::matlinsolve` to obtain the solution:

```B3 := [a, PI]: numeric::matlinsolve(A, B3, ShowAssumptions)```

Finally, we use `Dom::Matrix` objects to specify the system. Note that the results are returned as corresponding matrix objects:

```A := matrix([[1, 2, 3],[1, 1, 2]]): B4 := matrix([a, PI]): numeric::matlinsolve(A, B4)```

`delete A, B1, B2, B3, B4:`

### Example 2

We invert a matrix by solving AX = 1:

```A := hfarray(1..3, 1..3, [[1, 1, 0], [0, 1, 1], [0, 0, 1]]): B := matrix::identity(3, 3): InverseOfA := numeric::matlinsolve(A, B, Symbolic)[1]```

`delete A, B, InverseOfA:`

### Example 3

We solve an equation with a symbolic parameter x:

```A := matrix([[2, 2, 3], [1, 1, 2], [3, 3, 5]]): B := matrix([sin(x)^2, cos(x)^2, 0]): [X, Kernel, Constraints, Pivots] := numeric::matlinsolve(A, B, Symbolic, ShowAssumptions)```

This solution holds subject to the constraint sin(x)2 + cos(x)2 = 0 on the parameter x. `numeric::matlinsolve` does not investigate the `Constraints` and does not realize that they cannot be satisfied. We check the consistency of the "result" by inserting the solution into the original system. Since the input matrix A was of type `Dom::Matrix()`, the results `X` and `Kernel` were returned as corresponding matrices. The overloaded operators `*` and `-` for matrix multiplication and subtraction can be used:

`A*X - B, A*Kernel`

`delete A, B, X, Kernel, Constraints, Pivots:`

### Example 4

We give a further demonstration of the option `ShowAssumptions`. The following system does not have a solution for all values of the parameter `a`:

```A := array(1..2, 1..2, [[1, 1], [1, 1]]): B := array(1..2, 1..1, [[1], [a]]): numeric::matlinsolve(A, B, Symbolic)```

With `ShowAssumptions`, `numeric::matlinsolve` investigates under which conditions (on the parameter a `a`) there is a solution:

`numeric::matlinsolve(A, B, Symbolic, ShowAssumptions)`

We conclude that there is a 1-dimensional solution space for a = 1.

`delete A, B:`

### Example 5

We solve a sparse system with 3 diagonal bands:

```n := 100: A := matrix(n, n, [1, -2, 1], Banded): B := array(1..n, [1 \$ n]): numeric::matlinsolve(A, B)```

[matrix([[-50.0], [-99.0], [-147.0], [-194.0], [-240.0], [-285.0], [-329.0], [-372.0], [-414.0], [-455.0], [-495.0], [-534.0], [-572.0], [-609.0], [-645.0], [-680.0], [-714.0], [-747.0], [-779.0], [-810.0], [-840.0], [-869.0], [-897.0], [-924.0], [-950.0], [-975.0], [-999.0], [-1022.0], [-1044.0], [-1065.0], [-1085.0], [-1104.0], [-1122.0], [-1139.0], [-1155.0], [-1170.0], [-1184.0], [-1197.0], [-1209.0], [-1220.0], [-1230.0], [-1239.0], [-1247.0], [-1254.0], [-1260.0], [-1265.0], [-1269.0], [-1272.0], [-1274.0], [-1275.0], [-1275.0], [-1274.0], [-1272.0], [-1269.0], [-1265.0], [-1260.0], [-1254.0], [-1247.0], [-1239.0], [-1230.0], [-1220.0], [-1209.0], [-1197.0], [-1184.0], [-1170.0], [-1155.0], [-1139.0], [-1122.0], [-1104.0], [-1085.0], [-1065.0], [-1044.0], [-1022.0], [-999.0], [-975.0], [-950.0], [-924.0], [-897.0], [-869.0], [-840.0], [-810.0], [-779.0], [-747.0], [-714.0], [-680.0], [-645.0], [-609.0], [-572.0], [-534.0], [-495.0], [-455.0], [-414.0], [-372.0], [-329.0], [-285.0], [-240.0], [-194.0], [-147.0], [-99.0], [-50.0]]), 0]
`delete n, A, B:`

### Example 6

The option `Symbolic` should not be used for equations with floating-point coefficients, because the symbolic pivoting strategy favors efficiency instead of numerical stability.

```A := array(1..2, 1..2, [[1, 10^20], [1, 1]]): B := array(1..2, 1..1, [[10^20], [0]]):```

The float approximation of the exact solution is:

`map(numeric::matlinsolve(A, B, Symbolic)[1], float)`

We now convert the exact input data to floating-point approximations:

`A := map(A, float): B := map(B, float):`

The default pivoting strategy of the floating-point algorithm stabilizes floating-point operations. Consequently, one gets a correct result:

`numeric::matlinsolve(A, B)[1]`

With the option `Symbolic`, however, the pivoting strategy optimizes speed, assuming exact arithmetic. Numerical instabilities may occur if floating-point coefficients are involved. The following result is caused by internal round-off effects ("cancellation"):

`numeric::matlinsolve(A, B, Symbolic)[1]`

We need to increase `DIGITS` to obtain a better result:

```DIGITS := 20: numeric::matlinsolve(A, B, Symbolic)[1]```

`delete A, B, DIGITS:`

### Example 7

We demonstrate how a complete solution of the following linear system with symbolic parameters may be found:

```A := array(1..3, 1..2, [[1, 1], [a, b], [1, c]]): B := array(1..3, 1..1, [[1], [1], [1]]): numeric::matlinsolve(A, B, Symbolic, ShowAssumptions)```

This is the general solution assuming ab. We now set b = a to investigate further solution branches:

```A := subs(A, b = a): numeric::matlinsolve(A, B, Symbolic, ShowAssumptions)```

This is the general solution for a = b, assuming c ≠ 1. We finally set c = 1 to obtain the last solution branch:

```A := subs(A, c = 1): numeric::matlinsolve(A, B, Symbolic, ShowAssumptions)```

From the constraint on `a`, we conclude that there is a 1-dimensional solution space for the special values a = b = c = 1 of the symbolic parameters.

`delete A, B:`

### Example 8

Matrices from a domain such as `Dom::Matrix``(...)` are converted to arrays with numbers or expressions. Hence, `numeric::matlinsolve` finds no solution for the following system:

```M := Dom::Matrix(Dom::IntegerMod(7)): A := M([[1, 4], [6, 3], [3, 2]]): B := M([[9], [5], [0]]): numeric::matlinsolve(A, B)```

Use `linalg::matlinsolve` to solve the system over the coefficient field of the matrices. A solution does exist over the field `Dom::IntegerMod``(7)`:

`linalg::matlinsolve(A, B)`

`delete M, A, B:`

### Example 9

We demonstrate the difference between `Symbolic`, `HardwareFloats`, and `SoftwareFloats`. The following matrix A has a 1-dimensional kernel. Due to round-off, a further spurious kernel vector appears with `SoftwareFloats`. No kernel vector is detected with `HardwareFloats`:

```A := matrix([[2*10^14 + 2, 2*10^(-9), 2*10^(-4)], [3*10^15 + 3, 3*10^(-8), 3*10^(-3)], [4*10^16 + 4, 4*10^(-7), 4*10^(-2)] ]): b := matrix([2*10^(-9), 3*10^(-8), 4*10^(-7)]): float(numeric::matlinsolve(A, b, Symbolic)); numeric::matlinsolve(A, b, SoftwareFloats); numeric::matlinsolve(A, b, HardwareFloats)```

`delete A, b:`

## Parameters

 `A` An m×n matrix of domain type `DOM_ARRAY`, `DOM_HFARRAY`, or of category `Cat::Matrix` `B` An m×p matrix of domain type `DOM_ARRAY`, `DOM_HFARRAY`, or of category `Cat::Matrix`. Column vectors `B` may also be represented by a 1-dimensional ```array(1..m, [B1, B2, …] )```, a 1-dimensional ```hfarray(1..m, [B1, B2, …] )```, or by a list ```[B1, B2, …]```.

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

 Note:   For ill-conditioned matrices, the results returned with `HardwareFloats` and `SoftwareFloats` may differ significantly! See Example 9.

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

 Note:   This option should not be used for matrices with floating-point entries! Numerical instabilities may occur in floating-point operations. Cf. Example 6.

`ShowAssumptions`

Returns information on internal assumptions on symbolic parameters in `A` and `B`. With this option, either exact arithmetic or `SoftwareFloats` are used.

This option is only useful if the matrices contain symbolic parameters. Consequently, it should only be used in conjunction with the option `Symbolic`.

 Note:   This option changes the format of the return value to ```[X, KernelBasis, Constraints, Pivots]```.

`X` and `KernelBasis` represent the general solution subject to `Constraints` and `Pivots`.

`Constraints` is a list of equations for symbolic parameters in B which are necessary and sufficient for AX = B to be solvable.

Such constraints arise if Gaussian elimination leads to equations of the form 0 = c, where c is some expression involving symbolic parameters contained in B. All such equations are collected in `Constraints`; `numeric::matlinsolve` assumes that these equations are satisfied and returns a special solution `X`.

If no such constraints arise, the return value of `Constraints` is the empty list.

`Pivots` is a list of inequalities involving symbolic parameters in A. Internally, division by pivot elements occurs in the Gaussian elimination. The expressions collected in `Pivots` are the numerators of those pivot elements that involve symbolic parameters contained in A. If only numerical pivot elements are used, then the return value of `Pivots` is the empty list.

 Note:   `Constraints` usually is a list of non-linear equations for the symbolic parameters. It is not investigated by `numeric::matlinsolve`, i.e., solutions may be returned, even if the `Constraints` cannot be satisfied. See Example 3.
 Note:   This option changes the return strategy for "unsolvable" systems. Without the option `ShowAssumptions`, the result `[FAIL,NIL]` is returned, whenever Gaussian elimination produces an equation 0 = c with non-zero c. With `ShowAssumptions`, such equations are returned via `Constraints`, provided c involves symbolic parameters. If c is a purely numerical value, then `[FAIL, NIL, [], []]` is returned.

See Example 3, Example 4, and Example 7.

`NoWarning`

Suppresses warnings

If symbolic coefficients are found, `numeric::matlinsolve` automatically switches to the `Symbolic` mode with a warning. With this option, this warning is suppressed; `numeric::matlinsolve` still uses the symbolic mode for symbolic coefficients, i.e., exact arithmetic without floating-point conversions is used.

`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`, `DOM_HFARRAY`, `Dom::Matrix()`, or `Dom::DenseMatrix()`.

`Sparse`

Use a sparse internal representation for matrices.

This option only has an effect when used in conjunction with `HardwareFloats`. With the `Sparse` option, the linear solver uses a sparse representation of the matrices to save memory and increase efficiency. However, if the coefficient matrix is not sparse, this option will cost some additional memory and runtime.

## Return Values

Without the option `ShowAssumptions`, a list ```[X, KernelBasis]``` is returned. The (special) solution `X` is an n×p matrix. `KernelBasis` is an n×d matrix (d is the dimension of the kernel of A). Its columns span the kernel of A. If the kernel is trivial, `KernelBasis` is the integer 0.

`[FAIL, NIL]` is returned if the system is not solvable.

With `ShowAssumptions`, a list ```[X, KernelBasis, Constraints, Pivots]``` is returned. The lists `Constraints` and `Pivots` contain equations and inequalities involving symbolic parameters in `A` and `B`. Internally these were assumed to hold true when solving the system. ```[FAIL, NIL, [], []]``` is returned if the system is not solvable.