MATLAB Examples

ROT_DSVD2X2

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

Contents

Syntax

    [csl, snl, csr, cor] = rot_dsvd2x2 (A)
    [csl, snl, csr, cor, D] = rot_dsvd2x2 (A)
    D = rot_dsvd2x2 (A)
    [W, D, Z] = rot_dsvd2x2 (A)

Description

ROT_DSVD2X2 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 real 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

It happens that D(2,2) <0. It is because pure rotations cannot diagonalize the matrix with non-negative diagonal entries. If we want non-negative diagonals like SVD, we need an additional operation to change sign.

ROT_DSVD2X2 is derived from that A is written by a sum of rotation and reflection parts:

A = [ EE -HH ] + [FF  GG]
    [ HH  EE ]   [GG -FF]

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.

See also: SVD, Givens rotation, Householer reflection.

Examples

Generate a 2x2 matrix:

A = randn(2,2)
A =

         -1.08906429505224         0.552527021112224
        0.0325574641649735          1.10061021788087

Change format:

format short g

Compute rotation matrices using rot_svd2x2:

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

      0.78633


snl =

      -0.6178


csr =

      -0.6002


cor =

     -0.79985

Compute the complex diagonal D at the fifth return value:

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

      0.78633


snl =

      -0.6178


csr =

      -0.6002


cor =

     -0.79985


D =

       1.3933            0
            0      -0.8732

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

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

       1.3933
      -0.8732


ans =

       1.3933
       0.8732


ans =

       1.3933
       0.8732


ans =

   1.1102e-16

Compute W,D,Z with 3 return values:

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

      0.78633      -0.6178
       0.6178      0.78633


D =

       1.3933            0
            0      -0.8732


Z =

      -0.6002     -0.79985
      0.79985      -0.6002


err =

   2.4434e-16

Unlike SVD, D(2,2) may be negative. See the following 2x2 matrix:

A = [1 1; 1 -1]
A =

     1     1
     1    -1

Compute W,D,Z:

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

      0.92388     -0.38268
      0.38268      0.92388


D =

       1.4142            0
            0      -1.4142


Z =

      0.92388     -0.38268
      0.38268      0.92388


err =

   7.8505e-17

It is obviously that we can have positive diagonals by changing sign of columns like SVD:

