numeric::singularvectors
-- numerical singular value decomposition of a matrixnumeric::singularvectors
(A, ..)
returns
numerical singular values and singular vectors of the matrix
A
.
numeric::singularvectors(A <, NoLeftVectors> <, NoRightVectors> <, NoErrors>)
A |
- | a numerical matrix of domain type DOM_ARRAY or of category Cat::Matrix |
NoLeftVectors |
- | suppresses the computation of left singular vectors |
NoRightVectors |
- | suppresses the computation of right singular vectors |
NoErrors |
- | suppresses the computation of error estimates |
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.
The function is sensitive to the environment variable DIGITS
, which determines the
numerical working precision.
linalg::eigenvalues
, linalg::eigenvectors
,
numeric::eigenvalues
, numeric::eigenvectors
,
numeric::singularvalues
,
numeric::spectralradius
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)
.[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.U
and the
corresponding residues resU
. The return values for these
data are NIL
.V
and the
corresponding residues resV
. The return values for these
data are NIL
.resU
and
resV
. The return values for these data are
NIL
.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:
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:
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.