Math 291: Intro to Mathematica (for Linear Algebra)

with dusty wilson, mathematics instructor, highline community college
week 10

Contents

- SVD Defuses the Bomb
- SVD in terms of rotations
- SVD and image compression

$Post:=If[MatrixQ[#],MatrixForm[#],#]&

SVD Defuses the Bomb

This is a PowerPoint presentation

SVD in terms of rotations

Mathematica notebook which calculates the singular value decomposition of a matrix as well as generate an animation of the geometric interpretation for two by two case

Note: There is a small error in the code below which causes one of the rotations to go the wrong direction.  I will correct it when I find it.  Other than that, this code is pretty decent.  It is the result of four days of solid work, and so I include it even with its imperfections.

Off[General::spell1];
Off[ParametricPlot::pptr];

A=(({{1, 3}, {-5, 7}}));

m=Dimensions[A][[1]];
n=Dimensions[A][[2]];

V=Eigenvectors[Transpose[A].A];
diag=Eigenvalues[Transpose[A].A]//N;
sortdiag=Sort[diag,Greater]//N;
n=Dimensions[V,2][[2]];
temp={};
For[i=1,i≤n,i++,
AppendTo[temp,Normalize[V[[Position[diag,sortdiag[[i]]][[1]][[1]],All]]]]]
V=Transpose[temp]//N;

U={};
rank=Length[diag]-Count[diag,0.];

For[i=1,i≤rank,i++,
AppendTo[U,(1/(Sqrt[sortdiag[[i]]]))A.V[[All,i]]]];
For[i=rank+1,i≤m,i++,
AppendTo[U,Flatten[ZeroMatrix[1,m]]]];
U=Transpose[U];

If[m==rank && n==rank,
     (*Print[1];*)
    Λ=DiagonalMatrix[(Sqrt[sortdiag[[Table[i,{i,1,rank}]]]])],
    If[m==rank,
      (* Print[2];*)
       Λ=AppendRows[DiagonalMatrix[(Sqrt[sortdiag[[Table[i,{i,1,rank}]]]])],ZeroMatrix[rank,n-rank]],
       If[n==rank,
          (*Print[3];*)
          Λ=AppendColumns[DiagonalMatrix[(Sqrt[sortdiag[[Table[i,{i,1,rank}]]]])],ZeroMatrix[m-rank,rank]],
         (*Print[4];*)
          Λ=BlockMatrix[{{DiagonalMatrix[(Sqrt[sortdiag[[Table[i,{i,1,rank}]]]])],ZeroMatrix[rank,n-rank]}, 
 {ZeroMatrix[m-rank,rank],ZeroMatrix[m-rank,n-rank]}}]]]];

x01[t_]=t;
x02[t_]=t;
x03[t_]=Abs[t];
x04[t_]=0;
y01[t_]=(Sqrt[1-(t^2)]);
y02[t_]=-(Sqrt[1-(t^2)]);
y03[t_]=0;
y04[t_]=Abs[t];

ν=ArcTan[((Transpose[V][[2,1]])/(Transpose[V][[1,1]]))];
x11[t_,θ_]=x01[t] Cos[θ]-y01[t] Sin[θ];
x12[t_,θ_]=x02[t] Cos[θ]-y02[t] Sin[θ];
x13[t_,θ_]=x03[t] Cos[θ]-y03[t] Sin[θ];
x14[t_,θ_]=x04[t] Cos[θ]-y04[t] Sin[θ];
y11[t_,θ_]=x01[t] Sin[θ]+y01[t] Cos[θ];
y12[t_,θ_]=x02[t] Sin[θ]+y02[t] Cos[θ];
y13[t_,θ_]=x03[t] Sin[θ]+y03[t] Cos[θ];
y14[t_,θ_]=x04[t] Sin[θ]+y04[t] Cos[θ];
a=Λ[[1,1]];b=Λ[[2,2]];
x21[t_,α_,β_]=t;
x22[t_,α_,β_]=t;
x23[t_,α_,β_]=Abs[t];
x24[t_,α_,β_]=0;
y21[t_,α_,β_]=(β/α)(Sqrt[(α^2)-(t^2)]);
y22[t_,α_,β_]=-(β/α)(Sqrt[(α^2)-(t^2)]);
y23[t_,α_,β_]=0;
y24[t_,α_,β_]=(β/α) Abs[t];
x31[t_,θ_,α_,β_]=x21[t,α,β] Cos[θ]-y21[t,α,β] Sin[θ];
x32[t_,θ_,α_,β_]=x22[t,α,β] Cos[θ]-y22[t,α,β] Sin[θ];
x33[t_,θ_,α_,β_]=x23[t,α,β] Cos[θ]-y23[t,α,β] Sin[θ];
x34[t_,θ_,α_,β_]=x24[t,α,β] Cos[θ]-y24[t,α,β] Sin[θ];
y31[t_,θ_,α_,β_]=x21[t,α,β] Sin[θ]+y21[t,α,β] Cos[θ];
y32[t_,θ_,α_,β_]=x22[t,α,β] Sin[θ]+y22[t,α,β] Cos[θ];
y33[t_,θ_,α_,β_]=x23[t,α,β] Sin[θ]+y23[t,α,β] Cos[θ];
y34[t_,θ_,α_,β_]=x24[t,α,β] Sin[θ]+y24[t,α,β] Cos[θ];
υ=ArcTan[((U[[2,1]])/(U[[1,1]]))];
x41[t_,θ_,ω_,α_,β_]=x31[t,θ,α,β] Cos[ω]-y31[t,θ,α,β] Sin[ω];
x42[t_,θ_,ω_,α_,β_]=x32[t,θ,α,β] Cos[ω]-y32[t,θ,α,β] Sin[ω];
x43[t_,θ_,ω_,α_,β_]=x33[t,θ,α,β] Cos[ω]-y33[t,θ,α,β] Sin[ω];
x44[t_,θ_,ω_,α_,β_]=x34[t,θ,α,β] Cos[ω]-y34[t,θ,α,β] Sin[ω];
y41[t_,θ_,ω_,α_,β_]=x31[t,θ,α,β] Sin[ω]+y31[t,θ,α,β] Cos[ω];
y42[t_,θ_,ω_,α_,β_]=x32[t,θ,α,β] Sin[ω]+y32[t,θ,α,β] Cos[ω];
y43[t_,θ_,ω_,α_,β_]=x33[t,θ,α,β] Sin[ω]+y33[t,θ,α,β] Cos[ω];
y44[t_,θ_,ω_,α_,β_]=x34[t,θ,α,β] Sin[ω]+y34[t,θ,α,β] Cos[ω];
n=25;

RotateV={};
StretchΛ={};
RotateU={};

unitcircle=ParametricPlot[{{x11[t,0],y11[t,0]},{x12[t,0],y12[t,0]},{x13[t,0],y13[t,0]},{x14[t,0],y14[t,0]}},{t,-1,1},PlotStyle->{Thickness[0.01],Thickness[0.01],RGBColor[1,0,0],RGBColor[0,1,0]},AspectRatio->Automatic,Ticks->False,DisplayFunction->Identity, PlotLabel->"The Unit Circle",PlotRange->9];

For[i=0,i≤n,i++,
ϕ=(ν/n)i;
ψ=(υ/n)i;
fx=1+((a-1)/n)i;
fy=1+((b-1)/n)i;

AppendTo[RotateV,ParametricPlot[{{x11[t,ϕ],y11[t,ϕ]},{x12[t,ϕ],y12[t,ϕ]},{x13[t,ϕ],y13[t,ϕ]},{x14[t,ϕ],y14[t,ϕ]}},{t,-1,1},PlotStyle->{Thickness[0.01],Thickness[0.01],RGBColor[1,0,0],RGBColor[0,1,0]},AspectRatio->Automatic,Ticks->False,DisplayFunction->Identity, PlotLabel->"Rotation : (V^T)",PlotRange->9]];

AppendTo[StretchΛ,ParametricPlot[{{x31[t,ν,fx,fy],y31[t,ν,fx,fy]},{x32[t,ν,fx,fy],y32[t,ν,fx,fy]},{x33[t,ν,fx,fy],y33[t,ν,fx,fy]},{x34[t,ν,fx,fy],y34[t,ν,fx,fy]}},{t,-fx,fx},PlotStyle->{Thickness[0.01],Thickness[0.01],RGBColor[1,0,0],RGBColor[0,1,0]},AspectRatio->Automatic,Ticks->False,DisplayFunction->Identity, PlotLabel->"Componentwise  Stretch : Λ",PlotRange->9]];

AppendTo[RotateU,ParametricPlot[{{x41[t,ν,ψ,a,b],y41[t,ν,ψ,a,b]},{x42[t,ν,ψ,a,b],y42[t,ν,ψ,a,b]},{x43[t,ν,ψ,a,b],y43[t,ν,ψ,a,b]},{x44[t,ν,ψ,a,b],y44[t,ν,ψ,a,b]}},{t,-a,a},PlotStyle->{Thickness[0.01],Thickness[0.01],RGBColor[1,0,0],RGBColor[0,1,0]},AspectRatio->Automatic,Ticks->False,DisplayFunction->Identity, PlotLabel->"Rotation : U",PlotRange->9]]]

SVDMovie=Flatten[{RotateV,StretchΛ,RotateU}];

SVDecomp={{unitcircle,RotateV[[n+1]]},{StretchΛ[[n+1]],RotateU[[n+1]]}};

SVD Output

Print["A = ",A//MatrixForm]
Print["U = ",N[U]//MatrixForm]
Print["Λ = ",N[Λ]//MatrixForm]
Print["V = ",N[V]//MatrixForm]
Print["U·Λ·(V^T) = ",N[U.Λ.Transpose[V]]//MatrixForm]

Show[GraphicsArray[SVDecomp]]

Manipulate[Show[SVDMovie[[i]]],{i,1,Length[SVDMovie],1}]

SVD and image compression

PicPath="http://imagecache.allposters.com/images/pic/EUR/1700-61199~The-Last-Supper-c-1498-Posters.jpg";

(*"http://flightline.highline.edu/dwilson/images/wilson_03-06.jpg",
"http://imagecache.allposters.com/images/pic/HAD/4241~American-Gothic-Posters.jpg",
"http://imagecache.allposters.com/images/pic/PYR/PP30545~Mona-Lisa-Posters.jpg",
*)

PicData=ImageData[Import[PicPath]][[All,All,1]];
PicDim=Dimensions[PicData];
{u,Σ,v}=SingularValueDecomposition[PicData];

Manipulate[Print["Compression to: ",N[100*(n*(PicDim[[1]]+PicDim[[2]]+1))/(PicDim[[1]]*PicDim[[2]])]];Image[Sum[Σ[[i,i]] Take[u,All,{i}].Take[Transpose[v],{i},All],{i,1,n}]],Style["SVD: Image Compression",Bold,12],Delimiter,{n,1,PicDim[[1]],1},
Delimiter,ControlPlacement->Left]