## `numeric::singularvectors` -- numerical singular value decomposition of a matrix

#### Introduction

`numeric::singularvectors``(A, ..)` returns numerical singular values and singular vectors of the matrix `A`.

#### Call(s)

`numeric::singularvectors(A <, NoLeftVectors> <, NoRightVectors> <, NoErrors>)`

#### Parameters

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

#### Options

`NoLeftVectors` | - | suppresses the computation of left singular vectors
`NoRightVectors` | - | suppresses the computation of right singular vectors
`NoErrors` | - | suppresses the computation of error estimates

#### Returns

a list `[U,d,V,resU,resV]`. `U` is a unitary square float matrix of domain type `DOM_ARRAY`, whose columns are left singular vectors. `d` is a list of singular float values. `V` is a unitary square float matrix of domain type `DOM_ARRAY`, whose columns are right singular vectors. The lists of float residues `resU` and `resV` provide error estimates for the numerical data.

#### 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::eigenvectors`, `numeric::singularvalues`, `numeric::spectralradius`

#### Details

- All entries of `A` 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.
- `Cat::Matrix` objects, i.e., 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)`.
- The list `[U,d,V,resU,resV]` returned by `numeric::singularvectors` corresponds to the singular data of an m x n matrix A as described below.
- Let V^H denote the Hermitean transpose of the matrix V, i.e., the complex conjugate of the transpose. The singular value decomposition of an m x n matrix A is a factorization A=U*D*V^H.
- D is a sparse m x n matrix with real nonnegative ``diagonal'' entries D[i,i]=d[i], i=1,..,p, p=min(m,n).

  The list d=[d[1],..,d[p]] returned by `numeric::singularvectors` are the ``singular values'' of A. They are sorted by `numeric::sort`, i.e., 0.0 <= d[1] <= .. <= d[p].
- U is a unitary m x m matrix. Its i-th column is an eigenvector of A*A^H associated with the eigenvalue d[i]^2 (d[i]=0 for i>p). These are the ``left singular vectors'' of A. They are returned by `numeric::singularvectors` as an array of floating point numbers.
- V is a unitary n x n matrix. Its i-th column is an eigenvector of A^H*A associated with the eigenvalue d[i]^2 (d[i]=0 for i>p). These are the ``right singular vectors'' of A. They are returned by `numeric::singularvectors` as an array of floating point numbers. It is normalized such that its diagonal entries are real and nonnegative.
- `resU=[resU[1],..,resU[m]]` is a list of float residues associated with the left singular vectors:

      resU[i] = <A^H*u.i, A^H*u.i> - d[i]^2, i=1,..,m.

  Here u.i is the (normalized) i-th column of `U`, <.,.> is the usual complex Euclidean scalar product and d[i]=0 for p<i<=m.
- `resV=[resV[1],..,resV[n]]` is a list of float residues associated with the right singular vectors:

      resV[i] = <A*v.i, A*v.i> - d[i]^2, i=1,..,m.

  Here v.i is the (normalized) i-th column of `V`, d[i]=0 for p<i<=n.
- The residues `resU,resV` vanish for exact singular data U,d,V. Their size indicate the quality of the numerical data `U,d,V`.

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

- The singular values computed by `numeric::singularvectors` are identical to those computed by `numeric::singularvalues`.
- Singular data may also be computed via

      [d2, U, errU]:=numeric::eigenvectors(A*A^H);

  or

      [d2, V, errV]:=numeric::eigenvectors(A^H*A);

  respectively. The list `d2` is related to the singular values by

      d2 = [0, .. , 0, d[1]^2,d[2]^2, .. ,d[p]^2] .

  The use of `numeric::singularvectors` avoids the costs of the matrix multiplication. Further, the eigenvector routine requires about twice as many `DIGITS` to compute the data associated with small singular values with the same precision as `numeric::singularvectors`. Also note that the normalization of `U` and `V` may be different.

#### Option: NoLeftVectors

- If only right singular vectors are required, then this option may be used to suppress the computation of `U` and the corresponding residues `resU`. The return values for these data are `NIL`.

#### Option: NoRightVectors

- If only left singular vectors are required, then this option may be used to suppress the computation of `V` and the corresponding residues `resV`. The return values for these data are `NIL`.

#### Option: NoErrors

- If no error estimates are required, then this option may be used to suppress the computation of the residues `resU` and `resV`. The return values for these data are `NIL`.

#### Example 1

Numerical expressions are converted to floats:

    >> DIGITS := 5:

    >> A :=array(1..3, 1..2, [[1, PI], [2, 3], [3, exp(sqrt(2))]]):

    >> [U, d, V, resU, resV] := numeric::singularvectors(A):

The singular data are:

    >> U, d, V

          +-                             -+
          |  -0.88078, 0.45729,  0.12293  |
          |                               |
          |   0.14947, 0.51483, -0.84417  |, [0.89905, 6.9986],
          |                               |
          |   0.44932, 0.72515,  0.5218   |
          +-                             -+

             +-                   -+
             |   0.85215,  0.5233  |
             |                     |
             |  -0.52331, 0.85215  |
             +-                   -+

The residues indicate that these results are not severely affected by round-off within the working precision of 5 digits:

    >> resU, resV

            [2.0954e-9, 2.9802e-8, 1.7347e-18], [3.4925e-9, 7.4506e-8]

    >> delete DIGITS, A, U, d, V, resU, resV:

#### Example 2

We demonstrate how to reconstruct a matrix from its singular data:

    >> DIGITS := 3: A := array(1..2, 1..3, [[1.0, I, PI], [2, 3, I]]):

    >> [U, d, V, resU, resV] := numeric::singularvectors(A, NoErrors):

For convenience, the matrix domain `Dom::Matrix``()` is used to process the matrices:

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

                     +-                                     -+
                     |  - 0.789 + 0.336 I, 0.511 + 0.0665 I  |
                     |                                       |
                     |   0.487 - 0.168 I,    0.84 + 0.17 I   |
                     +-                                     -+

A ``diagonal'' matrix is built from the singular values:

    >> DD := M(2, 3, d, Diagonal)

                                +-              -+
                                |  3.27,  0,  0  |
                                |                |
                                |    0,  3.9, 0  |
                                +-              -+

    >> V := M(V)

            +-                                                      -+
            |       0.0568,      0.562 + 0.104 I, - 0.681 - 0.454 I  |
            |                                                        |
            |   0.55 + 0.0871 I,      0.663,       0.454 + 0.208 I   |
            |                                                        |
            |  - 0.81 + 0.174 I, 0.455 - 0.162 I,       0.283        |
            +-                                                      -+

We use the methods `conjugate` and `transpose` of the matrix domain to compute the Hermitean transpose of `V` and reconstruct `A` up to numerical round-off:

    >> VH := M::conjugate(M::transpose(V)):

    >> U*DD*VH

            +-                                                       -+
            |  1.0 + 2.62e-10 I,       1.0 I,      3.14 + 4.66e-10 I  |
            |                                                         |
            |  2.0 + 5.02e-10 I, 3.0 - 1.16e-10 I, - 3.73e-9 + 1.0 I  |
            +-                                                       -+

    >> delete DIGITS, A, U, d, V, resU, resV, M, DD, VH:

#### Changes

- Conversion of `Cat::Matrix` objects now uses the method `"expr"` of the matrix domain.
- The singular values are now sorted by `numeric::sort` from small values to large values. Matrices of singular vectors are now returned as arrays.

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