====== 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.