gibson:teaching:fall-2012:iam961:svddemo.m
matlab diary of SVD demo
% 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)
gibson/teaching/fall-2012/iam961/svddemo.m.txt · Last modified: 2012/09/27 14:24 by gibson
Except where otherwise noted, content on this wiki is licensed under the following license: CC Attribution-Share Alike 3.0 Unported
