MGM.07 Ill-conditioned Matrices and Roundoff
 BASICS

Needs["LinearAlgebra`MatrixManipulation`"];
Needs["Graphics`Colors`"];
Off[General::"spell"];
Off[General::"spell1"];
Off[ParametricPlot3D::"ppcom"];
Off[ParametricPlot::"ppcom"];
Off[Plot::"plnr"];
SetOptions[Limit, Analytic -> True];
$Pre = Chop;

CMView = {2.7, 1.6, 1.2};

Vector Drawer

\(\(\( (*\ \(\(:\)\(Date : \ \ Copyright\ 2001\)\), \ 
      Math\ Everywhere, \ \(\(Inc\)\(.\)\)\ 
*) \)\(
\)\( (*
Original\ Conception\ to\ create\ a\ vector\ drawer\ is\ due\ to\ Bill\ Davis . 
Current\ code\ is\ by\ Bruce\ \(\(Carpenter\)\(.\)\)
*) \)\(
\)\( (*\ \ \(\(\(:\)\(Drawbacks : \ \ Since\ the\ vectors\ are\ drawn\ out\ of\)\)\ 

      \ \ \ \ \ context\), \ 
      the\ shapes\ of\ the\ arrowheads\ are\ only\ 
\ \ \ \ \ correct\ with\ AspectRatio\ set\ to\ \(\(Automatic\)\(.\)\)
*) \)\)\)

BeginPackage["Vector3D`","Graphics`Colors`"];

Axes3D::usage = "Axes3D[a, b] creates a Graphics3D object of cartesian axes with x, y, and z running from -a/3 to a, and with axes labels b units beyond the tips of the axes.  Axes3D[a] is Axes3D[a, a/8].";

Perpend::usage = "Perpend[a] returns a unit vector which is perpendicular to the vector a.";

Arrow::usage = "Arrow[a] produces a 2 or 3 dimensional vector from the origin to a. Arrow[a, Tail -> tail] gives a vector from tail to a.";
	
Vector::usage = "Vector[a, b] gives a 2D or 3D vector from a to b.";

ArrowHead::usage = "ArrowHead[a, vec] produces an arrowhead with its tip placed at the point a, pointing in the direction of the vector vec.";

Tail::usage = "Tail -> point puts the tail of the vector at point.";

Norm::usage = "Norm[vec] returns the Euclidean length of vec.";

Normalize::usage = "Normalize[vec] returns a unit vector in the direction of vec, provided vec is nonzero.";

Aperture::usage = "Aperture is the ratio of the radius of the base to the length of the head of a vector.";

TipSize::usage = "TipSize specifies an absolute size for the inner tip length of the head of a vector.";

TipRatio::usage = "TipRatio is the ratio of the inner length of the arrowhead to the length of the shaft of a vector.";

HeadRatio::usage = "HeadRatio is the ratio of the outer length of the arrowhead to the length of the shaft of a vector.";

HeadSize::usage = "HeadSize specifies an absolute size for the head of a vector.";

ScaleFactor::usage = "ScaleFactor specifies the amount to scale a vector in length and may be either a number or a pure function.";

ZeroVectorPointSize::usage = "ZeroVectorPointSize is the size of the point used to represent a zero vector.";

Options[Vector] = Options[Arrow] = Options[ArrowHead] =
{HeadRatio -> 0.18, TipRatio -> 0.14, Aperture -> 0.3, ScaleFactor -> 1, HeadType -> Polygon, ShaftQ -> True, EdgesQ -> True, VectorColor -> RGBColor[0,0,1], ShaftWidth -> 0.005, ZeroVectorPointSize -> 0.01};
		
Begin["`Private`"];

SetOptions[ParametricPlot,AspectRatio -> Automatic];
SetOptions[Plot,AspectRatio -> Automatic];
SetOptions[Graphics,AspectRatio -> Automatic];

Axes3D[u_,v_] := 
      Graphics3D[{
            {Blue,Line[{{-u/3,0,0},{u,0,0}}]},
                    Text["x",{u + v,0,0}],
            {Blue,Line[{{0,-u/3,0},{0,u,0}}]},
                    Text["y",{0, u + v,0}],
            {Blue,Line[{{0,0,-u/3},{0,0,u}}]},
                    Text["z",{0,0,u + v}]}]; 