U = W
S = D; S(2,2) = -S(2,2)
V = Z; V(:,2) = -V(:,2)
err = norm(A-U*S*V.')/norm(A)
U =

      0.92388     -0.38268
      0.38268      0.92388


S =

       1.4142            0
            0       1.4142


V =

      0.92388      0.38268
      0.38268     -0.92388


err =

   7.8505e-17

Diagnostics

Derive the rotation matrices for a 2x2 real 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 diagonal, and W and Z are rotation matrices. ' denotes transpose. 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]'

Then we have

(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

where

dd1  = abs(D(1,1))^2
dd2  = abs(D(2,2))^2

Then dd1 and dd2 are calculated from A*A':

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

where

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

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]

Generate A:

A = randn(2)
A =

       1.5442      -1.4916
     0.085931      -0.7423

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
EE =

      0.40096


FF =

       1.1433


GG =

     -0.70283


HH =

      0.78876

Now, we try to write W, D, and Z in terms of EE, FF, GG, and HH. Note that A*A' or A'*A may suffer from round-off error, which is not stable.

From column uncorelatedness, the product with their own transposed matrices is diagonal:

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

where QQ and RR are scaling factors of each rotation matrix.

We define and compute QQ and RR:

QQ = sqrt(EE^2+HH^2)
RR = sqrt(FF^2+GG^2)
QQ =

      0.88482


RR =

        1.342

We write A*A' and A'*A again:

A*A'  = ( [ EE -HH ] + [FF  GG] ) * ( [ EE  HH ] + [FF  GG] )
          [ HH  EE ]   [GG -FF]       [-HH  EE ]   [GG -FF]
      = [ EE^2+HH^2         0 ]
        [         0 EE^2+HH^2 ]
      + ( [ EE -HH ] * [FF  GG] ) + ( [ FF  GG] * [ EE  HH ]   )
          [ HH  EE ]   [GG -FF]       [ GG -FF]   [-HH  EE ]
      + [ FF^2+GG^2         0 ]
        [         0 FF^2+GG^2 ]
      = (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]
      = [ EE^2+HH^2         0 ]
        [         0 EE^2+HH^2 ]
      + ( [ EE  HH ] * [FF  GG] ) + ( [ FF  GG] * [ EE -HH ]   )
          [-HH  EE ]   [GG -FF]       [ GG -FF]   [ HH  EE ]
      + [ FF^2+GG^2         0 ]
        [         0 FF^2+GG^2 ]
      = (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 ]

We have the components of A*A' in terms of QQ,RR,EE,FF,HH, and GG:

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

Note

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

Remark that (A*A')(1,1) - (A*A')(2,2) is reduced 2 product terms, where the previous 4 product term easily suffers from numerical errors.

AAt = A*A';
AAt(1,1) - AAt(2,2)
4*(EE*FF-HH*GG)
ans =

        4.051


ans =

        4.051

We can rewrite dd1 and dd2:

diff_sqr = 16*( EE*FF-HH*GG )^2 + 16 * abs( EE*GG+HH*FF )^2
         = 16*( (EE*FF)^2 + (HH*GG)^2 + (EE*GG)^2 + (HH*FF)^2 )
         = 16*( (EE^2+HH^2) * (FF^2+GG^2) )
         = 16*( QQ^2 * RR^2 )
dd1-dd2 = 4*sqrt( QQ^2 * RR^2 ) = 4*QQ*RR
dd1 = QQ^2 + RR^2 + 2*QQ*RR = (QQ + RR)^2
dd2 = (QQ - RR)^2

Compute D:

D = diag([QQ+RR, QQ-RR])
D =

       2.2268            0
            0     -0.45719

To compute the rotation matrices, recall:

(A*A')(1,1) + (A*A')(2,2) = dd1 + dd2
(A*A')(1,1) - (A*A')(2,2) = (dd1 - dd2)*cos(2*theta_l)
-2 * (A*A')(1,2) = (dd1 - dd2)*sin(2*theta_l)

where theta_l is such that

cos(theta_l) = csl, sin(theta_l) = snl.

Then we have:

tan(2*theta_l) = -2 * (A*A')(1,2) / ( (A*A')(1,1) - (A*A')(2,2) )
               = -2 * 2*(EE*GG+HH*FF) / ( 4*(EE*FF-HH*GG) )
               = (EE*GG+HH*FF) / (-EE*FF+HH*GG)
               = - (GG/FF+HH/EE) / (1-HH/EE*GG/FF)

Note tan(a+b) = ( tan(a) + tan(b) )/( 1 - tan(a) * tan(b) )

Therefore,

tan(2*theta_l) = - tan(a+b), tan(a)=GG/FF, tan(b)=HH/EE

or

theta_l = - ( arctan(GG/FF) + arctan(HH/EE) )/2

Compute the left rotation:

theta_l = - ( atan2(GG,FF) + atan2(HH,EE) )/2
csl = cos(theta_l)
snl = sin(theta_l)
W = [csl snl; -snl csl]
theta_l =

     -0.27465


csl =

      0.96252


snl =

     -0.27121


W =

      0.96252     -0.27121
      0.27121      0.96252

Similarly,

tan(2*theta_r) = -2 * (A'*A)(1,2) / ( (A'*A)(1,1) - (A'*A)(2,2) )
               = -2 * 2*(EE*GG-HH*FF) / ( 4*(EE*FF+HH*GG) )
               = - (EE*GG-HH*FF) / (EE*FF+HH*GG)
               = - (GG/FF-HH/EE) / (1+HH/EE*GG/FF)

Note tan(a-b) = ( tan(a) - tan(b) )/( 1 + tan(a) * tan(b) )

theta_r = - ( arctan(GG/FF) - arctan(HH/EE) )/2

Compute the right rotation:

theta_r = - ( atan2(GG,FF) - atan2(HH,EE) )/2
csr = cos(theta_r)
snr = sin(theta_r)
Z = [csr snr; -snr csr]
theta_r =

      0.82585


csr =

      0.67793


snr =

      0.73513


Z =

      0.67793      0.73513
     -0.73513      0.67793

Compute the WDZ decomposition error:

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

   2.8237e-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 based on A*A':

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 poor it is:

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

      7.08541996478277e-09

Now we try to compute the rotation-reflection sum of A:

A = [ EE -HH ] + [FF  GG]
    [ HH  EE ]   [GG -FF]
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
EE =

           1.4142135623731


FF =

     0


GG =

      3.72529029846191e-09


HH =

     -3.72529029846191e-09

Compute diagonals from EE,FF,GG, and HH:

QQ = sqrt(EE^2+HH^2)
RR = sqrt(FF^2+GG^2)
DD = diag([QQ+RR, QQ-RR])
QQ =

           1.4142135623731


RR =

      3.72529029846191e-09


DD =

          1.41421356609839                         0
                         0           1.4142135586478

Compute rotations:

alpha = atan2(GG,FF);
beta  = atan2(HH,EE);
phi   = -(alpha+beta)/2
theta = -(alpha-beta)/2
phi =

        -0.785398162080359


theta =

        -0.785398164714537

Note

tan(alpha) = GG/FF, cos(alpha) = FF/RR, sin(alpha) = GG/RR
tan(beta)  = HH/EE, cos(beta)  = EE/QQ, sin(beta)  = HH/QQ

Then we have the rotations:

csl=cos(phi);
snl=sin(phi);
csr=cos(theta);
snr=sin(theta);
WW = [csl snl; -snl csl]
ZZ = [csr snr; -snr csr]
WW =

          0.70710678211787        -0.707106780255225
         0.707106780255225          0.70710678211787


ZZ =

         0.707106780255225         -0.70710678211787
          0.70710678211787         0.707106780255225

Compute and compare error

err=norm(A-WW*DD*ZZ')/norm(A)
err =

      2.01093737165851e-16

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

References

See more comments on 2x2 matrix related topics from: