numeric::singularvectors
-- numerical singular value decomposition of a matrix
Introductionnumeric::singularvectors(A, ..) returns
numerical singular values and singular vectors of the matrix
A.
Call(s)numeric::singularvectors(A <, NoLeftVectors> <, NoRightVectors> <, NoErrors>)
ParametersA |
- | a numerical matrix of domain type DOM_ARRAY or of category Cat::Matrix |
OptionsNoLeftVectors |
- | suppresses the computation of left singular vectors |
NoRightVectors |
- | suppresses the computation of right singular vectors |
NoErrors |
- | suppresses the computation of error estimates |
Returnsa 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
EffectsThe function is sensitive to the environment variable DIGITS, which determines the
numerical working precision.
Related
Functionslinalg::eigenvalues, linalg::eigenvectors,
numeric::eigenvalues, numeric::eigenvectors,
numeric::singularvalues,
numeric::spectralradius
DetailsA 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).[U,d,V,resU,resV] returned by
numeric::singularvectors corresponds to the singular data
of an m x n matrix A as described below.
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].numeric::singularvectors as an array of floating point
numbers.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.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.
numeric::singularvectors are identical to those computed
by numeric::singularvalues.[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: NoLeftVectorsU and the
corresponding residues resU. The return values for these
data are NIL.
Option: NoRightVectorsV and the
corresponding residues resV. The return values for these
data are NIL.
Option: NoErrorsresU and
resV. The return values for these data are
NIL.
Example
1Numerical 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
2We 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:
Cat::Matrix objects now uses the method
"expr" of the matrix domain.numeric::sort from small values to
large values. Matrices of singular vectors are now returned as
arrays.