Complex SVD

[mbibby]

I have read all the previous contributions to this forum on the subject of complex svd. However, I have not seen anything that suggests that anyone has actually successfully converted the NR svdcmp routine to complex.

Has anyone done this? If so, would they be willing to share it?

Thanks.

M.

[evgeny]

Hi,

If the matrix is not too large and efficiency and accuracy is not a big problem you can process without rewriting of SVD from NR:

That is, given complex matrix A
1) Compute a real matrix Conjugate(Transpose(A)). A
2) Find the eigenvalues of above matrix by any method.
The square roots of these will be the singular values.
3) The (normalized) eigenvectors of above matrix will be the columns of V.
4) The normalized columns of A.V will be the columns of U.

Evgeny

[evgeny]

Hi,

If a memory and time is not critical, You don't need change svd from NR at all.

Indeed, let A = Ar+Ai*I, where I*I=-1.
Suppose, that U=Ur+Ui*I and V=Vr+Vi*I.
Than complex SVD, defined as A=Conjugate[U].W.V, may be written by block matrix form:
{{Ar,Ai},{-Ai,Ar}} =

{
{Transpose[Ur],Transpose[Ui]},
{Transpose[Ui],-Transpose[Ur]}
} .

{
{W,0},
{0,W}
} .

{
{Vr,Vi},
{Vi,-Vr}
}


Algorithm:
Step1: form an extended matrix Anew = {{Ar,Ai},{-Ai,Ar}} from a given matrix A=Ar+Ai*I
Step2: do svd for above extended matrix
Step3: extract from receved U and V matrices needed Ur,Ui,Vr,Vi parts to form U and V
Step4: take into attention order of vectors and their non uniqness

Evgeny.

[hjsmithh]

I implemented your second method, but I don't know how to extract U and V parts from the computed U and V. W was easy since the W vector is real and sorted, I just took every other entry.

-Harry

[evgeny]

Hi,

I checked this method via Mathematica only and didnt find any problem to exract vectors.

Of course, first, You need sort W values inconjuction with U and V matrices rows (or columns, depend on C of Fortran code style).

Then drop each second entry as for W as well as for U and V matrices (due to doubling of vectors and values).

Finally take corresponding entries of Re[V] and Im[V] and form complex vectors.

A problem is only non uniqness SVD (for example, up to a sign or phase).

Evgeny, Haifa.

[forumposter]

Hi evgeny,

I tried your 2nd method in matlab and was confused by the result. Here's what I got in matlab:

>> A

A =

58.0000 +33.0000i 68.0000 +52.0000i 77.0000 +41.0000i
27.0000 + 8.0000i 31.0000 +21.0000i 47.0000 +21.0000i
6.0000 +66.0000i 19.0000 +49.0000i 26.0000 +25.0000i

>> AA=[[real(A), -imag(A)];[imag(A), real(A)]]

AA =

58 68 77 -33 -52 -41
27 31 47 -8 -21 -21
6 19 26 -66 -49 -25
33 52 41 58 68 77
8 21 21 27 31 47
66 49 25 6 19 26

>> [U, S, V]=svd(A)

U =

-0.6404 - 0.4529i -0.0557 - 0.3126i -0.3408 - 0.4097i
-0.3312 - 0.1907i -0.1878 - 0.3888i 0.6102 + 0.5432i
-0.1577 - 0.4624i -0.3118 + 0.7845i -0.0215 + 0.2193i


S =

176.8697 0 0
0 37.1941 0
0 0 4.9700


V =

-0.5316 0.7576 0.3788
-0.6051 + 0.0141i -0.0407 + 0.1000i -0.7678 - 0.1802i
-0.5830 - 0.1057i -0.6075 - 0.2130i 0.3969 + 0.2776i

>> [UU, SS, VV]=svd(AA)

UU =

-0.6404 0.4529 -0.0583 -0.3121 -0.4113 0.3389
-0.3312 0.1907 0.0384 -0.4301 0.5460 -0.6077
-0.1577 0.4624 0.5688 0.6238 0.2192 0.0225
-0.4529 -0.6404 0.3121 -0.0583 0.3389 0.4113
-0.1907 -0.3312 0.4301 0.0384 -0.6077 -0.5460
-0.4624 -0.1577 -0.6238 0.5688 0.0225 -0.2192


SS =

176.8697 0 0 0 0 0
0 176.8697 0 0 0 0
0 0 37.1941 0 0 0
0 0 0 37.1941 0 0
0 0 0 0 4.9700 0
0 0 0 0 0 4.9700


VV =

-0.5316 0.0000 -0.7087 0.2676 0.0018 -0.3788
-0.6051 -0.0141 0.0734 0.0791 -0.1838 0.7670
-0.5830 0.1057 0.4931 -0.4138 0.2794 -0.3956
-0.0000 -0.5316 -0.2676 -0.7087 -0.3788 -0.0018
0.0141 -0.6051 -0.0791 0.0734 0.7670 0.1838
-0.1057 -0.5830 0.4138 0.4931 -0.3956 -0.2794

A is a randomly generated 3x3 complex matrix. AA is the extended matrix from A. You can see from the result it seems impossible to extract elements from UU and VV to form U and V because some elements in U and V don't exist at all in UU and VV. Did I do anything wrong? At first I thought I messed up the row/col order while extending A to AA. But I tried different combinations, all got similar results. I would really appreciate it if you can help me out.

Thx very much!

[hjsmithh]

Here is my Visual C# .NET code that works for me. This is part of my multiple precision complex matrix program XZCalc.

public static void SVdecompComplex(MultiZ a, MultiZ w, MultiZ v)
// Singular Value Decomposition, a = u * w * ConjTran(v), a replaced with u
// the diagonal matrix is output as a row vector, a.Rs must be >= a.Cs
// if smaller, a should be filled with zero rows
// this code is adapted from Complex SVD - Numerical Recipes Forum
{
int n = a.Rs;
int m = a.Cs;
MultiZ aNew = new MultiZ(Calc.fMC, n+n, m+m);
MultiZ wNew = new MultiZ(Calc.fMC, 1, m+m);
MultiZ vNew = new MultiZ(Calc.fMC, n+n, m+m);
for (int r = 0; r < n; r++)
for (int c = 0; c < m; c++)
{
aNew.E[r, c] &= (MultiF)a.E[r, c];
aNew.E[r, c+m] &= a.E[r, c].I;
aNew.E[n+r, c] &= -a.E[r, c].I;
aNew.E[n+r, m+c] &= (MultiF)a.E[r, c];
}
SVdecomp(aNew, wNew, vNew);
for (int r = 0; r < n; r++)
for (int c = 0; c < m; c++)
{
a.E[r, c] &= aNew.E[r, c+c];
a.E[r, c].I &= -aNew.E[r+n, c+c];
v.E[r, c] &= vNew.E[r, c+c];
v.E[r, c].I &= -vNew.E[r+n, c+c];
w.E[0, c] &= wNew.E[0, c+c];
}
} // SVdecompComplex

-Harry

[evgeny]

Hi,

First it is best do calculation (print) with double accuracy (not up to 4 digits).

Second, please take into attention that
first & second columns in VV give right V vector,
fifth & sixth columns in VV give slightly different values with V and
third & forth columns give very different values with V.
Why?

An answer is simple. SVD is not uniqness.
In your V matrix after MATLAB calculation, a first vector component in each column is calculated at assumption that it's real part equals null (zero phase). In VV matrix, calculated by proposed method, this isn't true.

A discrepancy especially large for a first component: 0.2676-0.7087*I in 3-4 columns and less for a first component: -0.3788+0.0018*I in 5-6 columns as "phase" is larger for a first one.

You can also check that for example components in 3-4 columns of VV matrix have same length (calculated as abs value of complex number) as in V matrix but are "rotated" at the same complex constant which is coincides with phase of a first component. Best to check this fact on components from 4,5,6 rows and 2,3 columns of VV matrix.

A choosing of the phase may be done if it is required by Your physical problem.

Evgeny. Haifa

[hjsmithh]

I changed my mind. I never could get the algorithm discussed here to work in all case. I have switched to ACM Algorithm 358 for Singular Value Decomposition of a complex matrix. See: http://www.scs.fsu.edu/~burkardt/f77...s358/toms358.f and http://www.scs.fsu.edu/~burkardt/f77.../toms358_prb.f for the Fortran code.

[evgeny]

You are right on 100 percent!

The source is good.

By the way, do you know a C-code for this algorithm?


Evgeny.

[hjsmithh]

I have published my latest version on my calculator program XZCalc that uses ACM Algorithm 358 for SVD of a complex matrix. The executable program and all source files are included in the install file XZCalc32L.exe. It can be downloaded from http://www.geocities.com/hjsmithh/download.html#XZCalc . This program is free software.

[hjsmithh]

I now have a working C++ version of CSVD that can be downloaded from http://www.geocities.com/hjsmithh/download.html#Algo358 The source and object files are included.

[hjsmithh]

If that is not enough, I now have a working C version of CSVD that can be downloaded from http://www.geocities.com/hjsmithh/download.html#Algo358 The source and object files are included.

I also used this to write a C program to compute the Eigenvalues and Eigenvectors of a square complex matrix. This can be downloaded from http://www.geocities.com/hjsmithh/download.html#Eigen The source and object files are included.