Axes3D[u_] := Axes3D[u,u/8];
				
Norm[a_?VectorQ] := N[Sqrt[a.a]];

Normalize[a_?VectorQ] := If[a.a == 0, a, a/Norm[a]];

Perpend[a_?VectorQ] := Normalize[{a[[2]],-a[[1]]}]/;Length[a]==2;

Perpend[a_?VectorQ]:=
    Normalize[If[a[[1]]==0,{1,0,0},{a[[2]],-a[[1]],0}]]/;
      Length[a]==3;
      
Base3D = Table[{Cos[2 k Pi/8.], Sin[2 k Pi/8.]},{k,0,8}]; 

Vector[from_?VectorQ, to_?VectorQ, opts___] := Module[
   {sq, sw, ht, hs, hr, tr, ap, sf, vc, zvps, diff = to - from, tip},
   {sq, sw, ht, hs, hr, ts, tr, ap, sf, vc, zvps} = {ShaftQ, ShaftWidth, HeadType, HeadSize, HeadRatio, TipSize, TipRatio, Aperture, ScaleFactor, VectorColor, ZeroVectorPointSize} /. {opts} /. Options[Vector];
   If[from==to, Return[Graphics[{PointSize[zvps],vc,Point[from]}]]];
   diff = If[NumberQ[N[sf]], sf diff, sf[diff]];
   tip = from + diff;
   If[NumberQ[hs], hr = hs/Norm[diff]; tr = 0.75 hr];
   If[NumberQ[ts], tr = ts/Norm[diff]];
   Graphics[{If[sq,{vc,Thickness[sw],Line[{from,tip - tr diff}]},{}],
      With[{trans = ap hr Norm[diff] Perpend[diff]},
      {vc, Thickness[sw], ht[{tip,
                   (tip - hr diff + trans),
                   (tip - tr diff),
                   (tip - hr diff - trans),tip}]}]}]
] /; Length[from]==Length[to]==2

