{VERSION 5 0 "IBM INTEL NT" "5.0" }
{USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 
1 0 0 0 0 1 }{CSTYLE "2D Math" -1 2 "Times" 0 1 0 0 0 0 0 0 2 0 0 0 0 
0 0 1 }{CSTYLE "2D Comment" 2 18 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 }
{CSTYLE "" -1 256 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 
257 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 258 "" 0 1 0 0 
0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 259 "" 0 1 0 0 0 0 1 0 0 0 0 0 
0 0 0 0 }{CSTYLE "" -1 260 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }
{CSTYLE "" -1 261 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 
262 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 263 "" 0 1 0 0 
0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 264 "" 0 1 0 0 0 0 1 0 0 0 0 0 
0 0 0 0 }{CSTYLE "" -1 265 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }
{CSTYLE "" -1 266 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{PSTYLE "Normal
" -1 0 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }1 
1 0 0 0 0 1 0 1 0 2 2 0 1 }{PSTYLE "Heading 1" -1 3 1 {CSTYLE "" -1 
-1 "Times" 1 18 0 0 0 1 2 1 2 2 2 2 1 1 1 1 }1 1 0 0 8 4 1 0 1 0 2 2 
0 1 }{PSTYLE "Title" -1 18 1 {CSTYLE "" -1 -1 "Times" 1 18 0 0 0 1 2 
1 1 2 2 2 1 1 1 1 }3 1 0 0 12 12 1 0 1 0 2 2 19 1 }{PSTYLE "Author" 
-1 19 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }3 
1 0 0 8 8 1 0 1 0 2 2 0 1 }{PSTYLE "Normal" -1 256 1 {CSTYLE "" -1 -1 
"Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }3 1 0 0 0 0 1 0 1 0 2 2 0 1 
}}
{SECT 0 {EXCHG {PARA 18 "" 0 "" {TEXT -1 22 "Graphical Examples of " }
}{PARA 18 "" 0 "" {TEXT -1 32 "the Singular Value Decomposition" }}
{PARA 19 "" 0 "" {TEXT -1 42 "Rob Corless, February 2001\nwww.orcca.on
.ca" }}{PARA 0 "" 0 "" {TEXT -1 296 "We here examine the singular valu
e decomposition of 2 by 2 and 3 by 3 matrices, by explicitly construct
ing the geometric picture of the image of the unit circle (sphere) und
er multiplication by a given matrix. The point of this worksheet is to
 help you to think about matrices in a geometric way." }}{PARA 0 "" 0 
"" {TEXT -1 0 "" }}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 24 "Two-dimension
al examples" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 8 "restart;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 
20 "with(LinearAlgebra):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 148 "We c
hoose a 2 by 2 nonsingular matrix as our first example.  We can return
 to this at our leisure and rework the example with other 2 by 2 examp
les." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 53 "A := Matrix(2,2,[[1
,2],[-1,-1]]); det=Determinant(A);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 
-1 110 "We examine what the matrix A does to unit vectors.  We can par
ameterize the family of unit vectors as follows:" }}}{EXCHG {PARA 0 ">
 " 0 "" {MPLTEXT 1 0 31 "b := Vector(2,[cos(t),sin(t)]);" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 157 "We plot the tips of all possible unit ve
ctors, using the \"scaling=CONSTRAINED\" option so as to make the plot
 axes 1 to 1 (so the circle looks like a circle)." }}}{EXCHG {PARA 0 "
> " 0 "" {MPLTEXT 1 0 90 "circplot := plot([cos(t),sin(t),t=0..2*Pi]):
 plots[display](circplot,scaling=CONSTRAINED);" }}}{EXCHG {PARA 0 "" 
0 "" {TEXT -1 56 "Now consider what the matrix A does to each unit vec
tor:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 10 "xy := A.b;" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 3 "We " }{TEXT 256 5 "claim" }{TEXT 
-1 216 " that this parameterizes an ellipse.  It's pretty easy to veri
fy this algebraically, but you have to remember that it isn't an ellip
se aligned with the coordinate axes; this makes an algebraic proof a l
ittle tedious." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 94 "ellipsepl
ot := plot([xy[1],xy[2],t=0..2*Pi]): plots[display](ellipseplot,scalin
g=CONSTRAINED);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 14 "The so-called \+
" }{TEXT 257 28 "Singular Value Decomposition" }{TEXT -1 28 " characte
rizes that ellipse:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 61 "U,si
gma,Vt := SingularValues(evalf(A),output=['U','S','Vt']);" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 213 "This is a very useful factoring of the m
atrix A.  Both U and V are unitary (so U*U = identity, and V*V = ident
ity; for a rectangular matrix A, the dimensions of U and V are not nec
essarily the same).  We have the " }{TEXT 264 29 "Singular Value Decom
position:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }{TEXT 265 0 "" }}{PARA 
256 "" 0 "" {TEXT -1 6 "A = U " }{XPPEDIT 18 0 "Sigma;" "6#%&SigmaG" }
{TEXT -1 3 " V*" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 31 "A - U.DiagonalMatrix(sigma).Vt;" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "HermitianTranspose(U).U;" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "Vt.HermitianTranspose(Vt);" 
}}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 140 "Now we plot the vectors determ
ined by the columns of U and by the columns of V (rows of Vt).  These \+
directions have meaning for the ellipse." }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 16 "with(plottools):" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 73 "u1v := arrow( [0,0], vector([U[1,1],U[2,1]]), .1, .2,
 .1, colour=green ):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 59 "u2v
 := arrow( [0,0], vector([U[1,2],U[2,2]]), .1, .2, .1 ):" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 58 "plots[display](\{u1v,u2v,ellipseplo
t\},scaling=CONSTRAINED);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 
"" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 56 "The columns of U and of V ar
e called the left and right " }{TEXT 262 16 "singular vectors" }{TEXT 
-1 20 " of A, respectively." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 
0 "" 0 "" {TEXT -1 6 "It is " }{TEXT 258 11 "always true" }{TEXT -1 
233 " that the first column of U will point in the direction of the lo
ngest semiaxis of the ellipse.  Likewise the last column of U will poi
nt in the direction of the shortest semiaxis.  What about the vectors \+
given by the columns of V?  " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 
0 "" 0 "" {TEXT -1 11 "It is also " }{TEXT 259 11 "always true" }
{TEXT -1 76 " that the lengths of the semiaxes are given by the sigmas
 --- the so-called " }{TEXT 260 15 "singular values" }{TEXT -1 17 " of
 the matrix A." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" 
{TEXT -1 121 "Exercise 1: Show that the singular values of A are not t
he eigenvalues of A, unless A is symmetric positive semidefinite." }}
{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 278 "Exercise
 2: Show that the eigenvalues/eigenvectors of A*A and of AA* can give \+
you the singular values and singular vectors of A.  Conclude that exac
t singular values are accessible for matrices of size 4 by 4 and less,
 but not, in general, for matrices of size 5 by 5 or higher." }}{PARA 
0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 135 "Exercise 3: Sh
ow that the eigenvalues of the block matrix [ [ 0 A*], [A 0]] can give
 you the singular values and singular vectors of A." }}{PARA 0 "" 0 "
" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 139 "Exercise 4: what is th
e ratio of the area of the ellipse to the area of the circle?  Does it
 have anything to do with the determinant of A?" }}{PARA 0 "" 0 "" 
{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 183 "Exercise 5: Consider all
 2 by 2 matrices with integer entries, each bounded by n.  Find the no
nsingular matrices with the largest condition number, i.e. ratio of si
gma[1] to sigma[2]." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "
" {TEXT -1 180 "Exercise 6: Derive the formula for the (hyper)volume o
f the n-unit hypersphere, and the formula for the hypervolume of the n
-unit hyperellipsoid.  Generalize the previous exercise." }}{PARA 0 "
" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 43 "Exercise 7: Show t
hat these ideas work for " }{TEXT 266 20 "rectangular matrices" }
{TEXT -1 82 ", so that we can factor A as before (but now U and V will
 be different sizes, and " }{XPPEDIT 18 0 "Sigma;" "6#%&SigmaG" }
{TEXT -1 22 " will be rectangular)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }
}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 75 "v1v := arrow( [0,0], vecto
r([Vt[1,1],Vt[1,2]]), .1, .2, .1, colour=green ):" }}}{EXCHG {PARA 0 "
> " 0 "" {MPLTEXT 1 0 61 "v2v := arrow( [0,0], vector([Vt[2,1],Vt[2,2]
]), .1, .2, .1 ):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 55 "plots[display](\{circplot,v1
v,v2v\},scaling=CONSTRAINED);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 190 
"Note that the green arrow is mapped, by multiplication by A, into the
 direction of the green arrow on the ellipse plot.  Likewise, the tran
sparent arrow is mapped to the corresponding arrow." }}}{EXCHG {PARA 
0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 
0 36 "v1 := HermitianTranspose(Row(Vt,1));" }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 36 "v2 := HermitianTranspose(Row(Vt,2));" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "u1 := Column(U,1);" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "u2 := Column(U,2);" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "A.v1-sigma[1]*u1;" }}}{EXCHG {PARA 
0 "> " 0 "" {MPLTEXT 1 0 17 "A.v2-sigma[2]*u2;" }}}{EXCHG {PARA 0 "" 
0 "" {TEXT -1 93 "That demonstrates that the green vectors correspond,
 as do the transparent arrows, as stated." }}{PARA 0 "" 0 "" {TEXT -1 
0 "" }}{PARA 0 "" 0 "" {TEXT -1 39 "What happens if the matrix is sing
ular?" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 54 "A := Matrix(2,2,[[
1,1],[1,1]]); det := Determinant(A);" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 10 "xy := A.b;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 
50 "plot([xy[1],xy[2],t=0..2*Pi],scaling=CONSTRAINED);" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 59 "We get a degenerate ellipse, one of whose
 semiaxes is zero." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "Singu
larValues(evalf(A));" }}}{PARA 0 "" 0 "" {TEXT -1 78 "Finally, note th
at we can derive a pair of inequalities from this computation:" }}
{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 48 "sigma[n]*
|| x || <= || Ax || <= sigma[1]*|| x ||" }}{PARA 0 "" 0 "" {TEXT -1 0 
"" }}{PARA 0 "" 0 "" {TEXT -1 129 "Note that for certain vectors x, al
igned with v(n) and v(1) respectively, equality is attained (so these \+
inequalities are tight)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "
" 0 "" {TEXT -1 27 "This allows us to define a " }{TEXT 261 11 "matrix
 norm" }{TEXT -1 12 " as follows:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT -1 42 "|| A || := max || Ax ||/|| x || = sigma[1
]" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 68 "and \+
we have the submultiplicative identity (so useful for analysis) " }}
{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 29 "|| Ax || \+
<= || A || || x || ." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "
" {TEXT -1 89 "Exercise: Show that if A is invertible, then sigma[n] >
 0, and || A^(-1) || = 1/sigma[n]." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}
}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 26 "Three-dimensional examples" }}
{PARA 0 "" 0 "" {TEXT -1 48 "Now let us consider a three-dimensional e
xample." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restart;" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "with(LinearAlgebra):" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 51 "sphere := [cos(t)*cos(psi),c
os(t)*sin(psi),sin(t)];" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 55 "
plot3d(sphere,t=0..2*Pi,psi=0..Pi,scaling=CONSTRAINED);" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 45 "A3 := Matrix(3,3,[[1,2,1],[2,1,1],[
-1,1,2]]);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "b3 := Vector(
3,sphere):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "xyz := A3.b3;
" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 85 "As we would expect, we get an
 ellipsoid (again not aligned with the coordinate axes)." }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 71 "plot3d([xyz[1],xyz[2],xyz[3]],t=0..
2*Pi,psi=0..Pi,scaling=CONSTRAINED);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 
-1 95 "Now there are three singular values, being the lengths of the t
hree semiaxes of this ellipsoid." }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 65 "U, sigma, Vt := SingularValues(evalf(A3), output=['U'
,'S','Vt']);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 294 "The ratio of the smallest singular value
 to the largest is a measure of how ``skinny'' the ellipsoid is; we wi
ll see that the skinnier the ellipsoid, the harder the matrix is to de
al with numerically.  This ratio, which an engineer might call the `as
pect ratio' of the ellipsoid, is called the " }{TEXT 263 28 "reciproca
l condition number." }{TEXT -1 136 "  A matrix with a `small' rcond is
 said to be ``ill-conditioned'', which phrase two centuries ago used t
o mean \"rude\" or \"ill-mannered\"." }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 18 "sigma[3]/sigma[1];" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 
-1 39 "Now let us again try a singular system." }}}{EXCHG {PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 13 "A3[1,1] := a;" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 3 "A3;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "Dete
rminant(A3);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "A3[1,1] := \+
7;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 65 "U, sigma, Vt := Singu
larValues(evalf(A3), output=['U','S','Vt']);" }}}{EXCHG {PARA 0 "" 0 "
" {TEXT -1 76 "We see that one of the singular values is, to floating-
point accuracy, zero." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "V \+
:= HermitianTranspose(Vt);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 
18 "v3 := Column(V,3);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 6 "A3
.v3;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 96 "Because every entry is ex
act, here, we can find what this singular vector v3 should be, exactly
:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "LinearSolve(A3,Vector(
3,0));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "subs(_t0[1]=1/sqr
t(1+25+9),%);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 9 "evalf(%);" 
}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{PARA 0 "" 0 "" 
{TEXT -1 0 "" }}}}{MARK "0 0 0" 0 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 }
{PAGENUMBERS 0 1 2 33 1 1 }