## `linalg::matlinsolve` -- solving systems of linear equations

#### Introduction

`linalg::matlinsolve``(A, b)` computes the general solution of the equation A*x = b.

#### Call(s)

`linalg::matlinsolve(A, b)`
`linalg::matlinsolve(A, b, list)`
`linalg::matlinsolve(A, b <, Special> <, Unique>)`
`linalg::matlinsolve(A, B)`
`linalg::matlinsolve(A)`

#### Parameters

`A` | - | an m x n matrix of a domain of category `Cat::Matrix`
`B` | - | an m x k matrix of a domain of category `Cat::Matrix`
`b` | - | an m-dimensional column vector, i.e., a m x 1 matrix of a domain of category `Cat::Matrix`
`list` | - | a list of n elements of the component ring of `A`

#### Options

`Special` | - | Computes one particular solution of the system A*x = 0.
`Unique` | - | Checks whether the system has a unique solution and returns the solution, or the value `NIL` otherwise.

#### Returns

a vector, a list `[s, kern]` (possibly empty), where `s` is a solution vector and `kern` is a list of basis vectors for the kernel of `A`, a matrix, or the value `NIL`.

The matrix and the vectors, respectively, are of the domain type `Dom::Matrix(R)`, where `R` is the component ring of `A`.

#### Related Functions

`linsolve`, `linalg::expr2Matrix`, `linalg::nullspace`, `linalg::matlinsolveLU`, `linalg::wiedemann`, `numeric::matlinsolve`

#### Details

- `linalg::matlinsolve``(A, b)` returns the solution vector x of the system A*x = b if it is a unique solution.
- `linalg::matlinsolve``(A, b)` returns a list [w, [v[1],...,v[r]]] if the system A*x = b has more than one solution, where w is one particular solution, i.e., A*w = b and v[1],...,v[r] form a basis of the kernel of `A`, i.e., the solution space of the homogenous system A*x = 0.

  Each solution x has the form w + s[1]*v[1] + ... + s[r]*v[r] (r <= n) with certain scalars s[1],...,s[r].

- A list of n scalars [s[1],...,s[n]] may be passed as the additional parameter `list`. This extracts the solution w + s[i[1]]*v[1] + ... + s[i[r]]*v[r] with {i[1],...,i[r]}= {1,...,n}minus {j[1],...,j[l]} from the solution space of the system A*x = b, where j[1],...,j[l] are the characteristic column indices of `A` (see `linalg::gaussJordan`).

  The entries of `list` are converted to elements of the component ring of `A` (an error message is returned if this is not possible).

- If the system A*x = b has no solution, then the empty list `[]` is returned.
- `linalg::matlinsolve``(A)` solves the matrix equation C*x = b, where b is the last column of `A` and C is `A` with the last column deleted.
- `linalg::matlinsolve``(A, B)` returns the solution X of the matrix equation A*X = B, if it has exactly one solution. Otherwise the empty list `[]` is returned.
- The vector `b` and the matrix `B` respectively, are converted into the domain `Dom::Matrix(R)`, where `R` is the component ring of `A`. Solution vectors also belong to this domain.
- The component ring of `A` must be an integral domain, i.e., a domain of category `Cat::IntegralDomain`.
- `linalg::matlinsolve` can compute the general solution for systems with more than one solution only over fields, i.e., component rings of category `Cat::Field`. If in this case the component ring of `A` does not have a canonical representation of the zero element, then it may happen that `linalg::matlinsolve` does not find a basis for the null space. In such a case, a wrong result is returned.
- `linalg::matlinsolve` does not exploit the structure of `A`, e.g., sparsity. A matrix is sparse if it has many zero components (see example 5).
- To get a floating-point approximation use the function `numeric::matlinsolve`. Also if the input matrices are defined over the component ring `Dom::Float`, we recommend for efficieny reasons to use `numeric::matlinsolve` instead of `linalg::matlinsolve`!

#### Option: Special

