# channelflow.org

### Site Tools

gibson:teaching:fall-2014:iam961:svddemo

# IAM 961 SVD demo

A numerical demo of how uniqueness works for the SVD of real-valued matrices, conducted in Matlab.

```% Form "random" matrix A with known SVD: A = U1 S1 V1'

% Use QR decomposition on random matrix to et random unitary matrix
% The QR decomposition of matrix A is a factorization A = QR with Q
% unitary Q and upper-triangular R. The QR decomp of a matrix with
% normally distributed random elements with produce a unitary matrix
% with its orthogonal columns in an isotropically random orientation.

[U1,tmp] = qr(randn(5,5));
[V1,tmp] = qr(randn(5,5));

>> U1

U1 =
-0.1735   -0.2978    0.4339    0.8304    0.0575
-0.5919   -0.1608   -0.7489    0.2185   -0.1230
0.7291    0.1574   -0.4449    0.4549   -0.1970
-0.2783    0.7273    0.1992    0.1386   -0.5785
-0.1029    0.5760   -0.1147    0.1910    0.7798

% Demonstrate that U1 is unitary, in a variety of ways
>> U1 * U1'

ans =
1.0000   -0.0000    0.0000   -0.0000   -0.0000
-0.0000    1.0000   -0.0000   -0.0000   -0.0000
0.0000   -0.0000    1.0000    0.0000    0.0000
-0.0000   -0.0000    0.0000    1.0000    0.0000
-0.0000   -0.0000    0.0000    0.0000    1.0000

>> U1 * U1' - eye(5)

ans =
1.0e-15 *

0   -0.0278    0.0798   -0.0347   -0.0069
-0.0278         0   -0.0798   -0.2914   -0.1249
0.0798   -0.0798   -0.1110    0.0833    0.0278
-0.0347   -0.2914    0.0833    0.4441    0.1665
-0.0069   -0.1249    0.0278    0.1665    0.2220

>> norm(U1 * U1' - eye(5))

ans =
6.9598e-16

% Show V1 is unitary
>> norm(V1 * V1' - eye(5))

ans =
3.0661e-16

% Now construct S1 with known diagonal entries
>> S1 = diag([127, 34, 3.5, 0.72, 0.00573])

S1 =
127.0000         0         0         0         0
0   34.0000         0         0         0
0         0    3.5000         0         0
0         0         0    0.7200         0
0         0         0         0    0.0057

>> A = U1*S1*V1'

A =
15.3394   15.7662   -2.5713    9.0749   -4.2474
50.8731   39.5276  -13.2333   13.7956  -34.2276
-61.9380  -46.1719   16.7968  -20.9126   43.7778
22.1138    1.4313  -12.9916   -1.3913  -34.6287
7.6613   -5.5310   -7.2904   -4.8075  -19.7186

% Notation:
%   use U1,S1,V1 for the constructed, known SVD
%   use U2,S2,V2 for the calculated SVD

% Now calculate the SVD of A
>> [U2,S2,V2] = svd(A);

% Verify that the calculated SVD actually equals A
>> U2 * S2 * V2'

ans =
15.3394   15.7662   -2.5713    9.0749   -4.2474
50.8731   39.5276  -13.2333   13.7956  -34.2276
-61.9380  -46.1719   16.7968  -20.9126   43.7778
22.1138    1.4313  -12.9916   -1.3913  -34.6287
7.6613   -5.5310   -7.2904   -4.8075  -19.7186

>> A

A =
15.3394   15.7662   -2.5713    9.0749   -4.2474
50.8731   39.5276  -13.2333   13.7956  -34.2276
-61.9380  -46.1719   16.7968  -20.9126   43.7778
22.1138    1.4313  -12.9916   -1.3913  -34.6287
7.6613   -5.5310   -7.2904   -4.8075  -19.7186

>> norm(A-U2*S2*V2')

ans =
7.6825e-14

>> norm(A-U1*S1*V1')

ans =
0

% Now compare the factors of the known and calculated SVDs

>> U1

U1 =
-0.1735   -0.2978    0.4339    0.8304    0.0575
-0.5919   -0.1608   -0.7489    0.2185   -0.1230
0.7291    0.1574   -0.4449    0.4549   -0.1970
-0.2783    0.7273    0.1992    0.1386   -0.5785
-0.1029    0.5760   -0.1147    0.1910    0.7798

>> U2

U2 =
0.1735   -0.2978    0.4339    0.8304    0.0575
0.5919   -0.1608   -0.7489    0.2185   -0.1230
-0.7291    0.1574   -0.4449    0.4549   -0.1970
0.2783    0.7273    0.1992    0.1386   -0.5785
0.1029    0.5760   -0.1147    0.1910    0.7798

>> U1-U2

ans =
1.0e-15 *

-0.1110    0.4441   -0.2776         0         0
0.1110   -0.0833    0.1110   -0.1110    0.0971
-0.1110    0.2498         0         0   -0.0555
-0.0555   -0.2220   -0.1388   -0.3886   -0.1110
-0.0278         0    0.0139   -0.1110         0

% Amazing: U1 and U2 are actually equal, even though the SVD
% is unique only up to +/- 1 in the columns of U and V.

% Will demonstrate this by changing the sign of 1st columns of U2
% and V2 and showing that the modified SVD still satisfies  A = U2 S2 V2'

>> U2(:,1)= -1*U2(:,1);
>> V2(:,1)= -1*V2(:,1);
>> V1

V1 =
-0.6683   -0.0590   -0.1033    0.2860    0.6763
-0.4695   -0.6019   -0.3716   -0.1836   -0.4955
0.1960   -0.2385   -0.1227   -0.8035    0.4939
-0.1898   -0.3528    0.9098   -0.1026   -0.0359
0.5085   -0.6730   -0.0920    0.4778    0.2277

>> V2
V2 =

0.6683   -0.0590   -0.1033    0.2860    0.6763
0.4695   -0.6019   -0.3716   -0.1836   -0.4955
-0.1960   -0.2385   -0.1227   -0.8035    0.4939
0.1898   -0.3528    0.9098   -0.1026   -0.0359
-0.5085   -0.6730   -0.0920    0.4778    0.2277

% Now V1 and V2 differ in sogn of 1st col (same w U1, U2), but still

>> norm(A-U2*S2*V2')

ans =
7.6825e-14

% and note that the diagonal matrices are the same.
>> S1

S1 =
127.0000         0         0         0         0
0   34.0000         0         0         0
0         0    3.5000         0         0
0         0         0    0.7200         0
0         0         0         0    0.0057

>> S2

S2 =
127.0000         0         0         0         0
0   34.0000         0         0         0
0         0    3.5000         0         0
0         0         0    0.7200         0
0         0         0         0    0.0057

% Ok, now let's really mess with the uniqueness of the SVD
% by setting two singular values to the same value.

>> S1 = diag([127, 34, 34, 0.72, 0.00573])

S1 =
127.0000         0         0         0         0
0   34.0000         0         0         0
0         0   34.0000         0         0
0         0         0    0.7200         0
0         0         0         0    0.0057

% Reform A from the new S1 and recompute the calculated SVD
>> A = U1*S1*V1';
>> [U2,S2,V2] = svd(A);

% By the way, what's the 2-norm of A?

>> norm(A)

ans =
127.0000

% The 2-norm of a matrix is its largest singular value!

% Compare U1 and U2
>> U1

U1 =
-0.1735   -0.2978    0.4339    0.8304    0.0575
-0.5919   -0.1608   -0.7489    0.2185   -0.1230
0.7291    0.1574   -0.4449    0.4549   -0.1970
-0.2783    0.7273    0.1992    0.1386   -0.5785
-0.1029    0.5760   -0.1147    0.1910    0.7798

>> U2

U2 =
-0.1735    0.2294   -0.4737   -0.8304    0.0575
-0.5919   -0.7302    0.2314   -0.2185   -0.1230
0.7291   -0.3084    0.3572   -0.4549   -0.1970
-0.2783    0.5334    0.5330   -0.1386   -0.5785
-0.1029    0.1858    0.5571   -0.1910    0.7798

diag(S1)'

ans =

127.0000   34.0000   34.0000    0.7200    0.0057

% Note that since sigma_2 = sigma_3, the 2nd and 3rd cols
% of U (and V) are indeterminate! The prescribed SVD and
% calculated SVD now differ in 2nd and 3rd cols of U and V

% Let's make three equal singular values

>> S1 = diag([127, 34, 34, 34, 0.00573])
S1 =
127.0000         0         0         0         0
0   34.0000         0         0         0
0         0   34.0000         0         0
0         0         0   34.0000         0
0         0         0         0    0.0057

>> A = U1*S1*V1';
>> [U2,S2,V2] = svd(A)
>> V1

V1 =
-0.6683   -0.0590   -0.1033    0.2860    0.6763
-0.4695   -0.6019   -0.3716   -0.1836   -0.4955
0.1960   -0.2385   -0.1227   -0.8035    0.4939
-0.1898   -0.3528    0.9098   -0.1026   -0.0359
0.5085   -0.6730   -0.0920    0.4778    0.2277

>> V2

V2 =
-0.6683         0   -0.3098   -0.0000   -0.6763
-0.4695   -0.7264   -0.0690    0.0405    0.4955
0.1960   -0.5303    0.6556   -0.0814   -0.4939
-0.1898    0.1663    0.3311    0.9085    0.0359
0.5085   -0.4042   -0.5999    0.4079   -0.2277

>> diag(S1)'
ans =
127.0000   34.0000   34.0000   34.0000    0.0057

% Notice now that the 2nd, 3rd, and 4th cols of V differ.

% But the 2nd, 3rd, and 4th cols of V1 span the same 3d subspace
% of R^5 and the same cols of V2. Can show this by expanding
% the 2nd (and 3rd and 4th) col of V1 in the 2,3,4 cols of V1

% Let x be the 2nd col of V1
>> x = V1(:,2);

% and v2,v3,v4 by the 2,3,4 cols of V2
>> v2 = V2(:,2);
>> v3 = V2(:,3);
>> v4 = V2(:,4);

% Show that x is in span(v2,v3,v4)
>> x
x =
-0.0590
-0.6019
-0.2385
-0.3528
-0.6730

% by expanding x in the 3d orthonormal set (v2,v3,v4)
>> (v2'*x)*v2 + (v3'*x)*v3 + (v4'*x)*v4
ans =
-0.0590
-0.6019
-0.2385
-0.3528
-0.6730

>> norm(x - ((v2'*x)*v2 + (v3'*x)*v3 + (v4'*x)*v4))
ans =
2.7336e-16

% Hooray! x (2nd col of V1) lies in the subspace spanned by the 2,3,4 cols of V2.
```