channelflow.org

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
====== IAM 961 SVD demo ======

A numerical demo of how uniqueness works for the SVD of real-valued matrices, conducted in Matlab.
<code>
% 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.

</code>
channelflow.org

User Tools

Site Tools

Page Tools