- Only one particular solution `w` of the system A*x = b is returned. This supresses the computation of a basis for the kernel of `A`.

#### Option: Unique

- Checks whether the system has a unique solution and returns it. The return value `NIL` means that the system has more than one solution.

#### Example 1

We solve the linear system:

    +-       -+       +-    -+
    |   1, 2  |       |   1  |
    |         | * x = |      |
    |  -1, 2  |       |  -1  |
    +-       -+       +-    -+

over the reals. First we enter the coefficient matrix and the right-hand side:

    >> MatR := Dom::Matrix(Dom::Real):
       A := MatR([[1, 2], [-1, 2]]); b := MatR([1, -1])

                                    +-       -+
                                    |   1, 2  |
                                    |         |
                                    |  -1, 2  |
                                    +-       -+

                                     +-    -+
                                     |   1  |
                                     |      |
                                     |  -1  |
                                     +-    -+

Next we call `linalg::matlinsolve` to solve the system:

    >> x := linalg::matlinsolve(A, b)

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

We see that the system has exactly one solution. The vector `x` satisfies the matrix equation given above:

    >> A * x

                                     +-    -+
                                     |   1  |
                                     |      |
                                     |  -1  |
                                     +-    -+

#### Example 2

The system:

    +-        -+       +-   -+
    |   1,  2  |       |  1  |
    |          | * x = |     |
    |  -1, -2  |       |  0  |
    +-        -+       +-   -+

does not have a solution over R (in fact, over no component domain):

    >> A := MatR([[1, 2], [-1, -2]]): b := MatR([1, 0]):
       linalg::matlinsolve(A, b)

                                        []

#### Example 3

We solve the linear system:

    +-                  -+       +-    -+
    |  1, 1, -4, -7, -6  |       |  30  |
    |                    | * x = |      |
    |  0, 1, -3, -5, -7  |       |  70  |
    +-                  -+       +-    -+

over the rational numbers. First we enter the coefficient matrix and the right-hand side:

    >> MatQ := Dom::Matrix(Dom::Rational):
       A := MatQ([[1, 1, -4, -7, -6], [0, 1, -3, -5, -7]]);
       b := MatQ([30, 17])

                              +-                  -+
                              |  1, 1, -4, -7, -6  |
                              |                    |
                              |  0, 1, -3, -5, -7  |
                              +-                  -+

                                     +-    -+
                                     |  30  |
                                     |      |
                                     |  17  |
                                     +-    -+

Next we call `linalg::matlinsolve` to solve the system:

    >> sol:= linalg::matlinsolve(A, b)

                 -- +-    -+  -- +-   -+  +-   -+  +-    -+ -- --
                 |  |  13  |  |  |  1  |  |  2  |  |  -1  |  |  |
                 |  |      |  |  |     |  |     |  |      |  |  |
                 |  |  17  |  |  |  3  |  |  5  |  |   7  |  |  |
                 |  |      |  |  |     |  |     |  |      |  |  |
                 |  |   0  |, |  |  1  |, |  0  |, |   0  |  |  |
                 |  |      |  |  |     |  |     |  |      |  |  |
                 |  |   0  |  |  |  0  |  |  1  |  |   0  |  |  |
                 |  |      |  |  |     |  |     |  |      |  |  |
                 |  |   0  |  |  |  0  |  |  0  |  |   1  |  |  |
                 -- +-    -+  -- +-   -+  +-   -+  +-    -+ -- --

The result is to be interpreted as follows: The first vector of the list `sol` is a particular solution of the linear system:

    >> A * sol[1]

                                     +-    -+
                                     |  30  |
                                     |      |
                                     |  17  |
                                     +-    -+

The second entry of the list contains a basis for the null space of A, i.e., the solution space of the corresponding homogenous system A * x = 0 (the kernel of A). The basis returned is given as a list of vectors.

The following input checks this fact by computing the product A * x for each vector x of the list `sol[2]`:

    >> map(sol[2], x -> A * x)

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

