MATLAB Examples

# ROT_SVD2X2

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

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

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