Vector[from_?VectorQ, to_?VectorQ, opts___] := Module[
   {sq, eq, sw, hs, hr, tr,ap, sf, vc, zvps, diff = to - from, tip},
   {sq, eq, sw, hs, hr, ap, sf, vc, zvps} = {ShaftQ, EdgesQ, ShaftWidth, HeadSize, HeadRatio, Aperture, ScaleFactor, VectorColor, ZeroVectorPointSize} /. {opts} /. Options[Vector];
   If[from==to, Return[Graphics3D[{PointSize[zvps],vc,Point[from]}]]];
   diff = If[NumberQ[N[sf]], sf diff, sf[diff]];
   tip = from + diff;
   If[NumberQ[hs], hr = hs/Norm[diff]];
   tr = hr/2;
   Graphics3D[{If[sq,{vc,Thickness[sw],Line[{from,tip-tr diff}]},{}], If[eq, {}, EdgeForm[]], {SurfaceColor[vc],
      With[{temp = Perpend[diff]},
      Polygon[Append[#,tip]]& /@ Partition[(tip - hr diff + #)&/@((ap hr Norm[diff] Base3D) . {temp,Normalize[Cross[temp,diff]]}),2,1]]}}]
] /; Length[from]==Length[to]==3


Arrow[a_, Tail->b_, opts___] := Vector[b,a+b,opts]

Arrow[a_, opts___] := Vector[Table[0,{Length[a]}],a,opts]

ArrowHead[a_, b_, opts___] := Vector[a-b,a,ShaftQ->False,opts];

End[];

EndPackage[];
Null

Graphics`Colors`GosiaGreen = 
        RGBColor[0,  0.392187,  0];

ThreeAxes[u,v]

ThreeAxes[u_, v_] := 
        Graphics3D[{{Blue, Line[{{-u, 0, 0}, {u, 0, 0}}]}, 
            Text["x", {u + v, 0, 0}], {Blue, 
              Line[{{0, -u, 0}, {0, u, 0}}]}, 
            Text["y", {0, u + v, 0}], {Blue, 
              Line[{{0, 0, 0}, {0, 0, u}}]}, 
            Text["z", {0, 0, u + v}]}];
ThreeAxes[u_] := ThreeAxes[u, (u/8)];
ThreeAxes::"usage" = "ThreeAxes[a,b] makes a standard cartesian \
axis graphics object with x, y, and z running from -a to a, and with axis \
labels b units beyond the tips of the axes.  ThreeAxes[a] is \
ThreeAxes[a,a/8].";

B.1) Rank one matrices.
 Using SVD to build matrices with sums of well-chosen rank one matrices

B.1.a.i)  A sample rank one matrix hitting on 2D and hanging in 3D

Give a sample of a rank one matrix hitting on 2D and hanging in 3D.
Show off the action of its hits.

Here's a rank one matrix hitting on 2D and hanging in 3D

A = (({{1.536, -1.024}, {1.92, -1.28}, {-1.728, 1.152}}));
MatrixForm[
      A]

Check:

rank= Length[SingularValues[A][[2]]]

See what a hit with this matrix does to these random points in 2D:

a = 2;
ranger  =   a;
points = 
        Table[{Random[Real, {-a, a}], Random[Real, {-a, a}]}, {k, 1, 
            600}];
pointcolor[k_] = 
        RGBColor[0.5 (Sin[((2 π k)/Length)[points]] + 1), 
          0.5 (Cos[((2 π k)/Length)[points]] + 1), 
          0.3];
pointplot = 
        Show[Table[
            Graphics[{PointSize[0.035], pointcolor[k], 
                Point[points[[
                    k]]]}], {k, 1, Length[points]}], 
          Axes -> True, AxesLabel -> {"x", "y"}, 
          PlotRange -> {{-ranger, ranger}, {-ranger, ranger}}, 
          PlotLabel -> "Before the hit with rank 1 matrix"];

Now hit all those points into 3D with A:

hitpointplot=Show[Table[Graphics3D[{PointSize[0.035],pointcolor[k],Point[A.points[[
                            k]]]}],{k, 1, Length[points]}],Axes -> True,AxesLabel -> {"x", "y", "z"},PlotRange -> All,PlotLabel -> "After the hit with rank 1 matrix",Axes -> True,ViewPoint -> CMView];

The rank 1 matrix A hits all those 2D points onto one 3D line - just as you would expect of a rank one matrix.

B.1.a.ii)  Making a rank one matrix

Make a rank one matrix A hitting on 6D and hanging in 4D

Go back to the SVD ingredients
          A = hanger.stretcher.aligner

Because A is a rank 1 , the SVD of A consists of 
▪ one 4D hangerframe (unit) vector
▪ one stretch factor,
▪ one 6D alignerframe (unit)vector
Make one of each and build the resulting hanger, stretcher and aligner:
▪The hanger:

hitdim = 6;
hangdim = 4;
hangerframevector = Normalize[Table[Random[Real, {-2, 2}],{i,1,hangdim}]];
hanger = Transpose[{hangerframevector}];
MatrixForm[
      hanger]

▪The stretcher:

Clear[stretch];
stretch =  
        Random[Real, {0.5, 4}];
stretcher= DiagonalMatrix[{stretch}];
MatrixForm[stretcher]

▪The aligner:

alignerframevector = Normalize[Table[Random[Real, {-2, 2}],{i,1,hitdim}]];
aligner ={alignerframevector};
MatrixForm[
      aligner]

Here is the resulting rank 1 matrix:

A = hanger.stretcher.aligner;
MatrixForm[A]

Check it out:

rank = Length[SingularValues[A][2]]

Had it all the way.

B.1.a.iii)  The action of a rank one matrix A is
              A.X = stretch (X.alignerframevector) hangerframevector

Here's what you do to build a random rank one matrix hitting on 5D and hanging in 7D:

hitdim = 5;
hangdim = 7;
hangerframevector = Normalize[Table[Random[Real, {-2, 2}],{i,1,hangdim}]];
hanger = Transpose[{hangerframevector}];
MatrixForm[
      hanger]

▪The stretcher:

Clear[stretch];
stretch =  
        Random[Real, {0.5, 4}];
stretcher= DiagonalMatrix[{stretch}];
MatrixForm[stretcher]

▪The aligner:

alignerframevector = Normalize[Table[Random[Real, {-2, 2}],{i,1,hitdim}]];
aligner ={alignerframevector};
MatrixForm[
      aligner]

The resulting rank 1 matrix is:

A = hanger.stretcher.aligner;
MatrixForm[A]

Check it out:

rank = Length[SingularValues[A][2]]

Talk about the action of A in terms of the hangerframe vector, the alignerframe vector and the stretch factor.

When you hit  A = hanger.stretcher.aligner on a 5D vector X, you get
                      A.X  =  stretch (alignerframevector.X) hangerframevector;
 Check this out:

Clear[x];
X  =  {x[1], x[2], x[3], x[4], x[5]};
A . X

Compare:

Expand[stretch(alignerframevector . X) hangerframevector]

B.1.b)  SVD tells you that any matrix is the sum of rank 1 matrices

Old matrix hands know that you can use SVD to see that any matrix is the sum of well-chosen rank 1 matrices.
Demonstrate what they mean by going with this random matrix:

hitdim =  Random[Integer, {3, 5}];
hangdim  =  Random[Integer, {4, 7}];
A  =  
        Table[Random[Real, {-4, 4}], {i, 1, hangdim}, {j, 1, 
            hitdim}];
MatrixForm[A]

And duplicating it with a sum of rank one matrices.

Do SVD analysis of A:

Clear[alignerframe, stretch, hangerframe, 
        k];
rank = Length[SingularValues[A][[2]]];
hangerframe[k_]:=SingularValues[A][[1,k]];
ColumnForm[
      Table[hangerframe[k], {k, 1, rank}]]

stretch[k_]:=SingularValues[A][[2,k]];
Table[
      stretch[k], {k, 1, rank}]

alignerframe[k_]:=SingularValues[A][[3,k]];
ColumnForm[
      Table[alignerframe[k], {k, 1, rank}]]

For each k with 1 ≤ k ≤ rank, build the rank 1 matrix rankOne[k] whose 
        ▪ hangerframe vector is hangerframe[k] 
        ▪  stretch factor is stretch[k]  and
        ▪  alignerframe vector is alignerframe[k].
        
 Here are rankOne[k] for k = 1,2,...,rank::

Clear[hanger, stretcher, aligner, rankOne]
hanger[k_] := Transpose[{hangerframe[k]}];
stretcher[k_]:= DiagonalMatrix[{stretch[k]}];
aligner[k_] :={alignerframe[k]};
rankOne[k_] := hanger[k].stretcher[k].aligner[k];

ColumnForm[
      Table[MatrixForm[rankOne[k]], {k, 1, rank}]]

These are all rank 1 matrices and their sum is:

sum  = (∑_(k = 1)^rank) rankOne[
            k];
MatrixForm[sum]

Compare:

MatrixForm[A]

Yessiree Bob.

A = (∑_(k = 1)^rank) rankOne[k]

Try it for other random matrices A:

hitdim =  Random[Integer, {3, 5}];
hangdim  =  Random[Integer, {4, 7}];
A  =  
        Table[Random[Real, {-4, 4}], {i, 1, hangdim}, {j, 1, 
            hitdim}];
MatrixForm[A]

Clear[alignerframe, stretch, hangerframe, k];
rank = Length[SingularValues[A][[2]]];
hangerframe[k_]:=SingularValues[A][[1,k]];
stretch[k_]:=SingularValues[A][[2,k]];
alignerframe[k_]:=SingularValues[A][[3,k]];

Clear[hanger, 
      stretcher, aligner, rankOne]
hanger[k_] := Transpose[{hangerframe[k]}];
stretcher[k_]:= DiagonalMatrix[{stretch[k]}];
aligner[k_] :={alignerframe[k]};
rankOne[k_] := hanger[k].stretcher[k].aligner[k];

ColumnForm[
      Table[MatrixForm[rankOne[k]], {k, 1, rank}]]

These are all rank 1 matrices and their sum is:

sum  = (∑_(k = 1)^rank) rankOne[
            k];
MatrixForm[sum]

Compare:

MatrixForm[A]

It works that way every time.

B.2) Image and data compression:
Creative rounding of matrices by zeroing out insignificant stretch factors