Any solution of the linear system can be represented as a sum of a particular solution (here: `sol[1]`) and a linear combination of the basis vectors of the kernel of A. Hence our input system has an infinite number of solutions.

For example, another solution of the system is given by:

    >> x := sol[1] + 1*sol[2][1] + 1/2*sol[2][2] - 2*sol[2][3]

                                    +-      -+
                                    |   17   |
                                    |        |
                                    |  17/2  |
                                    |        |
                                    |    1   |
                                    |        |
                                    |   1/2  |
                                    |        |
                                    |   -2   |
                                    +-      -+

    >> A * x

                                     +-    -+
                                     |  30  |
                                     |      |
                                     |  17  |
                                     +-    -+

If we identify the columns of the coefficient matrix A of our linear system with the variables x1,x2,x3,x4,x5, then we see from the general solution that the variables x3,x4,x5 act as free parameters. They can be assigned arbitrary rational values to obtain a unique solution.

By giving a list of values for these variables as a third parameter to `linalg::matlinsolve`, we can select a certain vector from the set of all solutions of the linear system. For example, to select the same vector `x` as chosen in the previous input, we enter:

    >> linalg::matlinsolve(A, b, [0, 0, 1, 1/2, -2])

                                    +-      -+
                                    |   17   |
                                    |        |
                                    |  17/2  |
                                    |        |
                                    |    1   |
                                    |        |
                                    |   1/2  |
                                    |        |
                                    |   -2   |
                                    +-      -+

If one is only interested in a particular solution and does not need the general solution of the linear system, one may enter:

    >> linalg::matlinsolve(A, b, Special)

                                     +-    -+
                                     |  13  |
                                     |      |
                                     |  17  |
                                     |      |
                                     |   0  |
                                     |      |
                                     |   0  |
                                     |      |
                                     |   0  |
                                     +-    -+

This call suppresses the computation of the kernel of A.

#### Example 4

If the linear system is given in form of equations the function `linalg::expr2Matrix` can be used to form the corresponding matrix equation:

    >> delete x, y, z:
       Ab := linalg::expr2Matrix(
         [x + y + z = 6, 2*x + y + 2*z = 10, x + 3*y + z = 10]
       )

                                 +-             -+
                                 |  1, 1, 1,  6  |
                                 |               |
                                 |  2, 1, 2, 10  |
                                 |               |
                                 |  1, 3, 1, 10  |
                                 +-             -+

The result here is the extended coefficient matrix of the input system, i.e., the right-hand side vector b is the 4th column vector of the matrix `Ab`. Since we did not specify a component ring for this matrix, the standard component ring for matrices, the domain `Dom::ExpressionField()`, was chosen.

To solve the linear system, we call:

    >> linalg::matlinsolve(Ab)

                           -- +-   -+  -- +-    -+ -- --
                           |  |  4  |  |  |  -1  |  |  |
                           |  |     |  |  |      |  |  |
                           |  |  2  |, |  |   0  |  |  |
                           |  |     |  |  |      |  |  |
                           |  |  0  |  |  |   1  |  |  |
                           -- +-   -+  -- +-    -+ -- --

We see that the system has an infinite number of solutions. The third variable z acts as a free parameter and therefore can have any (complex) value.

To get the general solution in parameter form, one may use parameters for the variables x, y, z of the input system:

    >> delete u, v, w:
       sol := linalg::matlinsolve(Ab, [u, v, w])

                                   +-         -+
                                   |  - w + 4  |
                                   |           |
                                   |     2     |
                                   |           |
                                   |     w     |
                                   +-         -+

This is possible here because we perform the matrix computations over `Dom::ExpressionField()` which allows to compute with symbolical (arithmetical) expressions.

To select a certain vector from the set of solutions, for example, the solution for w=1, we enter:

    >> x := subs(sol, w = 1)

                                      +-   -+
                                      |  3  |
                                      |     |
                                      |  2  |
                                      |     |
                                      |  1  |
                                      +-   -+

