SingularValueDecomp

SingularValueDecomp(a, i, j, j2)

SingularValueDecomp computes the singular value decomposition of a matrix. Singular value decomposition is often used with sets of equations or matrices that are singular or ill-conditioned (that is, very close to singular). It factors a matrix «a», indexed by «i» and «j», with IndexLength(i) >= IndexLength(i), into three matrices, U, W, and V, such that:

a = U . W . VT

where U and V are orthogonal matrices and W is a diagonal matrix. U is dimensioned by «i» and «j», W by «j» and «j2», and V by «j» and «j2». In Analytica notation:

`Variable A := Sum(Sum(U*W, J)*Transpose(V, J, J2), J2)`

The index «j2» must be the same size as «j» and is used to index the resulting W and V arrays. SingularValueDecomp returns an array of three elements indexed by a special system index named `SvdIndex` with each element, U, W, and V, being a reference to the corresponding array.

Use the # (dereference) operator to obtain the matrix value from each reference, as in:

`Index J2 := CopyIndex(J)`
`Variable SvdResult := SingularValueDecomp(A, I, J, J2)`
`Variable U := #SvdResult[SvdIndex = 'U']`
`Variable W := #SvdResult[SvdIndex = 'W']`
`Variable V := #SvdResult[SvdIndex = 'V']`

Matrix inverse

The inverse of a square matrix A, in Analytica syntax, is ``` ```

``` Var Winv := If J=J2 Then 1/W Else 0; Transpose(Sum(Sum(U*Winv, J)*Transpose(V, J, J2), J2),I,J) ```

Singular value decomposition can be used for matrix inverse when the matrix A is ill-conditioned, in which case the Invert function may encounter numeric instabilities. When the matrix is ill-conditioned (the Determinant is very close to zero), then some of the elements of the diagonal of `W` will be very close to zero. To avoid the numerical instabilities, the diagonal entries corresponding to the very small `W` can be replaced with 0 in `Winv`:

``` ```

``` Var Winv := If J=J2 And Abs(W)>1e-4 Then 1/W Else 0; Transpose(Sum(Sum(U*Winv, J)*Transpose(V, J, J2), J2),I,J) ```

To make this convenient to use, you can introduce a new User-Defined Function as follows:

Function `MatInvert( A : [I,J] ; I,J : Index )`
Definition:``` Index J2 := J; Local svd := SingularValueDecomp(A,I,J,J2); Local U := #svd[SvdIndex='U']; Local W := #svd[SvdIndex='W']; Local Winv := if J=J2 And W>1e-5 then 1/W else 0; Local V := #svd[SvdIndex='V']; Transpose(Sum(Sum(U*Winv, J)*Transpose(V, J, J2), J2),I,J)```
``` ```

You can then use `MatInvert(A,I,J)` in place of `Invert(A,I,J)`.