B.2.a.i) Data matrices as images.
Using SVD and rank 1 matrices to compress an image by zeroing out insignificant stretch factors

From the really great book by
David Kahaner, Cleve Moler and Stephen Nash - Numerical Methods and Software, Prentice-Hall, 1989:

If getting in on the ground floor of this topic interests you, go on.

Start by taking a look at this graphic which is based on a 2D matrix A whose entries are random numbers between 0.0 and 1.0:

dim = 2;
\(\(A = 
        Table[Random[Real, {0, 1}], {i, 1, dim}, {j, 1, dim}];\)
\)
\(\(jump = 1\/dim;\)
\)
Clear[box, boxlabel, x, y];
\(\(box[matrix_, x_, y_] := 
        Graphics[{GrayLevel[
              matrix[[dim\ y + 1, 
                dim\ x + 1]]], 
            Polygon[{{x, 1 - y}, {x + jump, 1 - y}, {x + jump, 
                  1 - y - jump}, {x, 1 - y - jump}, {x, 
                  1 - y}}]}];\)
\)
\(\(boxlabel[matrix_, i_, j_] := 
        Graphics[{LightCadmiumRed, 
            Text[matrix[[dim\ y + 1, 
                dim\ x + 1]], {x + 0.5\ jump, 
                1 - y - 0.5\ jump}]}];\)
\)
Show[
        Table[{box[A, x, y], boxlabel[A, x, y]}, {x, 0, 1 - jump, jump}, {y, 
            0, 1 - jump, jump}]];

Compare:

MatrixForm[A]

Kick up the dimension:

\(\(dim = 3;\)
\)
\(\(A = 
        Table[Random[Real, {0, 1}], {i, 1, dim}, {j, 1, dim}];\)
\)
\(\(jump = 1\/dim;\)
\)
Show[
        Table[{box[A, x, y], boxlabel[A, x, y]}, {x, 0, 1 - jump, jump}, {y, 
            0, 1 - jump, jump}]];

Compare:

MatrixForm[A]

Again:

\(\(dim = 4;\)
\)
\(\(A = 
        Table[Random[Real, {0, 1}], {i, 1, dim}, {j, 1, dim}];\)
\)
\(\(jump = 1\/dim;\)
\)
Show[
        Table[{box[A, x, y], boxlabel[A, x, y]}, {x, 0, 1 - jump, jump}, {y, 
            0, 1 - jump, jump}]];

Compare:

MatrixForm[A]

Get it?
Now look at this graphic which is based on a 16D matrix A whose entries are numbers between 0.0 and 1.0:

\(\(dim = 16;\)
\)
\(\(A = 
        Table[N[\(1 + Sin[1\/4\ π\ \((i + j)\)] + Random[Real, {0, 0.1}]\)\/2.5], {i, 1, dim}, {j, 1, dim}];\)
\)
\(\(jump = 1\/dim;\)
\)
original = 
        Show[Table[
            box[A, x, y], {x, 0, 1 - jump, jump}, {y, 0, 1 - jump, jump}], 
          PlotLabel -> "Original Image"];

To reproduce this exactly, you need to keep track of all of the 256 (= 16 × 16)  entries of A.
How do image compression folks use SVD go about trying to nearly reproduce this image keeping track of much less data?

Set up the rank one matrices whose sum is A:

Clear[alignerframe, stretch, hangerframe, k];
rank = Length[SingularValues[A][[2]]];
hangerframe[k_]:=SingularValues[A][[1,k]];
stretch[k_]:=SingularValues[A][[2,k]];
alignerframe[k_]:=SingularValues[A][[3,k]];

Clear[hanger, 
      stretcher, aligner, rankOne]
hanger[k_] := Transpose[{hangerframe[k]}];
stretcher[k_]:= DiagonalMatrix[{stretch[k]}];
aligner[k_] :={alignerframe[k]};
rankOne[k_] := hanger[k].stretcher[k].aligner[k];

At this stage you know
                A  = (∑_(k = 1)^rank) rankOne[k]
Trot out your SVD analysis of A and look at the SVD stretch factors for A:

SingularValues[A][[2]]

The only SVD stretch factors that pack much punch are the first three.

Data compression folks use this fact to creatively roundoff A so that only the first three significant SVD stretch factors are preserved and the others are zeroed out.

roundedA = (∑_(k = 1)^3)rankOne[k];
SingularValues[roundedA][[2]]

Then they plot the resulting image:

\(\(\(compressed\  = 
        Show[Table[
            box[roundedA, x, y], {x, 0, 1 - jump, jump}, {y, 0, 1 - jump, 
              jump}], PlotLabel -> "Compressed Image"];\)\(
\)\)\)

Cool.
Compare:

Show[original];

If anything, the compressed image is cleaner than the original!

B.2.a.ii) Why the compressed image can be sent with less data

Here are the two images from part i):

Show[GraphicsArray[{original, compressed}]];

Explain in detail why the compressed image can be communicated by sending much less data than is needed to communicate the original.

To send the data for the compressed image, you send only:

roundedA = (∑_(k = 1)^3)rankOne[k];

where rankOne[k] = hanger[k].stretcher[k].aligner[k]

To communicate this roundedA you need to send these data:

hangerdata = Table[hangerframe[k],{k,1,3}]

stretchdata  =Table[stretch[k],{k,1,3}]

alignerdata = Table[alignerframe[k],{k,1,3}]

The total count of the individual numbers needed to communicate roundedA  is:

totalcompressed  =  
      Length[Flatten[hangerdata]] +  Length[stretchdata]  +  
        Length[Flatten[alignerdata]]

The total count of the individual numbers needed to communicate the original data matrix A is

totaloriginal  =  16×16

The relative savings in decimals:

N[(totalcompressed/totaloriginal)]

You can send the cleaner compressed image with only about 39 % of the total data count needed to send the original.

B.2.a.iii) Why you don't want to multiply out  
           rankOne[k] = Transpose[{hangerframe[k]}].{stretch[k]}.{aligner[k]} 
before you send them out

Above it says that you would tell the recipient of the data to multiply everything out to get
                 rankOne[k] = Transpose[{hangerframe[k]}].{stretch[k]}.{aligner[k]} 
 and then put

roundedA = (∑_(k = 1)^3)rankOne[k]

Why shouldn't you be nice enough to multiply everything out before you send the data?

Look at what happens when you multiply out:
           rankOne[1] , rankOne[2] and rankone[3]

MatrixForm[rankOne[1]]

MatrixForm[rankOne[2]]

MatrixForm[rankOne[3]]

These three matrices each have as many entries as the original matrix A. So if you send them multiplied out, you will have to send 3 times as many numbers as need to communicate the original data matrix:

totalmultipliedout = (∑_(k = 1)^3)Length[
          Flatten[rankOne[k]]]

N[((totalmultipliedout )/totaloriginal)]

Multiplying out before sending is a big step backwards,

B.2.a.iv) Seeing the image emerge

Look at:

Table[Show[Table[box[ (∑_(k = 1)^p)rankOne[k],x,y],{x, 0, 1 - jump, jump},{y, 0, 1 - jump, jump}],PlotLabel -> p],{p,1,4}];

What do these depict?

The plot labeled "1" shows what image compression you get by rounding off A preserving only the first SVD stretch factor and zeroing out all the others.

p = 1;
Show[
        Table[box[rankOne[p], x, y], {x, 0, 1 - jump, jump}, {y, 0, 1 - jump, 
            jump}], PlotLabel -> p];

Terrible.
The plot labeled "2" shows what image compression you get by rounding off A preserving only the first two SVD stretch factors and zeroing out all the others:

p = 2;
Show[Table[box[(∑_(k = 1)^p)rankOne[k],x,y],{x, 0, 1 - jump, jump},{y, 0, 1 - jump, jump}],PlotLabel -> p];

Shaping up a bit.
The plot labeled "3" shows what image compression you get by rounding off A preserving only the first three SVD stretch factors and zeroing out all the others:

p = 3;
Show[Table[box[(∑_(k = 1)^p)rankOne[k],x,y],{x, 0, 1 - jump, jump},{y, 0, 1 - jump, jump}],PlotLabel -> p];

Excellent.
The plot labeled "4" shows what image compression you get by rounding off A preserving only the first four SVD stretch factors and zeroing out all the others:

p = 4;
Show[Table[box[(∑_(k = 1)^p)rankOne[k],x,y],{x, 0, 1 - jump, jump},{y, 0, 1 - jump, jump}],PlotLabel -> p];

Maybe that's overkill.