#### Example 5

Suppose that we have a system of linear equations with a sparse structure, i.e., with a coefficient matrix having many zero components. For example:

    >> eqs := {x1 + x5 = 0, x2 - x4 = 1, x3 + 2*x5 = 2, x4 - x5 = -1}:
       Ab := linalg::expr2Matrix(eqs)

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

As `linalg::matlinsolve` does not exploit the sparsity of the given coefficient matrix, we recommend not to solve the system with `linalg::matlinsolve`, but to use the function `linsolve` which directly works on the equations and therefore preserves the sparse structure of the input system:

    >> linsolve(eqs)

                [x3 = 2 x1 + 2, x4 = - x1 - 1, x2 = -x1, x5 = -x1]

If a system is given in matrix form, we recommend to use the function `numeric::matlinsolve` with option `Symbolic` instead of `linalg::matlinsolve` even for exact computations. This function works more efficiently than `linalg::matlinsolve` and is able to exploit sparsity:

    >> A := linalg::delCol(Ab, 5): b := linalg::col(Ab, 6):
       numeric::matlinsolve(A, b)

                           -- +-      -+  +-      -+ --
                           |  |   2.0  |  |  -2.0  |  |
                           |  |        |  |        |  |
                           |  |  -1.0  |  |   1.0  |  |
                           |  |        |  |        |  |
                           |  |    0   |, |    0   |  |
                           |  |        |  |        |  |
                           |  |    0   |  |    0   |  |
                           |  |        |  |        |  |
                           |  |    0   |  |    1   |  |
                           -- +-      -+  +-      -+ --

Note that the function `numeric::matlinsolve` always works over a subfield of the complex numbers and does not allow to specify the domain of computation. Without the option `Symbolic` `numeric::matlinsolve` converts input data to floating point numbers.

#### Example 6

Let us check whether the matrix equation

    +-       -+       +-      -+
    |   1, 2  |       |  4, 2  |
    |         | * X = |        |
    |  -2, 3  |       |  6, 3  |
    +-       -+       +-      -+

has a unique solution over the integers.

We start by entering the coefficient matrix and the right-hand side matrix:

    >> MatZ := Dom::Matrix(Dom::Integer):
       A := MatZ([[1, 2], [-2, 3]]); B := MatZ([[4, 2], [6, 3]])

                                    +-       -+
                                    |   1, 2  |
                                    |         |
                                    |  -2, 3  |
                                    +-       -+

                                    +-      -+
                                    |  4, 2  |
                                    |        |
                                    |  6, 3  |
                                    +-      -+

Next we solve the matrix equation:

    >> X := linalg::matlinsolve(A, B)

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

The equation indeed has a unique solution (otherwise the answer of `linalg::matlinsolve` would be the empty list `[]`). Let us check the result:

    >> A * X

                                    +-      -+
                                    |  4, 2  |
                                    |        |
                                    |  6, 3  |
                                    +-      -+

#### Background

- Let A be an m x n matrix with components from a field F and b an m-dimensional vector over F. Let (A,b) be the extended coefficient matrix of the linear system A x = b.

  Then the following holds:

  - The linear system A*x = b has a solution, if and only if rank(A,b) = rank(A).
  - It has exactly one solution, if and only if rank(A,b) = rank(A) = n.
  - If w is a solution of the system A x = b and {v[1],...,v[r]} a basis of the kernel of A, then

        L(A,b) = {w + s[1]*v[1] + ... + s[r]*v[r] | s[1],...,s[r] from F }

    is the set of all solutions of the linear system A x = b, the general solution of the (inhomogeneus) linear system.
- The kernel of the matrix A is defined as:

      ker(A):= {w | A*w = 0 }.

  The kernel of A is a vector space over F of dimension n - rang(A).

#### Changes

- used to be called `linalg::linearSolve`

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