3
$\begingroup$

I am working on Singular Value Decomposition for complex matrices. I implemented One Sided Jacobi algorithm. It gives exactly the same result as the svd function in Matlab for the real matrices. However, it does not give the same result in complex matrices. It confirms the formula A = UEV* however I am not sure whether this result is correct or not. Any idea why the result is different? Can I say that the result I obtained is absolutely true if it meets the relation A = USV*?

Do you have any suggestions about algorithms that work stably in complex matrices? Do you know which SVD algorithms using on libraries such as Matlab, Scipy, LAPACK are?

My and Numpy outputs:

A:

[(0.563612363551724+0.4358699131559025j), (0.2806793692533188+0.9930311421704819j), (0.5202215785270062+0.4481106844815471j)] [(0.0011901478455322856+0.6075877153377014j), (0.9997530325963236+0.464449485996173j), (0.8705886268046512+0.6556403407040208j)]
[(0.13985961755470633+0.2579120041586833j), (0.09334365980890702+0.0107473125009665j), (0.9198536550293809+0.6672661113638818j)]

U According to my result:

[(0.32263888225141446+0.49316146595502114j), (0.24956005697294448+0.42711101441813026j), (0.5117539922520571-0.39255845198236494j)]
[(0.5230826883717128+0.4215011626318952j), (0.28920165671330744-0.01699244263936318j), (-0.7012897455688379+0.17167125279675066j)]
[(0.3775020380130938+0.24271259233586367j), (-0.5162076548934359-0.6363169611154663j), (0.24545327028669411+0.04980039176599202j)]

S according to my result :

2.228961805416005 0.7603266403486695 0.6863428620403935

V according to my result:

[(0.4048700775539+0.15820544389995128j), (0.4186462523569478-0.1933392387422724j)(0.850657804557617+0.019852985890792853j)] [(0.5434747415164244+0.03930807376934979j), (0.6282038478106724-0.0013258064898205727j), (-0.559364893463487+0.03670246154675491j)]
[(0.7614928175839254-0.11216861617097869j), (-0.6938072848356582-0.11786215892893791j), (-0.019427233541509478-0.1874673042557197j)]

U according to numpy and Matlab:

[[-0.16336677-0.53547895j 0.66551039-0.11000351j -0.47255786-0.09087257j] [-0.37024417-0.60367276j -0.3400191 -0.14303676j 0.44455477-0.40595198j] [-0.22386647-0.36736338j -0.24426998+0.59097786j -0.10864474+0.62785478j]]

S according to numpy and MATLAB:

[2.26594193 0.91577243 0.47884249]

V according to numpy and MATLAB:

[[-0.36133177+0.j -0.55295939+0.19926312j -0.65938003+0.29864268j] [ 0.39102224+0.j -0.37701504+0.67597439j 0.08382929-0.49090992j] [-0.84648743+0.j 0.06187986+0.2271988j 0.32018676-0.35424718j]]

USV' :(conjugate transpose of V) according to my results ->

[(0.5636123635517245+0.43586991315590246j), (0.2806793692533189+0.9930311421704825j), (0.520221578527007+0.44811068448154784j)]
[(0.0011901478455318415+0.6075877153377023j), (0.9997530325963246+0.4644494859961733j), (0.8705886268046524+0.6556403407040218j)]
[(0.13985961755470683+0.25791200415868387j), (0.09334365980890735+0.010747312500966673j), (0.9198536550293819+0.6672661113638825j)]

A = USV' is confirmed but i could not understand why results are different from ready libraries.

Since my code in python is a bit complicated, I leave the implementation of the same method in matlab here.

function [U,S,V]=jacobi_svd1(A)
% Floating Point function equivalent to MATLAB function svd(A) implemented
% using 1-sided Jacobi algorithm

[m,n]=size(A);                      % Get size of matrix A
U = A;                              % Assign U as A
V=eye(n);                           % Assign V as identity matrix of size n
count=5;                            % Number of sweeps

%%
while(count>=1)

