## `numeric::matlinsolve` -- solve a linear matrix equation

#### Introduction

`numeric::matlinsolve``(A, B, ..)` returns the matrix solution X of the matrix equation A*X=B.

#### Call(s)

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

#### Parameters

`A` | - | an m x n matrix of domain type `DOM_ARRAY` or of category `Cat::Matrix`.
`B` | - | an m x 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,..]`.

#### Options

`Symbolic` | - | prevents conversion of input data to floating point numbers
`ShowAssumptions` | - | returns information on internal assumptions on symbolic parameters in `A` and `B`

#### Returns

Without the option ShowAssumptions a list `[X,KernelBasis]` is returned. The solution `X` is an n x p array. `KernelBasis` is an n x d array. Its columns span the kernel of A. `[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.

#### Side Effects

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

#### Related Functions

`linalg::matlinsolve`, `linsolve`, `numeric::inverse`, `numeric::linsolve`, `solve`

#### Details

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

- 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, then `numeric::matlinsolve` automatically switches to its symbolic mode, issuing a warning. Symbolic parameters in B are accepted without warning.

- The numerical working precision is set by the environment variable `DIGITS`.
- `X` is a special solution of A*X=B.
- The n x d array `KernelBasis` is the most general solution of A*X=0. Its columns span the d-dimensional kernel of A.
- If the kernel is empty, then the return value of `KernelBasis` is the null vector represented by `array(1..n,1..1,[0,..,0])`.
- If `KernelBasis=array(1..n,1..1,[0,..,0])`, then the dimension d of the kernel of A is `d:=0`. Otherwise it is `d:=op(KernelBasis,[0,3,2]])`.

- Due to roundoff errors some or all basis vectors in the kernel of `A` may be missed in the numeric mode.

- The special solution `X` in conjunction with `KernelBasis` provides the general solution of A*X=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. Cf. examples 3, 4, 7.
- `numeric::matlinsolve` internally uses a sparse representation of the matrices. It is suitable for solving large sparse systems. Cf. example 5. Note, however, that the input requires specification of entire matrices. For this reason you should consider using `numeric::linsolve` which allows to input the linear system as a sparse list of symbolic equations.

- `numeric::matlinsolve` does not react to any assumptions on symbolic parameters in `A,B` that are set via `assume`.

- Gaussian elimination with partial pivoting is used. Without 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! Cf. example 6.

- `Cat::Matrix` objects `A` from matrix domains such as `Dom::Matrix``(..)` or `Dom::SquareMatrix``(..)` are internally converted to arrays over expressions via `A::dom::expr(A)`. Note that `linalg::matlinsolve` must be used, when the solution is to be computed over the component domain. Cf. example 8. Note that the option Symbolic should be used, if the entries cannot be converted to numerical expressions.

#### Option: Symbolic

- This option prevents conversion of the input data to floats.
- This option must be used, if the matrix A contains symbolic parameters that cannot be converted to floating point numbers.

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

#### Option: ShowAssumptions

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

- 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, then the return value of `Constraints` is the empty list.

- `Constraints` usually is a non-linear list of equation for the symbolic parameters. It is not investigated by `numeric::matlinsolve`, i.e., solutions may be returned, even if the `Constraints` cannot be satisfied. Cf. example 3.

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

- `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.
- Cf. examples 3, 4, 7.

#### Example 1

We use equivalent input formats `B1`,...,`B4` to represent a vector with components [a,PI]. 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)

                    -- +-                     -+  +-      -+ --
                    |  |  6.283185307 - 1.0 a  |  |  -1.0  |  |
                    |  |                       |  |        |  |
                    |  |  1.0 a - 3.141592654  |, |  -1.0  |  |
                    |  |                       |  |        |  |
                    |  |           0           |  |    1   |  |
                    -- +-                     -+  +-      -+ --

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

    >> B2 := array(1..2, [a, PI]):

    >> numeric::matlinsolve(A, B2, Symbolic)

                          -- +-          -+  +-    -+ --
                          |  |  2 PI - a  |  |  -1  |  |
                          |  |            |  |      |  |
                          |  |   a - PI   |, |  -1  |  |
                          |  |            |  |      |  |
                          |  |      0     |  |   1  |  |
                          -- +-          -+  +-    -+ --

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)

                -- +-                     -+  +-      -+         --
                |  |  6.283185307 - 1.0 a  |  |  -1.0  |          |
                |  |                       |  |        |          |
                |  |  1.0 a - 3.141592654  |, |  -1.0  |, [], []  |
                |  |                       |  |        |          |
                |  |           0           |  |    1   |          |
                -- +-                     -+  +-      -+         --

Finally, we use `Dom::Matrix` objects to specify the system. Note that the result are still arrays:

    >> A := Dom::Matrix()([[1, 2, 3],[1, 1, 2]]):

    >> B4 := Dom::Matrix()([a, PI]):

    >> numeric::matlinsolve(A, B4)

                    -- +-                     -+  +-      -+ --
                    |  |  6.283185307 - 1.0 a  |  |  -1.0  |  |
                    |  |                       |  |        |  |
                    |  |  1.0 a - 3.141592654  |, |  -1.0  |  |
                    |  |                       |  |        |  |
                    |  |           0           |  |    1   |  |
                    -- +-                     -+  +-      -+ --

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

#### Example 2

We invert a matrix by solving A*X=1:

    >> A := array(1..3, 1..3, [[1, 1, 0], [0, 1, 1], [0, 0, 1]]):

    >> B := array(1..3, 1..3, [[1, 0, 0], [0, 1, 0], [0, 0, 1]]):

    >> InverseOfA := numeric::matlinsolve(A, B, Symbolic)[1]

                                  +-           -+
                                  |  1, -1,  1  |
                                  |             |
                                  |  0,  1, -1  |
                                  |             |
                                  |  0,  0,  1  |
                                  +-           -+

    >> delete A, B, InverseOfA:

#### Example 3

We solve an equation with a symbolic parameter x:

    >> M := Dom::Matrix():

    >> A := M(3, 3, [[2, 2, 3], [1, 1, 2], [3, 3, 5]]):

    >> B := M(3, 1, [sin(x)^2, cos(x)^2, 0]):

    >> [X, Kernel, Constraints, Pivots] :=
       numeric::matlinsolve(A, B, Symbolic, ShowAssumptions)

          -- +-             -+                                        --
          |  |           2   |  +-    -+                               |
          |  |   5 sin(x)    |  |  -1  |                               |
          |  |               |  |      |         2         2           |
          |  |       0       |, |   1  |, [cos(x)  + sin(x)  = 0], []  |
          |  |               |  |      |                               |
          |  |            2  |  |   0  |                               |
          |  |  - 3 sin(x)   |  +-    -+                               |
          -- +-             -+                                        --

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. For convenience we convert the arrays `X` and `Kernel` to `Dom::Matrix` objects and use the overloaded operators `*` and `-` for matrix multiplication and subtraction:

    >> A*M(X) - B, A*M(Kernel)

                        +-                     -+
                        |           0           |  +-   -+
                        |                       |  |  0  |
                        |          2         2  |  |     |
                        |  - cos(x)  - sin(x)   |, |  0  |
                        |                       |  |     |
                        |           0           |  |  0  |
                        +-                     -+  +-   -+

    >> delete M, 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)

                                    [FAIL, NIL]

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

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

                     -- +-   -+  +-    -+                  --
                     |  |  1  |  |  -1  |                   |
                     |  |     |, |      |, [a - 1 = 0], []  |
                     |  |  0  |  |   1  |                   |
                     -- +-   -+  +-    -+                  --

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

    >> delete A, B:

#### Example 5

We solve a large sparse system with 3 diagonal bands:

    >> n := 100: A := Dom::Matrix()(n, n, [1, -2, 1], Banded):

    >> B := array(1..n, [1 $ n]): numeric::matlinsolve(A, B)

                           -- +-         -+  +-   -+ --
                           |  |   -50.0   |  |  0  |  |
                           |  |           |  |     |  |
                           |  |   -99.0   |  |  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)

                                    +-      -+
                                    |  -1.0  |
                                    |        |
                                    |   1.0  |
                                    +-      -+

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

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

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

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

                                    +-      -+
                                    |  -1.0  |
                                    |        |
                                    |   1.0  |
                                    +-      -+

With the option Symbolic 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]

                                     +-     -+
                                     |   0   |
                                     |       |
                                     |  1.0  |
                                     +-     -+

    >> delete A, B:

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

          -- +-       -+                                               --
          |  |  b - 1  |                                                |
          |  |  -----  |  +-   -+                                       |
          |  |  b - a  |  |  0  |                                       |
          |  |         |, |     |, [a c - c - a + 1 = 0], [b - a <> 0]  |
          |  |  1 - a  |  |  0  |                                       |
          |  |  -----  |  +-   -+                                       |
          |  |  b - a  |                                                |
          -- +-       -+                                               --

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

    >> A := subs(A, b = a):

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

                 -- +-   -+  +-   -+                            --
                 |  |  1  |  |  0  |                             |
                 |  |     |, |     |, [1 - a = 0], [c - 1 <> 0]  |
                 |  |  0  |  |  0  |                             |
                 -- +-   -+  +-   -+                            --

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)

                     -- +-   -+  +-    -+                  --
                     |  |  1  |  |  -1  |                   |
                     |  |     |, |      |, [1 - a = 0], []  |
                     |  |  0  |  |   1  |                   |
                     -- +-   -+  +-    -+                  --

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)

                                    [FAIL, NIL]

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)

                                   +-         -+
                                   |  1 mod 7  |
                                   |           |
                                   |  2 mod 7  |
                                   +-         -+

    >> delete M, A, B:

#### Changes

- new function

Do you have questions or comments?
Copyright © SciFace Software GmbH & Co. KG 2000