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.
matlab diary of SVD demo

<code>
% demo existence of SVD for random matrices
N = 5;
A = rand(N,N)

A =
    0.8147    0.0975    0.1576    0.1419    0.6557
    0.9058    0.2785    0.9706    0.4218    0.0357
    0.1270    0.5469    0.9572    0.9157    0.8491
    0.9134    0.9575    0.4854    0.7922    0.9340
    0.6324    0.9649    0.8003    0.9595    0.6787

[U,D,V] = svd(A);

% check that U,V are unitary and D diagonal
% in matlab ' means hermitian conjugate and * means times

U' * U % calculate hermitian conj U times U

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
   -0.0000    0.0000   -0.0000         0    1.0000

eye(5,5)  % 5 x 5 identity matrix

ans =
     1     0     0     0     0
     0     1     0     0     0
     0     0     1     0     0
     0     0     0     1     0
     0     0     0     0     1

U' * U  - eye(5,5)

ans =
   1.0e-15 *

   -0.2220    0.0555   -0.1561   -0.0555   -0.1110
    0.0555    0.2220    0.1076    0.0139    0.0278
   -0.1561    0.1076    0.6661    0.1076   -0.3018
   -0.0555    0.0139    0.1076   -0.5551         0
   -0.1110    0.0278   -0.3018         0   -0.2220

