## `linalg::factorLU` -- LU-decomposition of a matrix

#### Introduction

`linalg::factorLU``(A)` computes an LU-decomposition of an m x n matrix A, i.e., a decomposition of the A into an m x m lower triangular matrix L and an n x m upper triangular matrix U such that P*A = L*U, where P is a permutation matrix.

#### Call(s)

`linalg::factorLU(A)`

#### Parameters

`A` | - | a matrix of a domain of category `Cat::Matrix`

#### Returns

a list `[L, U, pivindex]` with the two matrices L and U of the domain `Dom::Matrix(R)` and a list `pivindex` of positive integers. `R` is the component ring of `A`.

#### Related Functions

`linalg::factorQR`, `linalg::factorCholesky`, `linalg::inverseLU`, `linalg::matlinsolveLU`, `lllint`, `numeric::factorLU`

#### Details

- The diagonal entries of the lower triangular matrix L are equal to one (Doolittle-decomposition). The diagonal entries of U are the pivot elements used during the computation.

  The matrices L and U are unique.

- `pivindex` is a list `[r[1], r[2], ...]` representing the row exchanges of A in the pivoting steps, i.e., B = P*A = L*U, where B[i,j] = A[r[i],j].
- A floating-point approximation of the decomposition is computed using `numeric::factorLU`, if the matrix `A` is defined over the component ring `Dom::Float`. In this case it is recommended to call `numeric::factorLU` directly for a better efficiency.
- The algorithm also works for singular A. In this case either L or U is singular.
- L and U are nonsingular if and only if A is nonsingular.
- The component ring of the matrix `A` must be a field, i.e., a domain of category `Cat::Field`.

#### Example 1

We compute an LU-decomposition of the real matrix:

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

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

    >> [L, U, pivlist] := linalg::factorLU(A)

              -- +-             -+  +-              -+            --
              |  |   1,   0,  0  |  |  2,  -3,  -1   |             |
              |  |               |  |                |             |
              |  |  1/2,  1,  0  |, |  0, 5/2, -1/2  |, [1, 2, 3]  |
              |  |               |  |                |             |
              |  |   0,  2/5, 1  |  |  0,  0,  -4/5  |             |
              -- +-             -+  +-              -+            --

The lower triangular matrix L is the first element und the upper triangular matrix U is the second element of the list `LU`. The product of these two matrices is equal to the input matrix `A`:

    >> L * U

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

#### Example 2

An LU-decomposition of the 3 x 2 matrix:

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

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

gives a 3 x 3 lower triangular matrix and a 2 x 3 upper triangular matrix:

    >> [L, U, pivlist] := linalg::factorLU(A)

                 -- +-              -+  +-        -+            --
                 |  |   1,    0,  0  |  |  2,  -3  |             |
                 |  |                |  |          |             |
                 |  |  1/2,   1,  0  |, |  0, 7/2  |, [1, 2, 3]  |
                 |  |                |  |          |             |
                 |  |   1,  12/7, 1  |  |  0,  0   |             |
                 -- +-              -+  +-        -+            --

    >> L * U

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

#### Example 3

To compute the LU-decomposition of the matrix:

    >> A := matrix([[1, 2, -1], [0, 0, 3], [0, 2, -1]])

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

one row interchange is needed, and we therefore get a non-trivial permutation list:

    >> [L, U, pivlist] := linalg::factorLU(A)

                  -- +-         -+  +-          -+            --
                  |  |  1, 0, 0  |  |  1, 2, -1  |             |
                  |  |           |  |            |             |
                  |  |  0, 1, 0  |, |  0, 2, -1  |, [1, 3, 2]  |
                  |  |           |  |            |             |
                  |  |  0, 0, 1  |  |  0, 0,  3  |             |
                  -- +-         -+  +-          -+            --

The corresponding permutation matrix is the following:

    >> P := linalg::swapRow(matrix::identity(3), 3, 2)

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

Hence, we have a decomposition of A into the product of the three matrices P^(-1), L and U as follows:

    >> P^(-1) * L * U

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

#### Example 4

You may compute an LU-decomposition of a matrix with symbolic components, such as:

    >> delete a, b, c, d:
       A := matrix([[a, b], [c, d]])

                                    +-      -+
                                    |  a, b  |
                                    |        |
                                    |  c, d  |
                                    +-      -+

The diagonal entries of the matrix U are the pivot elements used during the computation. They must be non-zero, if the inverse of U is needed:

    >> [L, U, pivlist] := linalg::factorLU(A)

                    -- +-      -+  +-            -+         --
                    |  |  1, 0  |  |  a,    b     |          |
                    |  |        |  |              |          |
                    |  |  c     |, |         b c  |, [1, 2]  |
                    |  |  -, 1  |  |  0, d - ---  |          |
                    |  |  a     |  |          a   |          |
                    -- +-      -+  +-            -+         --

For example, if we use this decomposition to solve the linear system A*x=b for arbitrary vectors b=(b[1],b[2])^T, then the following result is only correct for a <> 0 and d-(b*c)/d <> 0:

    >> delete b1, b2:
       linalg::matlinsolveLU(L, U, matrix([b1, b2]))

                            +-                      -+
                            |         /      c b1 \  |
                            |       b | b2 - ---- |  |
                            |         \       a   /  |
                            |  b1 - ---------------  |
                            |               b c      |
                            |           d - ---      |
                            |                a       |
                            |  --------------------  |
                            |           a            |
                            |                        |
                            |             c b1       |
                            |        b2 - ----       |
                            |              a         |
                            |        ---------       |
                            |             b c        |
                            |         d - ---        |
                            |              a         |
                            +-                      -+

#### Background

- The following algorithm for solving the system A*x=b with a nonsingular matrix A uses LU-decomposition:
  1. Compute a LU-decomposition of A: A = L*U.
  2. Solve y=LMhat (-1)*b by forward substitution.
  3. Solve x=R^(-1)*y by backward substitution.
- The LU-decomposition of a matrix A is useful for solving several systems of linear equations A*x=b with the same coefficient matrix A and several right-hand side vectors b, because then step one of the algorithm above needs to be done only once.

#### Changes

- new function

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