MATLAB Examples

# ROT_DSVD2X2

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

## 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: