MGM.06 Beyond 3D
 LITERACY

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].";

Needs["Statistics`ContinuousDistributions`"];

Clear[matrixplucker];
matrixplucker[expr_, 
        vars_] := Simplify[
          Outer[D, {expr}, 
            vars]][[1]]

Here are the basic facts:

◼ The rank of a matrix is the number of non-zero SVD stretch factors for A. 

   This is also the number of hangerframe vectors for A.

◼ A given matrix A is of full rank if the rank of A = hitdim.

◼ Given a linear system A.X = Y,  put 
          Ytest = (∑_(k = 1)^rank)(Y . 
                  hangerframe[k]) hangerframe[k].

◼ If  Ytest = Y, and A is of full rank, then there is exactly one solution of A . X = Y and that solution is
          Xsol=(∑_(k = 1)^rank)(((Y . hangerframe[k])/(stretch[k]))) alignerframe[k]       (exactly determined).

◼ If  Ytest = Y, and A is of not of full rank, then 
          Xsol=(∑_(k = 1)^rank)(((Y . hangerframe[k])/(stretch[k]))) alignerframe[k] 
 is just one of infinitely many solutions  A . X = Y (under determined).

◼ If  Ytest = Y, then
          Xsol=(∑_(k = 1)^rank)(((Y . hangerframe[k])/(stretch[k]))) alignerframe[k] 
 is the solution of  A . X = Y of minimum norm.

◼ If Ytest ≠ Y, then the linear system A . X = Y has no solutions at all (over determined).

◼ If A is invertible, then
            (∑_(k = 1)^rank)(((Y . hangerframe[k])/(stretch[k]))) alignerframe[k]=(A^(-1)) . Y.

◼ For any matrix A (invertible or not)
            (∑_(k = 1)^rank)(((Y . hangerframe[k])/(stretch[k]))) alignerframe[k]=PseudoInverse[A] . Y.

◼ For any matrix A (invertible or not)
            Ytest =  (∑_(k = 1)^rank)(Y . hangerframe[k]) hangerframe[
              k]  
is as close to  Y as any other possible hit with A.

◼  For any matrix A (invertible or not)
A. (∑_(k = 1)^rank)(((Y . hangerframe[k])/(stretch[k]))) alignerframe[k]   = Ytest 
is as close to  Y as any other possible hit with A.

Indicate whether you agree or disagree with each of the following statements. 
You are not asked to explain your responses, but feel free to throw in an explanation if you want to.

L.a.i) Using  Ytest =  (∑_(k = 1)^rank)(Y . hangerframe[k]) hangerframe[
              k]

You are given a matrix A and use it to make a linear system 
A . X = Y.  You check the rank of A and learn A is of full rank. Then you do your SVD analysis of A and calculate               
             Ytest = (∑_(k = 1)^rank)(Y . 
                  hangerframe[k]) hangerframe[k]
and then you hit learn that
              Ytest = Y.
You now have enough information to announce that there is exactly one solution of the linear system A.X = Y.

L.a.ii) Using  Ytest =  (∑_(k = 1)^rank)(Y . hangerframe[k]) hangerframe[
              k]

You are given a matrix A and use it to make a linear system 
A . X = Y.  You check the rank of A and learn A is not of full rank. Then you do your SVD analysis of A and calculate               
             Ytest = (∑_(k = 1)^rank)(Y . 
                  hangerframe[k]) hangerframe[k]
and then you hit learn that
              Ytest = Y.
You now have enough information to announce that there is exactly one solution of the linear system A.X = Y.

L.a.iii) Using Ytest =  (∑_(k = 1)^rank)(Y . hangerframe[k]) hangerframe[
              k]

You are given a matrix A and use it to make a linear system 
A . X = Y. . Then you do your SVD analysis of A and calculate               
            Ytest =  (∑_(k = 1)^rank)(Y . hangerframe[k]) hangerframe[
              k]
and learn that
             Ytest ≠ Y.
You now have enough information to announce that there are no solutions of  the linear system A.X = Y.

L.d.vi) Using Ytest =  (∑_(k = 1)^rank)(Y . hangerframe[k]) hangerframe[
              k]

You are a given square A (so that hitdim = hangdim) and use it to make a linear system A.X = Y. 
 You do your SVD analysis of A and calculate               
       Ytest =  (∑_(k = 1)^rank)(Y . hangerframe[k]) hangerframe[
              k].
And then you find that
        Ytest  ≠  Y.
This tells you that
        Det[A] = 0.

L.b.i) Using rank

You are given a matrix A and use it to make a linear system A . X = Y.  You check the rank of A and learn A is of full rank. 
You now have enough information to announce that there is exactly one solution of the linear system A.X = Y.

L.b.ii) Xpprox and Ytest

You are given a matrix A and use it to make a linear system 
A . X = Y.  You do your SVD analysis of A and calculate  
            Xapprox = (∑_(k = 1)^rank)(((Y . hangerframe[k])/(stretch[k]))) alignerframe[k]

and