norm(U' * U  - eye(5,5))

ans =
   7.9587e-16

% so U is unitary

% check V is unitary
norm(V' * V  - eye(5,5))

ans =
   1.2082e-15

% yep!


D

D =
    3.3129         0         0         0         0
         0    0.9431         0         0         0
         0         0    0.8358         0         0
         0         0         0    0.4837         0
         0         0         0         0    0.0198

norm(A - U*D*V')

ans =
   1.7499e-15

% So, we have verified that U D V' is an SVD of A
% and demoed existence of SVD for a given random matrix A


% Next, look at uniqueness. Do this for just 3 x 3 matrix
% to be able to fit several matrices on a screen

N = 3;

% Approach: construct an A via SVD, i.e. start from U,D,V => A
% then calculate svd of A, see if we get the same U,D,V

% trick #1: how to construct a random unitary matrix

[U,R] = qr(randn(N,N)); % do a QR decomposition on a random matrix
[V,R] = qr(randn(N,N)); % do a QR decomposition on a random matrix

norm(U'*U-eye(N))

ans =
   2.5023e-16

norm(V'*V-eye(N))

ans =
   5.3881e-16

% so U and V are unitary

% Construct a diagonal matrix with positive order elements
d = sort(rand(N,1), 'descend')

d =
    0.9001
    0.4893
    0.3377

D = diag(d)

D =
    0.9001         0         0
         0    0.4893         0
         0         0    0.3377

% Set A = U D V'

A = U*D*V'


A =
   -0.0639    0.3319   -0.1336
    0.0952   -0.0322   -0.5831
   -0.7611   -0.1339    0.2899

U

U =
   -0.0868   -0.3192    0.9437
   -0.4657   -0.8244   -0.3217
    0.8807   -0.4674   -0.0771

V

V =

   -0.7878    0.6084   -0.0955
   -0.1464   -0.0343    0.9886
    0.5982    0.7929    0.1161

% Now compute SVD of A using Matlab svd function
% Do we get the same matrices we used to construct A?

[W,E,Z] = svd(A)

% Check that Matlab's svd works
norm(A-W*E*Z')

ans =
   2.2232e-16

% yep!

% Check that original and calculated diagonal matrices are same
diag(D)

ans =
    0.9001
    0.4893
    0.3377

diag(E)

ans =
    0.9001
    0.4893
    0.3377

norm(diag(E)-diag(D))

ans =
   5.6882e-16

% very close indeed!

format long
d = diag(D)

d =

   0.900053846417662
   0.489252638400019
   0.337719409821377

e = diag(E)

e =

   0.900053846417663
   0.489252638400019
   0.337719409821377

% very very close!

format short


% So A == U D V' == W E Z'

% norm of unitary matrix is 1!
norm(U)

ans =
     1

norm(W)

ans =
     1

norm(U-W)

ans =
     2

% Is original U same as W?
U

U =
   -0.0868   -0.3192    0.9437
   -0.4657   -0.8244   -0.3217
    0.8807   -0.4674   -0.0771

W

W =

    0.0868   -0.3192   -0.9437
    0.4657   -0.8244    0.3217
   -0.8807   -0.4674    0.0771

% No! There are sign changes in cols 1 and 3!

% How about original V versus calculated W? Are they the same
norm(V-Z)

ans =
     2

% no!

V

V =
   -0.7878    0.6084   -0.0955
   -0.1464   -0.0343    0.9886
    0.5982    0.7929    0.1161

Z

Z =
    0.7878    0.6084    0.0955
    0.1464   -0.0343   -0.9886
   -0.5982    0.7929   -0.1161

% No! There are sign changes in cols 1 and 3, the same cols as in U,W.

% This illustrates uniqueness of svd of real matrix:
% The singular values are unique, singular vectors are unique up to sign change

% Now do geometrical demo of SVD


d = [4 1 0.25]

d =
    4.0000    1.0000    0.2500

D = diag(d)

D =
    4.0000         0         0
         0    1.0000         0
         0         0    0.2500

A = U*D*V';
A

A =
    0.0569    0.2950   -0.4334
    0.9735    0.2214   -1.7773
   -3.0579   -0.5186    1.7347

U

U =
   -0.0868   -0.3192    0.9437
   -0.4657   -0.8244   -0.3217
    0.8807   -0.4674   -0.0771

V

V =
   -0.7878    0.6084   -0.0955
   -0.1464   -0.0343    0.9886
    0.5982    0.7929    0.1161

% observe action of A on unit ball
% construct matrix X of 400 pts randomly dist on unit sphere
X = randn(3,400);
figure(1); plot3(X(1,:), X(2,:), X(3,:),'k.'); axis equal

%normalize each colum to have unit length
for j=1:400; X(:,j) = X(:,j)/norm(X(:,j)); end
clf
sfigure(1); plot3(X(1,:), X(2,:), X(3,:),'k.'); axis equal
title('X, 400 pts on unit sphere')

% look at Y = AX, action of A on unit ball
Y = A*X;
figure(2); plot3(Y(1,:), Y(2,:), Y(3,:),'k.'); axis equal
title('Y = AX, unit sphere mapped through A')

% Note that norm(A) is equal to first singular value
d

d =
    4.0000    1.0000    0.2500

norm(A)

ans =
    4.0000



% estimate maximum stretching from figure
sqrt(3.5^2 + 2^2)

ans =
    4.0311

sqrt(3.4^2 + 2^2)

ans =
    3.9446

% maximum stretching is indeed 4 
d

d =
    4.0000    1.0000    0.2500


sfigure(2); j=1; plot3([0 U(1,j)], [0 U(2,j)], [0 U(3,j)], 'r-')
plot3(Y(1,:), Y(2,:), Y(3,:),'k.'); axis equal

hold on

% cols of U line up with directions of strechting

% 1st col of U is lined up with dir of max stretching
sfigure(2); j=1; plot3([0 U(1,j)], [0 U(2,j)], [0 U(3,j)], 'r-')

% 2nd col of U is lined up with dir of secondary stretching
sfigure(2); j=2; plot3([0 U(1,j)], [0 U(2,j)], [0 U(3,j)], 'g-')

% 3rd col of U is lined up with dir of minimal stretching
sfigure(2); j=3; plot3([0 U(1,j)], [0 U(2,j)], [0 U(3,j)], 'b-')

sfigure(2); j=1; plot3([0 U(1,j)], [0 U(2,j)], [0 U(3,j)], 'r-', 'Linewidth',2)
sfigure(2); j=3; plot3([0 U(1,j)], [0 U(2,j)], [0 U(3,j)], 'b-', 'Linewidth',2)
sfigure(2); j=2; plot3([0 U(1,j)], [0 U(2,j)], [0 U(3,j)], 'g-', 'Linewidth',2)

% Can do same with cols of V in plot of X unit sphere in figure(1)
</code>
channelflow.org

User Tools

Site Tools

Page Tools
