MATLAB Examples

ROT_SVD2X2

ROT_SVD2X2 computes the complex rotation matrices to obtain the singular value decomposition of a 2-by-2 complex matrix.

Contents

Syntax

    [csl, snl, csr, cor] = rot_svd2x2 (A)
    [csl, snl, csr, cor, D] = rot_svd2x2 (A)
    D = rot_svd2x2 (A)
    [W, D, Z] = rot_svd2x2 (A)
    [U, S, V] = rot_svd2x2 (A, UV_OPT)
    UV_OPT = {'u_comp','v_comp','v_row1_real','s_comp'}

Description

ROT_SVD2X2 can be used when we are solving large matrix SVD problem, which can be divided into many small matrices such as 2x2. As a subroutine, it can be applied to solve the subproblems with speed. It can compute SVD of 2x2 complex matrix at one step finding two rotation matrices, while iteration-based approach may take several steps.

Consider the complex rotation matrices W and Z of a 2x2 matrix A:

A = W*D*Z'
W = [csl snl; -snl' csl']
Z = [csr snr; -snr' csr']

where D is complex diagonal; csl, snl, csr and snr are complex such that

abs(D(1,1)) >= abs(D(2,2))
abs(csl)^2+abs(snl)^2 = 1
abs(csr)^2+abs(snr)^2 = 1

For real A, there exist the rotation angles, theta and phi, such that

csl = cos(theta), snl = sin(theta)
csr = cos(phi), snr = sin(phi)

Once we find W, D and Z, we can rewrite the decomposition of

A = W*sign(D)*S*Z'
U = W*sign(D)
V = Z

where signular values S = abs(D). Generally we can introduce a unit complex diagonal matrix, sign(X):

A = W*sign(D)*sign(X)'*S*(Z*sign(X)')'
U = W*sign(D)*sign(X)'
V = Z*sign(X)'

where sign(X) can be used to meet a phase constraint of U or V. For example, Lapack's SVD returns V whose first row is real, which is ontained by assigning X = diag(Z(1,:)). If we assign X = D, we have