Ytest =  (∑_(k = 1)^rank)(Y . hangerframe[k]) hangerframe[
              k]
At this stage, you are guaranteed that 
             A.Xapprox = Ytest.

L.b.iii) Using  Xapprox = (∑_(k = 1)^rank)(((Y . hangerframe[k])/(stretch[k]))) alignerframe[k]

You are given a matrix A and use it to make a linear system 
A . X = Y.  You check the rank of A and learn A is not of full rank. Then you do your SVD analysis of A and calculate               
             \(Xapprox = ∑\+\(k = 1\)\%rank\((Y . hangerframe[k]\/stretch[k])\)\ alignerframe[k]\)
and then you hit Xapprox with A and learn that
            A.Xapprox = Y.
You now have enough information to announce that the linear system A.X = Y has lots of solutions.

L.b.iv) Using  Xapprox = (∑_(k = 1)^rank)(((Y . hangerframe[k])/(stretch[k]))) alignerframe[k]

You are given a matrix A and use it to make a linear system 
A . X = Y.  You check the rank of A and learn A is not of full rank. Then you do your SVD analysis of A and calculate               
             \(Xapprox = ∑\+\(k = 1\)\%rank\((Y . hangerframe[k]\/stretch[k])\)\ alignerframe[k]\)
and then you hit Xapprox with A and learn that
             A.Xapprox ≠ Y.
You now have enough information to announce that the linear system A.X = Y no solution at all.

L.c.i) Full rank

You are given a linear system 
         A.X = Y.
You check the rank of A and find that A is of full rank. 
This guarantees that the given linear system has exactly one solution (and no more).

L.c.ii) Rank deficient

You are given a linear system 
         A.X = Y.
You check the rank of A and find that A is not of full rank. 
This guarantees that the given linear system has no solution at all.

L.c.iii) Rows and rank

Saying that a matrix A is of full rank is the same as saying that the rank of A is equal to the number of horizontal rows of A.

L.c.iv) Columns and rank

Saying that a matrix A is of full rank is the same as saying that the rank of A is equal to the number of vertical columns of A.

L.c.v)  hitdim > hangdim

You are given a linear system 
         A.X = Y.
with hitdim > hangdim
It is automatic that the coefficient matrix A is not of full rank.

L.c.vi)  Full rank

You are given a linear system 
         A.X = Y.
You check the rank of A and find that A is of full rank. This tells you that the system has exactly one solution or no solutions at all.

L.c.vii)  Same number of equations as variable x[i]'s

You are given a linear system 
         A.X = Y
with the same number of equations as variable x[i]'s
You check the rank of A and find that A is of full rank. This tells you that the system has exactly one solution.

L.c.viii) Rank deficient

You are given a linear system 
            A.X = Y. 
You check the rank of A and find that A is not of full rank. 
This tells you that the system either has no solution or has many solutions.

L.c.ix) Rank and hangdim

Given a matrix hitting on hitdimD and hanging in hangdimD, you can be sure that 
that the rank of A cannot be bigger than hangdim.

L.c.x) Rank and hitdim

Given a matrix hitting on hitdimD and hanging in hangdimD, you can be sure that 
that the rank of A cannot be bigger than hitdim.

L.d.i) Many solutions?

You are a given matrix A that hits on 8D and hangs in 5D.   If the rank of A is 5, then no matter what Y in 5D you go with, you are guaranteed that the corresponding linear system A.X = Y has many solutions.

L.d.ii) Rank of A is the same as the rank of (A^t)?

Given any matrix A, you are guaranteed that the rank of A is the same as the rank of (A^t).

L.d.iii)  Det[A] ≠ 0

Given any square matrix A (so that hitdim = hangdim), saying that A is invertible is the same as saying that Det[A] ≠ 0.

L.d.iv) Rank > hangdim?

It is possible to encounter a matrix A that hangs its hits in 4D and whose rank is 6.

L.d.v) Stretches

It is possible to encounter a matrix that hits on 6D and hangs its hits in 4D and has five non-zero SVD stretch factors.

L.d.vi) Hanger frame

It is possible to encounter a matrix that hangs its hits in 4D and whose SVD hanger frame consists of five mutually perpendicular unit vectors.

L.d.v) Stretches

It is possible to encounter a matrix that hits on 6D and hangs its hits in 4D and has exactly four non-zero SVD stretch factors.

L.d.vi) Hanger frame

It is possible to encounter a matrix that hangs its hits in 7D whose SVD hanger frame consists of six mutually perpendicular unit vectors.

L.d.vii) hitdim < hangdim

You are given a linear system 
         A.X = Y.
with hitdim < hangdim
It is automatic that the coefficient matrix A is of full rank.

L.d.viii) Hitting on 5D and hanging in 8D

You are a given matrix A that hits on 5D and hangs in 8D.   Regardless of the rank of A, there are guaranteed to be many Y's in 8D for which the linear system A.X = Y has no solution.

L.e.i) Stretches and determinants

Here is a 3D matrix shown with its nonzero SVD stretch factors:

