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:

[math]\displaystyle{ a = U W V^T }[/math]

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.

To get each component in a separate global variable, use:

Index J2 ::= CopyIndex(J)

Variable U ::= ( , W, V ) := SingularValueDecomp(A, I, J, J2)

Variable W ::= ComputedBy(U)

Variable V ::= ComputedBy(U)

An less convenient but alternative approach is to 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']

Note that SingularValueDecomp detects whether you are capturing multiple return values or not and changes its return value to one or the other behaviors shown above. In the first case, returning each matrix as a separate return value, and in the second case returning a single return value indexed by SvdIndex.

If you need the matrices in local values, use:

Local (u, w, v) := SingularValueDecomp(A, I, J, J2);

Matrix inverse

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

Local 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:

Local 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 (u, w, v) := SingularValueDecomp(a, I, J, J2);
Local wInv := if J=J2 And w>1e-5 then 1/w else 0;
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).

SingularValueDecomp

SingularValueDecomp(a, i, j, j2)

Matrix inverse

See Also