## `linalg::wiedemann` -- solving linear systems by Wiedemann's algorithm

#### Introduction

`linalg::wiedemann``(A, b, mult ...)` tries to find a vector x that satisfies the equation A*x = b by using Wiedemann's algorithm.

#### Call(s)

`linalg::wiedemann(A, b <, mult>)`
`linalg::wiedemann(A, b <, mult>, prob)`

#### Parameters

`A` | - | an n x n matrix of a domain of category `Cat::Matrix`
`b` | - | an n-dimensional column vector, i.e., an n x 1 matrix of a domain of category `Cat::Matrix`
`mult` | - | a matrix-vector multiplication method: function or functional expression; default: `_mult`
`prob` | - | `TRUE` or `FALSE` (default: `TRUE`)

#### Returns

either the list `[x, TRUE]` if a solution for the system A*x=b has been found, or the list `[x, FALSE]` if a non-zero solution for the corresponding homogeneous system A*x=0 has been found, or the value `FAIL` (see below).

#### Related Functions

`linalg::matlinsolve`, `linalg::vandermondeSolve`

#### Details

- The parameter `mult` must be a function such that the result of `mult(A,y)` equals A*y for every n-dimensional column vector y. The parameter `y` is of the same domain type as `A`. The argument `mult` does not need to handle other types of parameters, nor does it need to handle other matrices than `A`.
- `linalg::wiedemann` uses a probabilistic algorithm. For a deterministic variant enter `FALSE` for the optional parameter `prob`.
- If the system A*x=b does not have a solution, then `linalg::wiedemann` returns `FAIL`.
- If the system A*x=b has more than one solution, then a random one is returned.
- Due to the probabilistic nature of Wiedemann's algorithm, the computation may fail with small probability. In this case `FAIL` is returned. If the deterministic variant is chosen, then the algorithm may be slower for a small number of matrices.
- The vector `b` must be defined over the component ring of `A`.
- The coefficient ring of `A` must be a field, i.e., a domain of category `Cat::Field`.
- It is recommended to use `linalg::wiedemann` only if `mult` uses significantly less than O(n^2) field operations.

#### Example 1

We define a matrix and a column vector over the finite field with 29 elements:

    >> MatZ29 := Dom::Matrix(Dom::IntegerMod(29)):
       A := MatZ29([[1, 2, 3], [4, 7, 8], [9, 12, 17]]);
       b := MatZ29([1, 2, 3])

                       +-                                -+
                       |  1 mod 29,  2 mod 29,  3 mod 29  |
                       |                                  |
                       |  4 mod 29,  7 mod 29,  8 mod 29  |
                       |                                  |
                       |  9 mod 29, 12 mod 29, 17 mod 29  |
                       +-                                -+

                                  +-          -+
                                  |  1 mod 29  |
                                  |            |
                                  |  2 mod 29  |
                                  |            |
                                  |  3 mod 29  |
                                  +-          -+

Since `A` does not have a special form that would allow a fast matrix-vector multiplication, we simply use `_mult`. Wiedemann's algorithm works in this case, although it is less efficient than Gaussian elimination:

    >> linalg::wiedemann(A, b, _mult)

                            -- +-           -+       --
                            |  |  24 mod 29  |        |
                            |  |             |        |
                            |  |  21 mod 29  |, TRUE  |
                            |  |             |        |
                            |  |  17 mod 29  |        |
                            -- +-           -+       --

#### Example 2

Now let us define another matrix that has a special form:

    >> MatZ29 := Dom::Matrix(Dom::IntegerMod(29)):
       A := MatZ29([[1, 0, 0], [0, 1, 2], [0, 0, 1]]);
       b := MatZ29(3, 1, [1, 2, 3]):

                        +-                              -+
                        |  1 mod 29, 0 mod 29, 0 mod 29  |
                        |                                |
                        |  0 mod 29, 1 mod 29, 2 mod 29  |
                        |                                |
                        |  0 mod 29, 0 mod 29, 1 mod 29  |
                        +-                              -+

For this particular matrix, it is easy to define an efficient multiplication method:

    >> mult := proc(dummy, y)
       begin
           y[2]:=y[2]+2*y[3];
           y
       end:
       linalg::wiedemann(A, b, mult)

                            -- +-           -+       --
                            |  |   1 mod 29  |        |
                            |  |             |        |
                            |  |  25 mod 29  |, TRUE  |
                            |  |             |        |
                            |  |   3 mod 29  |        |
                            -- +-           -+       --

#### Background

- The expected running time for the probabilistic algorithm is O(n^2 + n*M), and the running time for the deterministic variant is O(n^2*M) in the worst case, but only O(n^2 + n*M) on average. Here, M is the number of field operations that the matrix-vector multiplication routine `mult` uses.
- The basic idea of the algorithm is to solve a linear system A*x = b by finding the minimal polynomial f(y) that solves f(A)*b = 0. If the constant coefficient c=f(0) is nonzero and g(y):= f(y) - c, the equality g(A)*b = -cb implies that x = (-1/c) ((g/y)(A)) is the solution.

  The polynomial f is found by looking for the minimal polynomial h satisfying u h(A) b = 0 for some randomly chosen row vector u. This may yield h ≠f in unlucky cases, but in general the probability for this is small.

- Reference: Douglas Wiedemann: Solving Sparse Linear equations over Finite Fields, IEEE Transactions on Information Theory, vol. 32, no.1, Jan. 1986.

#### Changes

- new function

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