A = (({{0.5, 0, 1.1}, {1., -0.2, 1.5}, {3.5, -0.8, 4.9}}));
SingularValues[A][[2]]

This is enough to tell you that the rank of A is 2 and that the determinant of A is 0

L.e.ii)  Stretches and determinants

Here is a 4D matrix shown with its SVD stretch factors:

A = (({{1., -0.4, 0.8, -1. }, {2., -0.4, 3., 0}, {-0.5, 0, -1.1, -0.5}, {3.5, -0.8, 4.9, -0.5}}));
SingularValues[A][[2]]

This is enough to tell you that the rank of A is 2 and that the determinant of A is:

7.30167   1.22704

L.e.iii) Four linear equations in four variables

Here is a linear system  consisting of four equations involving four variables x[1],x[2],x[3] and x[4]:

(Subscript[a, 1]) x[1] + (Subscript[b, 1]) x[2] + (Subscript[c, 1]) x[3] + (Subscript[d, 1]) x[4]  =  y[1]
(Subscript[a, 2]) x[1] + (Subscript[b, 2]) x[2] + (Subscript[c, 2]) x[3] + (Subscript[d, 2]) x[4]  =  y[2]
(Subscript[a, 3]) x[1] + (Subscript[b, 3]) x[2] + (Subscript[c, 3]) x[3] + (Subscript[d, 3]) x[4]  =  y[3]
(Subscript[a, 4]) x[1] + (Subscript[b, 4]) x[2] + (Subscript[c, 4]) x[3] + (Subscript[d, 4]) x[4]  =  y[4]

Given specific values for all the constants

(Subscript[a, 1]), (Subscript[a, 2]), (Subscript[a, 3]), (Subscript[a, 4]),  (Subscript[b, 1]),(Subscript[b, 2]),(Subscript[b, 3]),(Subscript[b, 4]),(Subscript[c, 1]),(Subscript[c, 2]),(Subscript[c, 3]),(Subscript[c, 4]),(Subscript[d, 1]),(Subscript[d, 2]),(Subscript[d, 3]),(Subscript[d, 4]) and y[1], y[2], y[3], y[4],

you are guaranteed that the resulting linear system has exactly one solution for each of
        x[1], x[2], x[3], x[4]
because you have four equations in four unknowns.

L.e.iv) Det[A] ≠ 0

You are given a square matrix A (so that hitdim = hangdim) A and use it to make a linear system A.X = Y.  You calculate Det[A] and find that Det[A] ≠ 0. You know have enough information to announce that the linear system has exactly one solution.

L.f.i) Approximate solution

You are given a matrix A and use it to make a linear system 
A . X = Y.  You do your SVD analysis of A and calculate               
             \(Xapprox = ∑\+\(k = 1\)\%rank\((Y . hangerframe[k]\/stretch[k])\)\ alignerframe[k]\)
and then you announce that Xapprox the best approximate solution you can get.

L.f.ii)  PseudoInverse

You are given a matrix A that hits on 10D and hangs in 13D.   Take a Y in 13D , do your SVD analysis of A and you  calculate               
             \(∑\+\(k = 1\)\%rank\((Y . hangerframe[k]\/stretch[k])\)\ alignerframe[k]\)
and then you announce that  this is the same as
             PseudoInverse[A] . Y.

L.f.iii) PseudoInverse

You are given an invertible matrix A that hits on 9D and hangs in 9D.   You take a Y in 9D , do your SVD analysis of A and you  calculate               
             \(∑\+\(k = 1\)\%rank\((Y . hangerframe[k]\/stretch[k])\)\ alignerframe[k]\).
And then you announce that  this is the same as
             PseudoInverse[A] . Y.

L.f.iv) Inverse and PseudoInverse

If A is an invertible matrix, then
             (A^(-1))  =  PseudoInverse[A].

L.f.v) PseudoInverse

A square matrix has in inverse if its determinant is not 0, but all other matrices have a pseudo inverse.

L.f.vi) Ytest and invertibility

You are a given square A (so that hitdim = hangdim) and use it to make a linear system A.X = Y. 
 You do your SVD analysis of A and calculate               
        Ytest =   (∑_(k = 1)^rank)(Y . 
                hangerframe[k]) hangerframe[k].
And then you find that
        Ytest  ≠  Y.
This tells you that A is not invertible.

L.f.vii) Solution of smallest norm

You are given a matrix A and use it to make a linear system 
A . X = Y.  You check the rank of A and learn A is not of full rank. Then you do your SVD analysis of A and calculate               
             \(Xapprox = ∑\+\(k = 1\)\%rank\((Y . hangerframe[k]\/stretch[k])\)\ alignerframe[k]\)
and then you hit Xapprox with A and learn that
             A.Xapprox= Y.
You now have enough information to announce that the linear system A.X = Y is underdetermined (lots of solutions) and to announce that if Xother is any solution, then
       (Sqrt[Xtest . Xtest]) < (Sqrt[Xother . Xother])