## `numeric::eigenvectors` -- numerical eigenvalues and eigenvectors of a matrix

#### Introduction

`numeric::eigenvectors``(A, ..)` returns numerical eigenvalues and eigenvectors of the matrix `A`.

#### Call(s)

`numeric::eigenvectors(A <, NoErrors>)`

#### Parameters

`A` | - | a numerical square matrix of domain type `DOM_ARRAY` or of category `Cat::Matrix`

#### Options

`NoErrors` | - | suppresses the computation of error estimates

#### Returns

a list `[d, X, res]`. The sorted list `d=[d[1],d[2],..]` contains the numerical eigenvalue. The array `X` is the matrix of eigenvectors, i.e., the i-th column of `X` is the eigenvector associated with the eigenvalue `d[i]`. The list of residues `res=[res[1],res[2],..]` provides error estimates for the numerical eigenvalues.

#### Side Effects

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

#### Related Functions

`linalg::eigenvalues`, `linalg::eigenvectors`, `numeric::eigenvalues`, `numeric::singularvalues`, `numeric::singularvectors`, `numeric::spectralradius`

#### Details

- All entries of the matrix must be numerical. Numerical expressions such as `exp(PI), sqrt(2)` etc. are accepted and converted to floats. Non-numerical symbolic entries lead to an error.
- The eigenvalues are sorted by `numeric::sort`.

- Eigenvalues are approximated with an absolute precision of 10^(-DIGITS)*r, where r is the spectral radius of `A` (i.e., r is the maximum of the absolute values of the eigenvalues). Consequently, large eigenvalues should be computed correctly to `DIGITS` decimal places. The numerical approximations of the small eigenvalues are less accurate.

- The array `X` is the matrix of eigenvectors, i.e., the i-th column of `X` is a numerical eigenvector corresponding to the eigenvalue `d[i]`. Each column is either zero or normalized to the Euclidean length 1.0.
- For matrices with multiple eigenvalues and an insufficient number of eigenvectors some of the eigenvectors may coincide or may be zero, i.e., `X` is not necessarily invertible.
- The list of residues res=[res[1],res[2],..] provides some control over the quality of the numerical spectral data. The residues are given by

      res[i]=norm(A*x[i]-d[i]*x[i], 2),

  where x[i] is the normalized eigenvector (the i-th column of `X`) associated with the numerical eigenvalue d[i]. For Hermitean matrices res[i] provides an upper bound for the absolute error of d[i].
- With the option NoErrors the computation of the residues is suppressed, the returned value is `NIL`.

- Matrices `A` of a matrix domain such as `Dom::Matrix``(..)` or `Dom::SquareMatrix``(..)` are internally converted to arrays over expressions via `A::dom::expr(A)`. Note that `linalg::eigenvectors` must be used, when the eigenvalues/vectors are to be computed over the component domain. Cf. example 3.

- For a numerical algorithm it is not possible to distinguish between badly separated distinct eigenvalues and multiple eigenvalues. For this reason `numeric::eigenvectors` and `linalg::eigenvectors` use different return formats: the latter can provide information on the multiplicity of eigenvalues due to its internal exact arithmetic.

- Use `numeric::eigenvalues`, if only eigenvalues are to be computed.

#### Example 1

We compute the spectral data of the 2x2 Hilbert matrix:

    >> A := linalg::hilbert(2)

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

    >> [d, X, res] := numeric::eigenvectors(A):

The eigenvalues:

    >> d

                           [0.06574145409, 1.267591879]

The eigenvectors:

    >> X

                        +-                              -+
                        |   0.4718579255, -0.8816745988  |
                        |                                |
                        |  -0.8816745988, -0.4718579255  |
                        +-                              -+

Hilbert matrices are Hermitean, i.e., computing the spectral data is a numerically stable process. This is confirmed by the small residues:

    >> res

                        [3.965706585e-20, 5.421010863e-20]

For further processing we convert the data to matrices of the domain `Dom::Matrix``()`:

    >> M := Dom::Matrix(): X := M(X): d := M(2, 2, d, Diagonal):

We reconstruct the matrix from its spectral data:

    >> X*d/X

                              +-                   -+
                              |  1.0,      0.5      |
                              |                     |
                              |  0.5, 0.3333333333  |
                              +-                   -+

We extract an eigenvector from the matrix `X`:

    >> eigenvector1 := X::dom::col(X, 1)

                                +-               -+
                                |   0.4718579255  |
                                |                 |
                                |  -0.8816745988  |
                                +-               -+

    >> delete A, d, X, res, M, eigenvector1:

#### Example 2

We demonstrate some numerically ill-conditioned cases. The following matrix has only one eigenvector and cannot be diagonalized. Only one numerical eigenvector is found:

    >> A := array(1..2, 1..2, [[5, -1], [4, 1]]):

    >> numeric::eigenvectors(A)

          --             +-                 -+                         --
          |              |  0, 0.4472135955  |                          |
          |  [3.0, 3.0], |                   |, [0.0, 1.084202173e-19]  |
          |              |  0,  0.894427191  |                          |
          --             +-                 -+                         --

Dividing `A` by 3 leads to slightly increased internal round-off. This time two badly separated eigenvalues are computed. The two corresponding numerical eigenvectors are almost collinear. Both represent the same exact eigenvector:

    >> B := map(A, _divide, 3):

    >> numeric::eigenvectors(B)

          --                      +-                             -+
          |                       |  0.4472135954, -0.4472135957  |
          |  [0.9999999997, 1.0], |                               |,
          |                       |   0.894427191, -0.8944271909  |
          --                      +-                             -+

                                                --
                                                 |
             [1.084202173e-19, 6.060874398e-20]  |
                                                 |
                                                --

    >> delete A:

#### Example 3

The following matrix has domain components:

    >> A := Dom::Matrix(Dom::IntegerMod(7))([[6, -1], [0,  3]])

                              +-                  -+
                              |  6 mod 7, 6 mod 7  |
                              |                    |
                              |  0 mod 7, 3 mod 7  |
                              +-                  -+

Note that `numeric::eigenvectors` computes the spectral data of the following matrix:

    >> A::dom::expr(A)

                                    +-      -+
                                    |  6, 6  |
                                    |        |
                                    |  0, 3  |
                                    +-      -+

    >> numeric::eigenvectors(A, NoErrors)

                  --             +-                   -+      --
                  |              |  -0.894427191, 1.0  |       |
                  |  [3.0, 6.0], |                     |, NIL  |
                  |              |  0.4472135955, 0.0  |       |
                  --             +-                   -+      --

If the spectral data are to be computed over the component domain `Dom::IntegerMod``(7)`, then `linalg::eigenvectors` should be used:

    >> linalg::eigenvectors(A)

          -- --             -- +-         -+ -- --
          |  |              |  |  1 mod 7  |  |  |
          |  |  6 mod 7, 1, |  |           |  |  |,
          |  |              |  |  0 mod 7  |  |  |
          -- --             -- +-         -+ -- --

             --             -- +-         -+ -- -- --
             |              |  |  5 mod 7  |  |  |  |
             |  3 mod 7, 1, |  |           |  |  |  |
             |              |  |  1 mod 7  |  |  |  |
             --             -- +-         -+ -- -- --

    >> delete A:

#### Background

- The function implements standard numerical algorithms from the Handbook of Automatic Computation by Wilkinson and Reinsch.

#### Changes

- Conversion of `Cat::Matrix` objects now uses the method `"expr"` of the matrix domain. Triangular matrices are now processed numerically.
- The eigenvector matrix is now returned as an array.

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