A = W*S*(Z*sign(D)')'
U = W
V = Z*sign(D)'

When called with 4 return values, it returns the rotation matrix entries (csl,snl,csr,snr).

When called with 5 return values, it returns the rotation matrix entries and the complex diagonal D.

When called with one return value, it returns the diagonal D as a column vector.

When called with 3 return value, it returns the rotation matrices W and Z, and the complex diagonal D. If UV_OPT argument is used, it returns U, S, and V in the singular value decomposition.

UV_OPT argument selects a method to obtain singular values S by assigning 'u_comp','v_comp','v_row1_real', or 's_comp':

  • 'u_comp' : the complex sign of D is applied to W, or V = Z;
  • 'v_comp' : the complex sign of D is applied to Z, or U = W;
  • 'v_row1_real' : properly chozen to have the first row of V real;
  • 's_comp' : U, S and V return W, D and Z without changing complex sign.

See also: SVD, Givens rotation, Householer reflection.

Examples

Generate a 2x2 complex matrix:

A = randn(2,2)+i*randn(2,2)
A =

  Column 1

          0.135174942099456 -     0.162337672803828i
          0.515246335524849 -     0.146054634331526i

  Column 2

          0.261406324055383 -     0.532011376808821i
         -0.941485770955434 +      1.68210359466318i

Change format:

format short g

Compute rotation matrices using rot_svd2x2:

[csl, snl, csr, cor] = rot_svd2x2 (A)
csl =

      0.26012


snl =

      0.96557 +  0.0025183i


csr =

      0.22844


cor =

       0.6392 -    0.73432i

Compute the complex diagonal D at the fifth return value:

[csl, snl, csr, cor, D] = rot_svd2x2 (A)
csl =

      0.26012


snl =

      0.96557 +  0.0025183i


csr =

      0.22844


cor =

       0.6392 -    0.73432i


D =

      -2.0255 +    0.42681i            0              
            0                   0.026995 -    0.33621i

Compute D only with one return value. Note that abs(D) is the vector of singular values of A.

D = rot_svd2x2 (A)
abs(D)
svd(A)
norm(abs(D)-svd(A))
D =

      -2.0255 +    0.42681i
     0.026995 -    0.33621i


ans =

         2.07
      0.33729


ans =

         2.07
      0.33729


ans =

   5.5511e-16

Compute W,D,Z with 3 return values:

[W, D, Z] = rot_svd2x2 (A)
err = norm(A-W*D*Z')/norm(A)
W =

      0.26012                    0.96557 +  0.0025183i
     -0.96557 +  0.0025183i      0.26012              


D =

      -2.0255 +    0.42681i            0              
            0                   0.026995 -    0.33621i


Z =

      0.22844                     0.6392 -    0.73432i
      -0.6392 -    0.73432i      0.22844              


err =

   2.5068e-15

Apply the second argument 's_comp', and it returns the complex diagonal:

[W, D, Z] = rot_svd2x2 (A, 's_comp')
err = norm(A-W*D*Z')/norm(A)
W =

      0.26012                    0.96557 +  0.0025183i
     -0.96557 +  0.0025183i      0.26012              


D =

      -2.0255 +    0.42681i            0              
            0                   0.026995 -    0.33621i


Z =

      0.22844                     0.6392 -    0.73432i
      -0.6392 -    0.73432i      0.22844              


err =

   2.5068e-15

Apply the second argument 'u_comp', and it returns SVD with U=W*sign(D):

[U, S, V] = rot_svd2x2 (A, 'u_comp')
err = norm(A-U*S*V')/norm(A)
U =

     -0.25453 +   0.053634i      0.07979 -    0.96227i
      0.94431 -    0.20155i     0.020819 -    0.25929i


S =

         2.07            0
            0      0.33729


V =

      0.22844                     0.6392 -    0.73432i
      -0.6392 -    0.73432i      0.22844              


err =

   2.5032e-15

Apply the second argument 'v_comp', and it returns SVD with V=Z*sign(D'):

[U, S, V] = rot_svd2x2 (A, 'v_comp')
err = norm(A-U*S*V')/norm(A)
U =

      0.26012                    0.96557 +  0.0025183i
     -0.96557 +  0.0025183i      0.26012              


S =

         2.07            0
            0      0.33729


V =

     -0.22353 -   0.047102i      0.78313 +    0.57838i
      0.47406 +    0.85034i     0.018283 +    0.22771i


err =

   2.4767e-15

Apply the second argument 'v_row1_real', and it returns SVD with the real values of the first row V:

[U, S, V] = rot_svd2x2 (A, 'v_row1_real')
err = norm(A-U*S*V')/norm(A)
U =

     -0.25453 +   0.053634i       0.7782 -    0.57161i
      0.94431 -    0.20155i      0.20924 -    0.15454i


S =

         2.07            0
            0      0.33729


V =

      0.22844                    0.97356              
      -0.6392 -    0.73432i      0.14999 +    0.17231i


err =

   2.5032e-15

Derivation based on products of A*A' and A'*A

Derive the rotation matrices for a 2x2 complex matrix. Suppose that A is decomposed such that

A = W*D*Z'
W = [csl snl; -snl' csl']
Z = [csr snr; -snr' csr']

where D is complex diagonal, and W and Z are complex rotation matrices. Each rotation matrix can be computed by digonalizing A*A' or A'*A:

A*A' = [csl snl; -snl' csl']*D*D'*[csl snl; -snl' csl']'
A'*A = [csr snr; -snr' csr']*D'*D*[csr snr; -snr' csr']'

First, we diagonalize A*A' to compute W:

A*A' = [csl snl; -snl' csl']*abs(D)^2*[csl snl; -snl' csl']'
     = [csl snl; -snl' csl']*[dd1 0; 0 dd2]*[csl snl; -snl' csl']'
dd1  = abs(D(1,1))^2
dd2  = abs(D(2,2))^2
A*A' = [ csl  snl ]*[dd1  0]*[ csl  snl ]'
       [-snl' csl'] [0  dd2] [-snl' csl']
     = [ csl*dd1  snl*dd2 ]*[ csl' -snl]
       [-snl'*dd1 csl'*dd2] [ snl'  csl]
(A*A')(:,1)  = [ csl*dd1  snl*dd2 ]*[ csl']
               [-snl'*dd1 csl'*dd2] [ snl']
             = [  csl*dd1*csl'+snl*dd2*snl' ]
               [-snl'*dd1*csl'+csl'*dd2*snl']
(A*A')(:,2)  = [ csl*dd1  snl*dd2 ]*[-snl]
               [-snl'*dd1 csl'*dd2] [ csl]
             = [ -csl*dd1*snl+snl*dd2*csl ]
               [ snl'*dd1*snl+csl'*dd2*csl]
A*A' = [  csl*dd1*csl'+snl*dd2*snl'  , -csl*dd1*snl+snl*dd2*csl ]
       [-snl'*dd1*csl'+csl'*dd2*snl' , snl'*dd1*snl+csl'*dd2*csl]
     = [  dd1*abs(csl)^2+dd2*abs(snl)^2  , -(dd1-dd2)*csl*snl           ]
       [ -(dd1-dd2)*csl'*snl'            , dd1*abs(snl)^2+dd2*abs(csl)^2]
(A*A')(1,1) = dd1*abs(csl)^2+dd2*abs(snl)^2
(A*A')(2,2) = dd1*abs(snl)^2+dd2*abs(csl)^2
(A*A')(1,2) = -(dd1-dd2)*csl*snl

Recall that

abs(csl)^2 + abs(snl)^2)^2 = 1.

The sum and difference of the diagonal of A*A' are:

(A*A')(1,1) + (A*A')(2,2) = dd1 + dd2
(A*A')(1,1) - (A*A')(2,2) = dd1*(abs(csl)^2-abs(snl)^2)
                           +dd2*(abs(snl)^2-abs(csl)^2)
                          = (dd1-dd2)*(abs(csl)^2-abs(snl)^2)

Note that

(abs(csl)^2-abs(snl)^2)^2 + 4*abs(csl*snl)^2 = 1.

From

((A*A')(1,1) - (A*A')(2,2))^2 = (dd1-dd2)^2 * (abs(csl)^2-abs(snl)^2)^2
         abs( (A*A')(1,2) )^2 = (dd1-dd2)^2 * abs(csl*snl)^2,

We can rewrite dd1-dd2 in terms of entries of A*A' without trigometric terms:

(dd1-dd2)^2 = ((A*A')(1,1) - (A*A')(2,2))^2 + 4 * abs( (A*A')(1,2) )^2

or

dd1 - dd2 = sqrt(...)
          = sqrt( ((A*A')(1,1) - (A*A')(2,2))^2 + 4 * abs( (A*A')(1,2) )^2 )

Recall

dd1 + dd2 = (A*A')(1,1) + (A*A')(2,2).

If we assuem dd1 >= dd2, we have

dd1 = ( (A*A')(1,1) + (A*A')(2,2) + sqrt(...) )/2
dd2 = ( (A*A')(1,1) + (A*A')(2,2) - sqrt(...) )/2

Then we can compute the diagonals of S, singular values:

AAt = A*A';
diff_sqr = (AAt(1,1) - AAt(2,2))^2 + 4 * abs( AAt(1,2) )^2; % (dd1-dd2)^2
dd1 = ( AAt(1,1) + AAt(2,2) + sqrt(diff_sqr) )/2;
dd2 = ( AAt(1,1) + AAt(2,2) - sqrt(diff_sqr) )/2;
S = diag([sqrt(dd1) sqrt(dd2)])
S =

         2.07            0
            0      0.33729

Compute the trigometric terms, csl and snl, from dd1 and dd2:

(A*A')(1,1) = dd1*abs(csl)^2+dd2*abs(snl)^2
            = dd1*abs(csl)^2+dd2*(1-abs(csl)^2)
            = (dd1-dd2)*abs(csl)^2 + dd2
(A*A')(2,2) = dd1*abs(snl)^2+dd2*abs(csl)^2
            = dd1*(1-abs(csl)^2)+dd2*abs(csl)^2
            = dd1 - (dd1-dd2)*abs(csl)^2
(A*A')(1,2) = -(dd1-dd2)*csl*snl
abs(csl)^2  =  ( (A*A')(1,1) - dd2 ) / (dd1-dd2)
            =  ( (A*A')(1,1) - (A*A')(2,2) + sqrt(...) ) /2/(dd1-dd2)
            =  ( (A*A')(1,1) - (A*A')(2,2) + (dd1-dd2) ) /2/(dd1-dd2)
            =  ( (A*A')(1,1) - (A*A')(2,2) )/(dd1-dd2)/2 + 1/2
abs(csl)    =  sqrt( ( (A*A')(1,1) - (A*A')(2,2) )/(dd1-dd2)/2 + 1/2 )
snl         =  -(A*A')(1,2)/(dd1-dd2)/csl

We may choose any phase of csl. For simplicity, we use positive real csl:

csl         =  sqrt( ( (A*A')(1,1) - (A*A')(2,2) )/(dd1-dd2)/2 + 1/2 )

Also, note that if dd1==dd2, then csl=1 and snl=0.

We can compute csr and snr similarly using the product A'*A:

abs(csr)    =  sqrt( ( (A'*A)(1,1) - (A'*A)(2,2) )/(dd1-dd2)/2 + 1/2 )
snr         =  -(A'*A)(1,2)/(dd1-dd2)/csr

Note that if (A*A')(1,1) - (A*A')(2,2) == 0, csl is fixed to sqrt(1/2) unstably. We need to find csl without using A*A', but using A! Recall (A*A')(1,1) = (dd1-dd2)*abs(csl)^2 + dd2

abs(csl)^2  =  ( (A*A')(1,1) - dd2 ) / (dd1-dd2)

Compute rotations W and Z:

AtA = A'*A; % for Z
if dd1==dd2,
	csl = 1;
	snl = 0;
	csr = 1;
	snr = 0;
else
	abs_csl = sqrt( ( AAt(1,1) - AAt(2,2) )/(dd1-dd2)/2 + 1/2 );
	csl = abs_csl
	snl =  -AAt(1,2)/(dd1-dd2)/csl
	abs_csr = sqrt( ( AtA(1,1) - AtA(2,2) )/(dd1-dd2)/2 + 1/2 );
	csr = abs_csr
	snr =  -AtA(1,2)/(dd1-dd2)/csr
end
%
W = [csl snl; -snl' csl']
Z = [csr snr; -snr' csr']
csl =

      0.26012


snl =

      0.96557 +  0.0025183i


csr =

      0.22844


snr =

       0.6392 -    0.73432i


W =

      0.26012                    0.96557 +  0.0025183i
     -0.96557 +  0.0025183i      0.26012              


Z =

      0.22844                     0.6392 -    0.73432i
      -0.6392 -    0.73432i      0.22844              

If A is real, we can calculate the rotation angles directly:

See

Alternatively, we may use Givens rotation G to compute csr and snr:

(W'*A) * G = low-triangular

It turns out that Z = G because W is chosen such that W'*A*Z = diagonal. Recall Givens rotation (Lapack style, pseudo code)

function [cs,sn,rr] = rot_givens (f,g);
r = hypot(f,g);
cs = abs(f)/r;
sn = sign(f)*g'/r;
rr = r*sign(f);

Compute csr and snr using Givens rotation G, and compare G and Z:

WtA = W'*A;
[csrr,snrr] = rot_givens (WtA(1,1),WtA(1,2));
G = [csrr -snrr'; snrr csrr']
err_G = norm(G-Z)
G =

      0.22844                     0.6392 -    0.73432i
      -0.6392 -    0.73432i      0.22844              


err_G =

   7.7716e-16

Now we can compute D from the rotation matrices W and Z; We can determine the sign of D:

D11 = [csl; -snl']'*A*[csr; -snr'];
D22 = [snl;  csl']'*A*[snr;  csr'];
D = diag([D11 D22]) % assign D
sign_D = sign(D)
D =

      -2.0255 +    0.42681i            0              
            0                   0.026995 -    0.33621i


sign_D =

     -0.97851 +    0.20619i            0              
            0                   0.080035 -    0.99679i

Note that the singular values are

S = diag(sqrt([dd1 dd2])) = abs(D).

Recall that we need A'*A or W'*A for computing Z. The product W'*A is also used for computing D. Thus, we may prefer to use W'*A for both computing Z and D simultaneously.

Find the right-rotation matrix Z and the diagonal D from the left-rotation matrix W:

WtA = W'*A;
% Givens rotation
f = WtA(1,1);
g = WtA(1,2);
r = hypot(f,g);
if r==0
    csrr = 1;
    snrr = 0;
    rr = 0;
else
    csrr = abs(f)/r;
    snrr = sign(f)*g'/r;
    rr = r*sign(f);
end
% complete Z
Z = [csrr -snrr'; snrr csrr']
% finish D
D11 = rr;
D22 = WtA(2,:)*[-snrr'; csrr'];
D = diag([D11 D22]) % assign D
% compute the complex sign of D:
sign_D = sign(D)
Z =

      0.22844                     0.6392 -    0.73432i
      -0.6392 -    0.73432i      0.22844              


D =

      -2.0255 +    0.42681i            0              
            0                   0.026995 -    0.33621i


sign_D =

     -0.97851 +    0.20619i            0              
            0                   0.080035 -    0.99679i

Finally, we have W, D, and Z as a decomposition of A = W*D*Z'. Speaking of the existence of this decomposition, it is straightforward from the derivation of the rotation matrices.

Compute the decomposition error:

err = norm(A-W*D*Z')/norm(A)
err =

   4.5201e-16

From A = W*D*Z', we can compute the 2x2 SVD: We can rewrite A for any complex diagonal X:

A = W*sign(D)*sign(X)'*S*(Z*sign(X)')'
S = abs(D)
U = W*sign(D)*sign(X)'
V = Z*sign(X)'

There are many U and V candidates for SVD from WDZ decomposition. Among them, the three cases are easily obtained as below:

  • case1) X=eye(2) or U_comp:
U = W*sign_D, V = Z
S = abs(D)
err_svd1 = norm(A-U*S*V')/norm(A)
U =

     -0.25453 +   0.053634i      0.07979 -    0.96227i
      0.94431 -    0.20155i     0.020819 -    0.25929i


V =

      0.22844                     0.6392 -    0.73432i
      -0.6392 -    0.73432i      0.22844              


S =

         2.07            0
            0      0.33729


err_svd1 =

   4.5201e-16

  • case2) X=eye(D^-1) or V_comp:
U = W, V = Z*sign_D'
err_svd2 = norm(A-U*S*V')/norm(A)
U =

      0.26012                    0.96557 +  0.0025183i
     -0.96557 +  0.0025183i      0.26012              


V =

     -0.22353 -   0.047102i      0.78313 +    0.57838i
      0.47406 +    0.85034i     0.018283 +    0.22771i


err_svd2 =

   4.5943e-16

  • case3) X=diag(V(1,:)) or V_row1_real:
X=diag(V(1,:))
U = W*sign_D*sign(X)', V = Z*sign(X)'
err_svd3 = norm(A-U*S*V')/norm(A)
X =

     -0.22353 -   0.047102i            0              
            0                    0.78313 +    0.57838i


U =

        0.238 -    0.10496i      -0.5075 -    0.82145i
     -0.88246 +    0.39193i     -0.13729 -    0.22094i


V =

     -0.22353 +   0.047102i     0.077919 -    0.97043i
      0.77688 +    0.58675i      0.18376 -    0.13572i


err_svd3 =

   4.8479e-16

Compare the errors with SVD's:

[UU,SS,VV] = svd(A)
err_svd = norm(A-UU*SS*VV')/norm(A)
UU =

      0.25453 -   0.053634i       0.7782 -    0.57161i
     -0.94431 +    0.20155i      0.20924 -    0.15454i


SS =

         2.07            0
            0      0.33729


VV =

     -0.22844                    0.97356              
       0.6392 +    0.73432i      0.14999 +    0.17231i


err_svd =

   2.6618e-16

Derivation based on column-uncorelated matrices

Note that algorithm using a product of A and A' with a small off-diagonal quantity produces round-off error when the product is computed. For robust algorithm for such extreme case, we need a different approach. See http://math.stackexchange.com/questions/861674/decompose-a-2d-arbitrary-transform-into-only-scaling-and-rotation

ROT_SVD2X2 is implemented as follows in order to get robust algorithm. ROT_SVD2X2 allows the complex 2x2 matrix, whereas the original robust algorithm can handle real matrix only.

Suppose A is a sum of two column uncorelated matrices: rotation and reflection.

A = [A(1,1) A(1,2)] = [ EE -HH' ] + [FF  GG']
    [A(2,1) A(2,2)]   [ HH  EE' ]   [GG -FF']
EE+FF = A(1,1)
EE-FF = A(2,2)'
HH+GG = A(2,1)
-HH+GG = A(1,2)'

Generate A:

A = randn(2)+i*randn(2)
A =

     -0.87573 -    0.19224i       -0.712 +     1.5301i
     -0.48382 -    0.27407i      -1.1742 -    0.24902i

Compute the components:

EE = (A(1,1)+A(2,2)')/2
FF = (A(1,1)-A(2,2)')/2
GG = (A(2,1)+A(1,2)')/2
HH = (A(2,1)-A(1,2)')/2
QQ = sqrt(abs(EE)^2+abs(HH)^2)
RR = sqrt(abs(FF)^2+abs(GG)^2)
EE =

       -1.025 +   0.028393i


FF =

      0.14924 -    0.22063i


GG =

     -0.59791 -    0.90207i


HH =

      0.11409 +      0.628i


QQ =

       1.2078


RR =

       1.1145

With these terms, re-write A*A' and A'*A:

A*A'  = ( [ EE -HH' ] + [ FF  GG' ] ) * ( [ EE' HH'] + [FF'  GG'] )
          [ HH  EE' ]   [ GG -FF' ]       [-HH  EE]    [GG  -FF ]
      =  (QQ^2+RR^2) * [1 0]
                       [0 1]
      + ( [ EE -HH' ] * [FF'  GG'] ) + ( [ FF  GG'] * [ EE'  HH'] )
          [ HH  EE' ]   [GG  -FF ]       [ GG -FF']   [-HH   EE ]
      = (QQ^2+RR^2) * [1 0]
                      [0 1]
      + [ EE*FF'-HH'*GG EE*GG'+HH'*FF] + [ FF*EE'-GG'*HH FF*HH'+GG'*EE]
        [ HH*FF'+EE'*GG HH*GG'-EE'*FF]   [ GG*EE'+FF'*HH GG*HH'-FF'*EE]
A'*A  = ( [ EE'  HH'] + [FF'  GG'] ) * ( [ EE -HH'] + [ FF  GG'] )
          [-HH   EE]    [GG  -FF ]       [ HH  EE']   [ GG -FF']
      = (QQ^2+RR^2) * [1 0]
                      [0 1]
      + ( [ EE'  HH' ] * [ FF  GG'] ) + ( [ FF'  GG'] * [ EE -HH'] )
          [-HH   EE  ]   [ GG -FF']       [ GG  -FF ]   [ HH  EE']
      = (QQ^2+RR^2) * [1 0]
                      [0 1]
      + [ EE'*FF+HH'*GG  EE'*GG'-HH'*FF'] + [ FF'*EE+GG'*HH -FF'*HH'+GG'*EE']
        [ -HH*FF+EE*GG   -HH*GG'-EE*FF' ]   [  GG*EE-FF*HH   -GG*HH'-FF*EE' ]

Note that Z + Z' = 2 * real(Z).

Then we have:

(A*A')(1,1) = QQ^2+RR^2 + 2 * real(EE*FF'-HH'*GG)
(A*A')(2,2) = QQ^2+RR^2 + 2 * real(HH*GG'-EE'*FF)
(A*A')(1,2) = 2 * (EE*GG'+HH'*FF)

More,

(A*A')(1,1) + (A*A')(2,2) = 2*(QQ^2 + RR^2)
(A*A')(1,1) - (A*A')(2,2) = 4*real(EE*FF'-HH*GG')
dd1 = (QQ + RR)^2
dd2 = (QQ - RR)^2
diff_sqr = ((A*A')(1,1) - (A*A')(2,2))^2 + 4 * abs( (A*A')(1,2) )^2
         = 16 * (real(EE*FF'-HH*GG'))^2 + 16 * abs( (EE*GG'+HH'*FF) )^2
         = 16 * abs( (EE*FF'-HH*GG') )^2 + 16 * abs( (EE*GG'+HH'*FF) )^2
          -16 * (imag(EE*FF'-HH*GG'))^2
         = 16 * ( abs(EE*FF')^2 + abs(HH*GG')^2 + abs(EE*GG')^2 + abs(HH'*FF)^2 )
          - EE*FF'*HH'*GG - EE'*FF*HH*GG' + EE*GG'*HH*FF' + EE'*GG*HH'*FF
          -16 * (imag(EE*FF'-HH*GG'))^2
         = 16 * ( abs(EE*FF')^2 + abs(HH*GG')^2 + abs(EE*GG')^2 + abs(HH'*FF)^2 )
          +16 * (- EE*FF'*HH'*GG - EE'*FF*HH*GG' + EE*FF'*HH*GG' + EE'*FF*HH'*GG)
          -16 * (imag(EE*FF'-HH*GG'))^2
         = 16 * ( abs(EE)^2 + abs(HH)^2 ) * ( abs(FF)^2 + abs(GG)^2 )
          -16 * 4*imag(EE*FF')*imag(HH*GG') - 16 * (imag(EE*FF'-HH*GG'))^2
         = 16 * QQ^2 * RR^2 - 16 * (imag(EE*FF'+HH*GG'))^2
dd1-dd2 = 4* sqrt( (QQ*RR)^2 - (imag(EE*FF'+HH*GG'))^2 )
dd1 = (QQ^2 + RR^2) + sqrt(diff_sqr)/2;
dd2 = (QQ^2 + RR^2) - sqrt(diff_sqr)/2;

Recall

abs(csl)    =  sqrt( ( (A*A')(1,1) - (A*A')(2,2) )/(dd1-dd2)/2 + 1/2 )
snl         =  -(A*A')(1,2)/(dd1-dd2)/csl

Then we have

abs(csl)    =  sqrt( 4*real(EE*FF'-HH*GG')/(dd1-dd2)/2 + 1/2 )
snl         =  -4*(EE*GG'+HH'*FF)/(dd1-dd2)/2/csl

Simliarly,

abs(csr)    =  sqrt( 4*real(EE*FF'+HH*GG')/(dd1-dd2)/2 + 1/2 )
snr         =  -4*(EE'*GG'-HH'*FF')/(dd1-dd2)/2/csr

Note that ROT_SVD2X2 never calls trigonometric functions such as atan2, cos, sin, and so on. It uses sqrt and 4 basic operations.

Compute W,D,Z:

diff_DD = 4 * sqrt( (QQ*RR)^2 - (imag(EE*FF'+HH*GG'))^2 );
dd1 = (QQ^2 + RR^2) + diff_DD/2;
dd2 = (QQ^2 + RR^2) - diff_DD/2;
abs_D = diag([sqrt(dd1) sqrt(dd2)]);
%abs_D = diag(abs([QQ+RR, QQ-RR])) % for real case only
if diff_DD==0
	csl = 1
	snl = 0
	csr = 1
	snr = 0
else
	csl = sqrt( 4*real(EE*FF'-HH*GG')/diff_DD/2 + 1/2 )
	csr = sqrt( 4*real(EE*FF'+HH*GG')/diff_DD/2 + 1/2 )
	if csl == 0
		snl = 1
	else
		snl = -4*(EE*GG'+HH'*FF)/diff_DD/2/csl
	end
	if csr == 0
		snr = 1
	else
		snr = -4*(EE'*GG'-HH'*FF')/diff_DD/2/csr
	end
end
W = [csl snl; -snl' csl]
Z = [csr snr; -snr' csr]
sign_D = diag(sign(diag(W'*A*Z)))
D = abs_D*sign_D
csl =

      0.83059


csr =

      0.42771


snl =

     -0.22391 +    0.50988i


snr =

     -0.45086 +    0.78345i


W =

      0.83059                   -0.22391 +    0.50988i
      0.22391 +    0.50988i      0.83059              


Z =

      0.42771                   -0.45086 +    0.78345i
      0.45086 +    0.78345i      0.42771              


sign_D =

     -0.99965 +   0.026283i            0              
            0                    -0.1882 -    0.98213i


D =

      -2.2807 +   0.059964i            0              
            0                  -0.083511 -    0.43582i

Compute the decomposition error and see how it is better:

err = norm(A-W*D*Z')/norm(A)
err =

   6.5573e-16

Accuracy with a small off-diagonal

Change format:

format long g

Consider an example with small off-diagonal, which is hard to converge with QR iterations.

A = [ sqrt(2) sqrt(eps(1))/2; 0 sqrt(2) ]
A =

           1.4142135623731      7.45058059692383e-09
                         0           1.4142135623731

Compute the diagonals of S, singular values:

AAt = A*A'
diff_sqr = (AAt(1,1) - AAt(2,2))^2 + 4 * abs( AAt(1,2) )^2 % (dd1-dd2)^2
dd1 = ( AAt(1,1) + AAt(2,2) + sqrt(diff_sqr) )/2;
dd2 = ( AAt(1,1) + AAt(2,2) - sqrt(diff_sqr) )/2;
S = diag([sqrt(dd1) sqrt(dd2)])
AAt =

                         2      1.05367121277235e-08
      1.05367121277235e-08                         2


diff_sqr =

      4.44089209850063e-16


S =

          1.41421356609839                         0
                         0           1.4142135586478

Compute rotations:

abs_csl = sqrt( ( AAt(1,1) - AAt(2,2) )/(dd1-dd2)/2 + 1/2 );
csl = abs_csl
snl =  -AAt(1,2)/(dd1-dd2)/csl
AtA = A'*A
abs_csr = sqrt( ( AtA(1,1) - AtA(2,2) )/(dd1-dd2)/2 + 1/2 );
csr = abs_csr
snr =  -AtA(1,2)/(dd1-dd2)/csr
W = [csl snl; -snl' csl']
Z = [csr snr; -snr' csr']
csl =

         0.707106781186548


snl =

        -0.707106785837584


AtA =

                         2      1.05367121277235e-08
      1.05367121277235e-08                         2


csr =

         0.707106781186548


snr =

        -0.707106785837584


W =

         0.707106781186548        -0.707106785837584
         0.707106785837584         0.707106781186548


Z =

         0.707106781186548        -0.707106785837584
         0.707106785837584         0.707106781186548

Compute sign

D11 = [csl; -snl']'*A*[csr; -snr'];
D22 = [snl;  csl']'*A*[snr;  csr'];
sign_D = sign(diag([D11 D22])) % assign D
D = S*sign_D
sign_D =

     1     0
     0     1


D =

          1.41421356609839                         0
                         0           1.4142135586478

Compute error and see how it is poor:

err=norm(A-W*D*Z')/norm(A)
err =

      7.08541996478277e-09

Compare the errors with SVD's:

[UU,SS,VV] = svd(A)
err_svd = norm(A-UU*SS*VV')/norm(A)
UU =

          0.70710678211787        -0.707106780255225
         0.707106780255225          0.70710678211787


SS =

          1.41421356609839                         0
                         0           1.4142135586478


VV =

         0.707106780255225         -0.70710678211787
          0.70710678211787         0.707106780255225


err_svd =

     0

Compute W,D,Z based on robust algorithm:

[W,D,Z] = rot_svd2x2(A)
err = norm(A-W*D*Z')/norm(A)
W =

          0.70710678211787        -0.707106780255225
         0.707106780255225          0.70710678211787


D =

          1.41421356609839                         0
                         0           1.4142135586478


Z =

         0.707106780255225         -0.70710678211787
          0.70710678211787         0.707106780255225


err =

      1.57009245454787e-16

Speed up: SVD2X2 vs. QR-LQ iteration

Compute QR-LQ iteration

A = randn(2)+i*randn(2)
R = A;
for ii=1:100
  [Q,R] = qr(R');
  [Q,R] = qr(R');
  % check R(1,2)
  if abs(R(1,2))<eps(sum(abs(diag(R)))),
      disp([ '>> stopped at step ' num2str(ii)])
      break
  end
end
% Show the converged diagonal
abs(diag(R))
err_iter = norm(abs(diag(R))-abs(rot_svd2x2(A)))
A =

  Column 1

          -1.06421341288933 -      1.50615970397972i
           1.60345729812004 -     0.444627816446985i

  Column 2

           1.23467914689078 -     0.155941035724769i
          -0.22962645096318 +     0.276068253931536i

>> stopped at step 12

ans =

          2.74669405506671
         0.550636361113653


err_iter =

      9.93013661298909e-16

Although iteration-based approach may take several steps or more (one rotation each step), SVD2X2 computes only two rotation steps.

References

See more comments on 2x2 matrix related topics from: