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.

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

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.

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

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.

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!

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

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

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_src/toms358/toms358.f and http://www.scs.fsu.edu/~burkardt/f77_src/toms358/toms358_prb.f for the Fortran code.

You are right on 100 percent!

The source is good.

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


Evgeny.

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.

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.

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.

Hi Harry (hjsmithh),

Thank you for sharing the C# CSVD code. While testing the code, I found out that it does not work for some cases. After SVD, the decomposed matrices should reproduce the same original matrix (i.e., [A] = [U]*[S]*[V^H]). However, for the below matrix, it does not produces the same matrix after decomposition. Can you please take a close look at this and let me know what is the problem? Any input on this would be greatly appreciated.

a[0, 0] = ComplexF.Set(4.48113367485427, 0.0F);
a[0, 1] = ComplexF.Set(-4.4811336748542656, 0.0F);
a[0, 2] = ComplexF.Set(9868.5271491812455, 0.0F);
a[0, 3] = ComplexF.Set(9868.5271491812473, 0.0F);
a[0, 4] = ComplexF.Set(0.0F, 0.0F);
a[0, 5] = ComplexF.Set(0.0F, 0.0F);
a[0, 6] = ComplexF.Set(0.0F, 0.0F);
a[0, 7] = ComplexF.Set(0.0F, 0.0F);
a[1, 0] = ComplexF.Set(-0.0070243880579161271, 0.0F);
a[1, 1] = ComplexF.Set(-0.007024388057916128, 0.0F);
a[1, 2] = ComplexF.Set(-0.44008432487412769, 0.0F);
a[1, 3] = ComplexF.Set(0.44008432487412791, 0.0F);
a[1, 4] = ComplexF.Set(0.0F, 0.0F);
a[1, 5] = ComplexF.Set(0.0F, 0.0F);
a[1, 6] = ComplexF.Set(0.0F, 0.0F);
a[1, 7] = ComplexF.Set(0.0F, 0.0F);
a[2, 0] = ComplexF.Set(-3.3706538071975187E+22, 0.0F);
a[2, 1] = ComplexF.Set(5.9574670555112882E-22, 0.0F);
a[2, 2] = ComplexF.Set(-69629.691066082887, 0.0F);
a[2, 3] = ComplexF.Set(-1398.6537438705602, 0.0F);
a[2, 4] = ComplexF.Set(3.1286590272065163E+22, 0.0F);
a[2, 5] = ComplexF.Set(-6.475415127685551E-22, 0.0F);
a[2, 6] = ComplexF.Set(83735.285287348073, 0.0F);
a[2, 7] = ComplexF.Set(198.99031727012286, 0.0F);
a[3, 0] = ComplexF.Set(-5307078354158971, 0.0F);
a[3, 1] = ComplexF.Set(-9.3800034843110651E-29, 0.0F);
a[3, 2] = ComplexF.Set(-0.19776167652692014, 0.0F);
a[3, 3] = ComplexF.Set(0.0039724448727768294, 0.0F);
a[3, 4] = ComplexF.Set(11748278475306892, 0.0F);
a[3, 5] = ComplexF.Set(2.4315522881120115E-28, 0.0F);
a[3, 6] = ComplexF.Set(0.36692418290149531, 0.0F);
a[3, 7] = ComplexF.Set(-0.00087196645140804414, 0.0F);
a[4, 0] = ComplexF.Set(5.283658571380994E+19, 0.0F);
a[4, 1] = ComplexF.Set(9.3386101546107548E-25, 0.0F);
a[4, 2] = ComplexF.Set(3.1051174223656561, 0.0F);
a[4, 3] = ComplexF.Set(-0.0623725890701952, 0.0F);
a[4, 4] = ComplexF.Set(-4.8927292292290486E+19, 0.0F);
a[4, 5] = ComplexF.Set(-1.0126527880191323E-24, 0.0F);
a[4, 6] = ComplexF.Set(-5.7784548297779725, 0.0F);
a[4, 7] = ComplexF.Set(0.013732043259454074, 0.0F);
a[5, 0] = ComplexF.Set(2.1467341077085055E+21, 0.0F);
a[5, 1] = ComplexF.Set(-3.7942483729140873E-23, 0.0F);
a[5, 2] = ComplexF.Set(6.9938324225786337, 0.0F);
a[5, 3] = ComplexF.Set(0.14048532676325157, 0.0F);
a[5, 4] = ComplexF.Set(-1.9866567983179102E+21, 0.0F);
a[5, 5] = ComplexF.Set(4.1118023323984237E-23, 0.0F);
a[5, 6] = ComplexF.Set(-20.106272055444567, 0.0F);
a[5, 7] = ComplexF.Set(-0.047780973596764469, 0.0F);
a[6, 0] = ComplexF.Set(0.0F, 0.0F);
a[6, 1] = ComplexF.Set(0.0F, 0.0F);
a[6, 2] = ComplexF.Set(0.0F, 0.0F);
a[6, 3] = ComplexF.Set(0.0F, 0.0F);
a[6, 4] = ComplexF.Set(-8.1661911579361312E+37, 0.0F);
a[6, 5] = ComplexF.Set(-3.498148996946759E-50, 0.0F);
a[6, 6] = ComplexF.Set(-7.5268734749779691, 0.0F);
a[6, 7] = ComplexF.Set(4.2507101888191323E-05, 0.0F);
a[7, 0] = ComplexF.Set(0.0F, 0.0F);
a[7, 1] = ComplexF.Set(0.0F, 0.0F);
a[7, 2] = ComplexF.Set(0.0F, 0.0F);
a[7, 3] = ComplexF.Set(0.0F, 0.0F);
a[7, 4] = ComplexF.Set(1.3809188481849999E+43, 0.0F);
a[7, 5] = ComplexF.Set(-5.9154381647664128E-45, 0.0F);
a[7, 6] = ComplexF.Set(412.44860073843813, 0.0F);
a[7, 7] = ComplexF.Set(0.0023292532754155336, 0.0F);

Hi Harry (hjsmithh)

Interested to see that you have written a C version of complex SVD and also C software to determine the eigen values/vectors of a complex matrix using the SVD routine. The links you gave previously (to 'geocities') no longer work. Is it possible you have placed this software somewhere else accessible to all? The routines sound exactly what I am looking for. Many thanks, Simon

Hello
I encountered a similar problem in my project.
I need to solve a complex linear system of equations, Ax=b.
A is an hermitian square complex matix of dimesion up to
100.
Could you provide me with complex versions of SVD and the associated files for solving the equations system, similar to the real version in the book ? I need C files.
It would be greatly appreciated.
Regards,
Amir.