for i = 1:n-1

    for j = i+1:n

         a = norm(U(:,i),2);    % Calculate norm of ith column
         b = norm(U(:,j),2);    % Calculate norm of jit column
         
            % Assure the singular values will appear in decreasing order in S
            % swap columns i and j of U and V
            if a < b

            temp(:,j) = U(:,j);
            U(:,j) = U(:,i);
            U(:,i) = temp(:,j);

            temp1(:,j) = V(:,j);
            V(:,j) = V(:,i);
            V(:,i) = temp1(:,j);

            end
            
            %% 
            % Compute submatrix of U'U
            x=0;
            y=0;
            w=0;

            for k=1:m
            x=x+(U(k,i))^2; 
            end

            for k=1:m
            y=y+(U(k,j))^2; 
            end

            for k=1:m
            w=w+(U(k,i))*(U(k,j)); 
            end
            %% 
            % Compute the Jacobi rotation that diagonalizes the
            % submatrix
            if w ~= 0
                alpha=(y-x)/(2*w);

                if alpha>=0
                    t = 1/(abs(alpha)+sqrt(1+alpha^2));
                else
                    t = -(1/(abs(alpha)+sqrt(1+alpha^2)));
                end

                c=1/sqrt(1+t^2);
                 s=c*t;
            else
                 c=1;
                 s=0;
            end
            %%   
                % update columns i and j of U
                T = U(:,i);
                U(:,i)=c*T-s*U(:,j);
                U(:,j)=s*T+c*U(:,j);
                
                % update matrix V of right singular vectors
                T = V(:,i);
                V(:,i)=c*T-s*V(:,j);
                V(:,j)=s*T+c*V(:,j);


    end 

end
count = count - 1;
end
%%
%singular values are the norms of the columns of U
for j=1:n
singvals(j)=norm(U(:,j),2);
if (singvals(j) > 0)
 U(:,j) = U(:,j)/singvals(j);
end 
end

S=diag(singvals);           % Arrange singular values along the diagonal of S
end

Code reference : https://github.com/jaysmitjadhav/Singular-Value-Decomposition/blob/master/SVD/jacobi_svd1.m

$\endgroup$
1
  • 1
    $\begingroup$ You could check if all the array-creation-by-assignment operations actually produce a copy and not just a reference. In the Jacobi rotations if T is just a reference it could contain/point to, against the obvious intent, the newly computed values after the first part of the rotation, giving a wrong result in the second line. This is just speculation, and I believe it's more likely to happen in python than in matlab. $\endgroup$ Jul 5, 2022 at 5:16

2 Answers 2

5
$\begingroup$

This should not be possible. $U$ and $V$ may be non-unique in the case where there are repeated singular values, but $s$ must be unique, since it is the sorted list of eigenvalues of $A^*A$ and eigenvalues are unique. And max(s) must be equal to the 2-norm of the matrix, for a further check.

So there should be something wrong in what you compute.

One check that you do not mention is that the $U$ and $V$ produced by your algorithm are unitary; are they so? That would be the first additional check I suggest.

$\endgroup$
4
  • 2
    $\begingroup$ you are right i tried now, 2-norm of the matrix is max(s) for (correct singular values). I check also orthagonality, they are not orthagonal :( Somethings are totally wrong for complex case on algorithm. Thank you for interest :) $\endgroup$
    – Emir evcil
    Jul 4, 2022 at 20:59
  • 1
    $\begingroup$ Do you mean the best algorithm? The algorithms for real SVD should work with minor changes; probably you just forgot a bunch of conjugations in your work. For instance in your Matlab implementation I see a loop that produces a sum of squares, and Givens transforms without conjugations. $\endgroup$ Jul 5, 2022 at 6:28
  • $\begingroup$ U and V are generally not unique, even if there are no repeated singular values. See, e.g., math.stackexchange.com/a/644397 $\endgroup$ Jul 11, 2022 at 8:58
  • $\begingroup$ @AmitHochman True, but in that case the non-uniqueness is "tame", as they may only differ by a diagonal scaling. Good point though. $\endgroup$ Jul 11, 2022 at 10:42
1
$\begingroup$

I ran your code in MATLAB to reproduce your $U$,$S$,$V$. You stated "A = USV' is confirmed*", but that is not true. What is true is "A = USV.'", where the .' is the MATLAB non-conjugate transpose!

I cannot quickly pinpoint where your code goes wrong. Personally, I use this:

Algo358 - The singular value decomposition of a complex matrix in C, C++ and C#

$\endgroup$
1
  • $\begingroup$ Thank you so much ! $\endgroup$
    – Emir evcil
    Jul 8, 2022 at 14:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.