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)