MGM.06 Beyond 3D
TUTORIALS

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]]

T.1) Quick first impressions based on comparing hitdim and hangdim
If hitdim ≠ hangdim, then A is not invertible

T.1.a.i) What to look for when hitdim > hangdim (i.e. fewer equations than unknown x[k]'s)

What are the possibilities when you go with a linear system
A.X = Y and hitdim > hangdim (fewer equations than unknown x[k]'s)?

In the case of fewer equations than unknown x[k]'s, most folks react initially ( because there are fewer constraints than degrees of freedom) with the view that the linear system A.X = Y is likely to have more than one solution (underdetermined).  This is a real possibility:

Possibility 1: More than one solution if  Ytest  =  Y.

Reason: 
If  Ytest  =  Y., then A.Xsol = Y.
Because hitdim > hangdim, hits with A have to squash at least 
        ( hitdim - hangdim) dimensions of hitdimD 
onto 
         {0, 0, . . ., 0}. 
This tells you that if you take any nonzero vector XX with
             A . XX = {0, 0, . . ., 0},
then 
             A.(Xsol + XX) = A.Xsol + A.XX = A.Xsol + {0,0,. . .,0} = Y.
            
This gives you more than one solution.

Possibility 2: No solution when Ytest  ≠  Y .

There are no other possibilities.

T.1.a.ii) What to look for when hitdim = hangdim  (same number of equations as unknown x[k]'s)

What are the possibilities when you have a linear system
        A.X = Y 
with 
        hitdim = hangdim (i.e. same number of equations as unknown x[k]'s)?

In the case of the same number equations than unknown x[k]'s, most folks react initially ( because there are the same number of constraints as degrees of freedom) with the view that the linear system A.X = Y is likely to have exactly one solution(exactly determined). This is not always correct. In fact, there are three possibilities:

Possibility 1: Exactly one solution when Det[A] ≠ 0

This does happen if
             Det[A] ≠ 0 (ensuring that A is of full rank).
In the case Det[A] ≠ 0, you don't need to check whether Ytest = Y because the lone solution is
             X = (A^(-1)).Y.

Possibility 2: More than one solution when Det[A] = 0

This happens when Ytest = Y but Det[A] = 0  (ensuring that rank of A < hit dim so that A is not of full rank).

Possibility 3: No solution when Det[A] = 0

This happens when Ytest ≠ Y.  In the case that Det[A] = 0, then you are guaranteed that A is not of full rank so hits with A cannot fill up all of hangdim.

T.1.a.iii) What to look for when hitdim < hangdim (i.e. more equations than unknown x[k]'s)

What are the possibilities when you have a linear system
         A.X = Y 
with 
         hitdim < hangdim (more equations than unknown x[k]'s)?

In the case of more equations than unknown x[k]'s, most folks react initially ( because there are more constraints than degrees of freedom) with the view that the linear system A.X = Y is likely to have no solutions at all.

This reaction is more-or-less natural because hits on hitdimD cannot fill up all of hangdimD. And this outcome is a possibility:

Possibility 1: No solution when Ytest ≠Y.

Possibility 2: Exactly one solution  when Ytest = Y and rank of A = hitdim (so that A is of full rank).

Possibility 3: More than one solution  when Ytesy = Y and rank of A < hitdim (so that A is not of full rank).

T.1.b.i) If hitdim > hangdim, then A is not invertible.

Here is a random matrix with hitdim > hangdim:

hitdim = 5; hangdim = 3;
Clear[a, i, j];
\(\(a[i_, j_] := Random[Real, {\(-2\), 2}];\)
\)
\(\(A = Table[a[i, j], {i, 1, hangdim}, {j, 1, hitdim}];\)
\)
MatrixForm[A]

Here's Mathematica's attempt at a calculation of (A^(-1)):

MatrixForm[Inverse[A]]

Garbage.
Explain how you could have known in advance that A is not invertible.

This matrix A hits 5D vectors into 3D.  
When you hit all of 5D with A, you squash at least two dimensions to {0,0,0}. 
You cannot undo the hits on the squashed dimensions.

T.1.b.ii) If hitdim < hangdim, then A is not invertible.

Here is a random matrix with hitdim < hangdim:

hitdim = 4; hangdim = 7;
Clear[a, i, j];
\(\(a[i_, j_] := Random[Real, {\(-2\), 2}];\)
\)
\(\(A = Table[a[i, j], {i, 1, hangdim}, {j, 1, hitdim}];\)
\)
MatrixForm[A]

Here's Mathematica's attempt at a calculation of (A^(-1)):

MatrixForm[Inverse[A]]

Useless.
Explain how you could have known in advance that A is not invertible.

This matrix A hits 4D vectors into 7D.  
When you hit all of 4D with A, you cover at most four dimensions of 7D. 
That leaves three full dimensions of 7D that are not accessible via hits with A.
When you try to undo A inside these three extra dimensions, you are up the old creek. 
And that's why A is not invertible.

T.2) The absolute value of determinant for square matrices A in higher dimensions. 
If Det[A] ≠ 0, then the corresponding linear system A . X = Y has exactly one solution (exactly determined). 
If Det[A] = 0, then the linear system A . X = Y either has no solution (overdetermined) or has many solutions (underdetermined)

T.2.a.i) For a square matrix, |Det[A]| is the product of A's stretch SVD factors

Look at these matrices that hit and hang in the same dimension:

Clear[a, i, j];
dim = 4;
\(a[i_, j_] := Random[Real, {\(-7\), 7}]
\)
\(\(A = Table[a[i, j], {i, 1, dim}, {j, 1, dim}];\)
\)
MatrixForm[A]

The rank of this matrix is:

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

Here is the product of the four stretch factors and a calculation of the absolute value of the determinant of A:

Clear[stretch, k];
stretch[
            k_] := SingularValues[
                A][[2]][[k]];

(∏_(k = 1)^rank)stretch[k]

Abs[Det[A]]

More:

Clear[a, i, j];
\(\(dim = Random[Integer, {4, 8}];\)
\)
\(a[i_, j_] := Random[Real, {\(-7\), 7}]
\)
\(\(A = Table[a[i, j], {i, 1, dim}, {j, 1, dim}];\)
\)
MatrixForm[A]

The rank of this matrix is:

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

Here is the product of the five stretch factors and a calculation of the absolute value of the determinant of A:

Clear[stretch, k];
stretch[
            k_] := SingularValues[
                A][[2]][[k]];

(∏_(k = 1)^rank)stretch[k]

Abs[Det[A]]

Another:

Clear[a, i, j];
\(\(dim = Random[Integer, {4, 7}];\)
\)
\(a[i_, j_] := Random[Real, {\(-7\), 7}]
\)
\(\(A = Table[a[i, j], {i, 1, dim}, {j, 1, dim}];\)
\)
MatrixForm[A]

The rank of this matrix is:

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

Here is the product of the five stretch factors and a calculation of the absolute value of the determinant of A:

Clear[stretch, k];
stretch[
            k_] := SingularValues[
                A][[2]][[k]];

(∏_(k = 1)^rank)stretch[k]

Abs[Det[A]]

What's the message?

The message is that even in higher dimensions, the absolute value of the determinant of a square matrix A is the same as the product of A's stretch factors. 

It's quite legal to take this as the definition of the absolute value of the determinant of A.

T.2.a.ii) Non-zero determinants and exactly determined linear systems

You are given a linear system with the same number of equations as unknowns. You take the coefficient matrix A and use Mathematica to calculate the determinant of A. 
It turns out to be non-zero. 
You announce that the linear system is exactly determined in the sense that it has exactly one solution. 
You are right. Why are you right?

Just for the sake of explanation, agree that you have 6 equations in 6 unknowns. This means the coefficient matrix A hits on 6D and hangs in 6D.
You are given
       |Det[A]| ≠ 0
and you know that
       |Det[A]| is the product of all of A's 6 stretch factors.
Because 
       |Det[A]| ≠ 0, 
you are left with the valid conclusion that none of A's stretch factors is 0. 
And this tells you that hits with A on 6D fill up all of 6D. 
This also tells you that 
       rank of A = hitdim = 6,
so that A is of full rank.
This means that no matter what 6D Y you go with, you are guaranteed that A.X = Y has at least one solution.
And because A is of full rank, there is exactly one solution.

T.2.a.iii) Zero determinants and trouble

You are given a linear system with the same number of equations as unknowns. You take the coefficient matrix A and use Mathematica to calculate the determinant of A. 
It turns out to be zero. 
You announce that the linear system is either under determined (has more than one solution) or is overdetermined (has no solution).
You are right. Why are you right?

Just for the sake of explanation, agree that you have 7 equations in 7 unknowns. This means the coefficient matrix A hits on 7D and hangs in 7D.
You are given
       |Det[A]| = 0
and you know that
       |Det[A]| is the product of A's 7 stretch factors.
Because 
       |Det[A]| = 0, 
you are left with the valid conclusion that at least one of A's stretch factors is 0. 
And this tells you that hits with A on 7D fill up at most six of the seven dimensions of 7D. 
This allows for plenty of Y's for which A.X = Y has no solution.

This also allows for plenty of Y's for which A.X = Y does have a solution. But in this case, the zero stretch factor tells you that
       rank of A < hitdim
so that A is not of full rank. This tells you that any one solution is accompanied by many more solutions.

T.2.b) Vandermonde matrices

Given distinct (different) numbers t[1], t[2], t[3],t[4] and t[5] make this matrix:

dim = 5;
Clear[t, x, y, i, j];
A=Table[((t[i])^(j - 1)),{i, 1, dim},{j, 1, dim}];

MatrixForm[A]

Use this matrix as the coefficient matrix to make this linear system:

X = Table[x[k], {k, 1, dim}];
\(\(Y = Table[y[k], {k, 1, dim}];\)
\)
ColumnForm[Thread[  A .  X ==  Y]]

Explain this:
This linear system is exactly determined in the sense that there is exactly one solution {x[1],y[2],x[3],x[4],x[5]} corresponding to any given {y[1],y[2],y[3],y[4],y[5]}.

Let Mathematica calculate the absolute value of the determinant of A in factored form:

Abs[Factor[Det[A]]]

Because the t[i]'s are all different, this is the product of non-zero numbers, so there is no way for this to be zero. 
This is enough to tell you that the linear system is exactly determined  in the sense that there is exactly one solution {x[1],y[2],x[3],x[4],x[5]} corresponding to any given {y[1],y[2],y[3],y[4],y[5]}.

T.3) Columns and rank; rows and rank: How to make matrices that are not of full rank

T.3.a.i) Columns and rank

Here is a matrix entered via its vertical columns:

Clear[column, j];
hitdim = 
        4;
column[1] = {0.0, 0.0, 3.3, 0.0, 
          0.2};
column[2] = {0.0, -2.0, -1.3, 1.0, 
          0.0};
column[3] = {1.4, 0.0, 0.0, 
          2.0, -3.0};
column[4] = 1.2 column[1] - 3.5 column[2];

A = Transpose[Table[column[j], {j, 1, 4}]];

MatrixForm[A]

Note that column[4] is a combination of some of other columns.
Now calculate the rank of this matrix:

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

This matrix hits on 4D; so 
         rank < hitdim 
and the matrix is not of full rank.
The parts below attempt to give an explanation about why this approach was guaranteed to produce a matrix not of full rank.

Agree that 
       {col[1], col[2], col[3], . . ., 
        col[hitdim]}
are the vertical columns of a given matrix A.
Explain this statement:
Saying that matrix A is not of full rank is the same as saying that there is a vector 
       {c[1], c[2], . . ., c[hitdim]} 
in hitdimD with at least one non-zero entry such that
        (∑_(i = 1)^hitdim) c[i] column[i] = {0, 
          0, . . ., 0}.

Saying that A is not of full rank is the same as saying that a hit with A squashes at least one non-zero vector in hitdimD onto the vector of all zeroes.

This is the same as saying there is a vector  
       {c[1], c[2], . . ., c[hitdim]} 
in hitdimD with at least one non-zero entry such that
       A . {c[1], c[2], . . ., 
            c[hitdim]} =  {0, 0,  ... , 0}.
       
Because
       A . {c[1], c[2], . . ., c[hitdim]}  =  (∑_(i = 1)^hitdim) c[i] column[i], 
 this is the same as saying that there is a vector 
           {c[1], c[2], . . ., c[hitdim]}
in hitdimD with at least one non-zero entry such that
       (∑_(i = 1)^hitdim) c[i] column[i] = {0, 
          0, . . ., 0}.

T.3.a.ii) More on columns and rank

Expound on this statement: 
Saying that matrix A is not of full rank is the same as saying that you can express at least one of the vertical columns of A in terms of some of the others.

Saying that matrix A is not of full rank is the same as saying that there is a vector 
       {c[1], c[2], . . ., c[hitdim]} 
in hitdimD with at least one non-zero entry such that
        (∑_(i = 1)^hitdim) c[i] column[i] = {0, 
          0, . . ., 0}.

Saying that matrix A is not of full rank is the same as saying that there is a vector 

Take (Subscript[i, 0]) with  c[(Subscript[i, 0])] ≠ 0 and note that 

        -c[(Subscript[i, 0])] column[(Subscript[i, 0])] 
        
           =    (∑_(i = 1)^hitdim) c[i] column[i] - c[(Subscript[i, 0])] column[(Subscript[i, 0])]
   
           = (∑_(i = 1)^((Subscript[i, 0]) - 1))c[i] column[i] +     (∑_(i = (Subscript[i, 0]) + 1)^hitdim) c[i] column[i].
This gives
      -c[(Subscript[i, 0])] column[(Subscript[i, 0])] = (∑_(i = 1)^((Subscript[i, 0]) - 1))c[i] column[i] +     (∑_(i = (Subscript[i, 0]) + 1)^hitdim) c[i] column[i].
so that all the influence of column[(Subscript[i, 0])] on the right is gone.
Divide through by  -c[(Subscript[i, 0])] to get

 column[(Subscript[i, 0])] = (∑_(i = 1)^((Subscript[i, 0]) - 1))(-c[i]/c[(Subscript[i, 0])]) 
          column[i] +     (∑_(i = (Subscript[i, 0]) + 1)^hitdim)(-c[i]/c[(Subscript[i, 0])]) 
          column[i]

 which exhibits column[(Subscript[i, 0])] as a combination of the other columns.

T.3.a.iii) Example: Columns and rank

Here is a matrix A entered via its vertical columns:

Clear[column, j];
hitdim = 5;
column[1] = {0.43, -1.32, -2.26, 1.47, 4.13, 
          2.32};
column[2] = {-0.34, 0.34, 0.0, 7.89, -2.96, 
          0.89};
column[3] = {2.58, 1.34, 3.82, 2.84, -1.44, 0.78};
column[4] = {0.61, 1.35, 1.91, 17.20, -6.64, 2.17};
\(\(column[5] = {1.72, \(-0.65\), \(-0.35\), 2.89, 3.41, 
          2.71};\)
\)
\(\(A = Transpose[Table[column[j], {j, 1, hitdim}]];\)
\)
MatrixForm[A]

Calculate the rank of this matrix:

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

This matrix hits on 5D; its rank is 3. 
So this matrix is not of full rank.
Exhibit at least one column that is a combination of at least some of the other columns.

Look at:

NullSpace[A]

Put:

Clear[c];
{c[1], c[2], c[3], c[4], 
        c[5]} = NullSpace[
          A][[1]]

This tells you:

A . {c[1], c[2], c[3], c[4], c[5]}

This is the same as:

(∑_(i = 1)^5)c[i] column[i]

This tells you how to exhibit column[3] as a combination of the other columns.

column[3]

(∑_(i = 1)^2)-((c[i] column[i])/(c[3]))+(∑_(i = 4)^5)-((c[i] column[i])/(c[3]))

It also tells you tells you how to exhibit column[1] as a combination of the other columns.

column[1]

(∑_(i = 2)^5)-((c[i] column[i])/(c[1]))

T.3.b) Rows and rank: How to make matrices that are not of full rank

Here is a matrix entered via its horizontal rows:

Clear[row, 
        j];
row[1] = {0.0, 0.0, 3.3, 0.0, 0.2};
row[
          2] = {0.0, -2.0, -1.3, 1.0, 0.0};
row[3] = {1.4, 
          0.0, 0.0, 2.0, -3.0};
row[4] = {2.4, -1.0, 1.0, 
          0.0, -1.0};
row[5] = 1.2 row[1] - 3.5 row[2];

A = Table[row[j], {j, 1, 5}];

MatrixForm[
      A]

Examine the code to see that horizontal 
        row[5] = 1.2 row[1] - 3.5 row[2]
 is a combination of some of other horizontal rows.
Explain why this fact guarantees that A is not of full rank.

Look at A:

MatrixForm[A]

Look at (A^t):

MatrixForm[Transpose[A]]

The horizontal rows of (A^t) are the vertical columns of A.
Because
     row[5] of A = 1.2 row[1] of A  - 3.5 row[2] of A
you know that
      column[5] of (A^t)=1.2  column[1]of (A^t) -3.5 column[2] of (A^t)
Part a) above tells you that 
         rank of (A^t) < 5 
 and because 
        rank of (A^t) = rank of A, 
you are guaranteed that
        rank of A < 5.

T.4) For under determined linear systems A.X = Y, 
Xapprox = PseudoInverse[A].Y is the solution of smallest norm

T.5.a) Xapprox = PseudoInverse[A].Y is the solution of smallest norm

Go with this linear system:

hitdim = 5; hangdim = 3;
Clear[a, i, j];
\(\(a[i_, j_] := Random[Real, {\(-2\), 2}];\)
\)
A = Table[a[i, j], {i, 1, hangdim}, {j, 1, hitdim}];
\(\(Y = {1.0, 0.0, 2.0};\)

\)
Clear[x, k];
\(\(X = Table[x[k], {k, 1, hitdim}];\)
\)
linearsystem = ColumnForm[Thread[A . X == Y]]

To look for solutions, determine the rank of A:

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

Compare with:

hitdim

rank < hitdim

Bad news. The coefficient matrix A is of not of full rank.
Now build
     \(Xapprox = ∑\+\(k = 1\)\%rank\((Y . hangerframe[k]\/stretch[k])\)\ alignerframe[k]\) = PseudoInverse[A].Y:

Xapprox=PseudoInverse[A].Y

Try it out:

A . Xapprox
Y

Xapprox is a solution and because A is of not of full rank, Xapprox is not the only solution.
 
Explain why Xapprox is the solution of smallest norm.
In other words, explain why
                ||Xapprox|| ⩽ ||any other solution||.

To get a hold of all the solutions, look at the alignerframe vectors that are squashed to {0,0,0} when they are hit with A:

NullSpace[A]

All the solutions are parameterized by:

Clear[s, t, solution];
solution[s_,t_]=Xapprox+s NullSpace[
              A][[1]]+ t NullSpace[
              A][[2]]

Try them out:

Expand[A . solution[s, t]]
Y

Good
Now look at
\(\(\(\(||\)\(\ \)\(\(solution[s, t]\)\( || \^2\)\)\)\(=\)\)\) solution[s, t] . solution[s, t]:

Expand[solution[s, t] . solution[s, t]]

Xapprox . Xapprox

This tells you that when you go with
     solution[s, t] = 
        Xapprox + 
          s NullSpace[A][[1]] + 
          t NullSpace[A][[2]],
you get
     \(\(\(||\)\(\ \)\(\(solution[s, t]\)\( || \^2\)\)\) = 
        Xapprox . Xapprox + s\^2 + t\^2\).
So the solution solution[s,t] with smallest norm is 
     solution[0, 0]=Xapprox .
The upshot:
    Xapprox
is the solution of smallest norm.

T.5) Plucking the coefficient matrix off a linear system

T.5.a) Going from linear formulas to matrices

Here's a common way of delivering a linear system:
First, you are given a list of linear formulas:

Clear[x];
\(\(linearformulas = 
{2.1\ x[1] + 2.3\ x[2] + 
            3.4\ x[3] - 5.4\ x[4] - 0.5\ x[5] - 2.5\ x[6], \(-1.5\)\ x[1] - 
            0.4\ x[2] + 13.5\ x[4] - 0.7\ x[5] + 3.7\ x[6], \(-0.7\)\ x[1] + 
            3.4\ x[3] + 4.5\ x[4] - 2.1\ x[5] + 
            4.2\ x[6], 
3.4\ x[2] - 0.6\ x[3] - 1.2\ x[4] - 
            0.6\ x[5], 
1.6\ x[1] + 2.8\ x[3] + 3.6\ x[4] + 
            2.1\ x[5] - 0.7\ x[6], 
2.3\ x[1] + 3.2\ x[3] + 
            4.9\ x[5]};\)
\)
ColumnForm[linearformulas]

Then you are given values 
             y[1], . . . ,y[6],
one for each linear formula, and you set 
             linearformula[k] = y[k]:

Clear[y];
Y = {2.4, -0.4, 2.7, -1.9, 0.8, 1.6};
\(\(y[k_] := Y[[k]];\)
\)
\(linearsystem\  = \ 

      ColumnForm[
        Thread[Table[
            linearformulas[\([k]\)]\  == \ y[k], {k, 1, 
              Length[linearformulas]}]]]\)

Once you have the coefficient matrix A for this linear system, then going about solving the linear system is not bad at all.
But reading off and manually typing A into Mathematica code is a real pain in the butt.
How can you get Mathematica  pluck out that coefficient matrix for you?

Here's one way:

xvariables = Table[x[k], {k, 1, 6}];
A=Outer[D, {linearformulas}, 
                xvariables][[1]];

MatrixForm[
      A]

Check it out:

hitdim = 6; hangdim = 6;
X = Table[x[k], {k, 1, hitdim}];
\(\(Y = Table[y[k], {k, 1, hangdim}];\)
\)
ColumnForm[Thread[A . X  ==   Y]]

Compare:

linearsystem

Worked like a charm. 
The solution of this linear system is:

X = Inverse[A] . Y

Check:

Y
A . X

Had it all the way.

T.5.a.ii) The weird instruction Outer[D,{linearformulas}, xvariables][[1]];

Go back to:

ColumnForm[linearformulas]

xvariables = Table[x[k], {k, 1, 6}];
\(\(A = \(Outer[D, {linearformulas}, 
            xvariables]\)[[1]];\)
\)
MatrixForm[A]

How did that weird instruction deliver the goods?

The instruction 
      Outer[D,{linearformulas},xvariables][[1]]
just took each linear formula and replaced it by its gradient.

T.6) Tips on entering matrices in Mathematica

T.6.a.i) Entering big matrices in Mathematica by rows

Here is a matrix

A =  (({{0.43`, -0.34`, 2.58`}, {-1.32`, 0.34`, 1.34`}, {-2.26`, 0, 3.82`}, {1.47`, 7.89`, 2.84`}, {4.13`, -2.96`, -1.44`}, {2.32`, 0.89`, 0.78`}}))

How do you enter this matrix in Mathematica by rows?

Here's one way:

Clear[row];
row[1] = {0.43, 0.34, 2.58};
row[2] ={-1.32,  0.34, 1.34};
row[3] = {-2.26, 0, 3.82};
row[4] = {1.47, 7.89, 2.84};
row[5] = {4.13, -2.96, -1.44};
row[6] = {2.32, 0.89, 0.78};
A = Table[row[k],{k,1,6}];
MatrixForm[A]

Compare:

A = (({{0.43`, -0.34`, 2.58`}, {-1.32`, 0.34`, 1.34`}, {-2.26`, 0, 3.82`}, {1.47`, 7.89`, 2.84`}, {4.13`, -2.96`, -1.44`}, {2.32`, 0.89`, 0.78`}}))

T.6.a.ii) Entering big matrices in Mathematica by columns

Here is a matrix

A = (({{0.43`, -0.34`, 2.58`}, {-1.32`, 0.34`, 1.34`}, {-2.26`, 0, 3.82`}, {1.47`, 7.89`, 2.84`}, {4.13`, -2.96`, -1.44`}, {2.32`, 0.89`, 0.78`}}))

How do you enter this matrix in Mathematica by columns?

Here's one way:

Clear[col];
col[1] = {0.43, -1.32, -2.26, 1.47, 4.13, 2.32};
col[2] ={-0.34, 0.34, 0, 7.89, -2.96, 0.89};
col[3] = {2.58, 1.34, 3.82, 2.84, -1.44, 0.78};

A =Transpose[ Table[col[k],{k,1,3}]];
MatrixForm[A]

Compare:

A = (({{0.43`, -0.34`, 2.58`}, {-1.32`, 0.34`, 1.34`}, {-2.26`, 0, 3.82`}, {1.47`, 7.89`, 2.84`}, {4.13`, -2.96`, -1.44`}, {2.32`, 0.89`, 0.78`}}))

T.6.a.iii) Entering big matrices in Mathematica by template

Here is a matrix

A = (({{0.43`, -0.34`, 2.58`}, {-1.32`, 0.34`, 1.34`}, {-2.26`, 0, 3.82`}, {1.47`, 7.89`, 2.84`}, {4.13`, -2.96`, -1.44`}, {2.32`, 0.89`, 0.78`}}))

How do you enter this matrix in Mathematica by template?

Here's one way:
Notice that A hits on 3D and hangs in 6D:

Go to the Input menu and select Create Table/Matrix Palette.
Check the Matrix button and select 6 rows and 3 columns and hit the OK button.
The result will be this:
                               (({{□, □, □}, {□, □, □}, {□, □, □}, {□, □, □}, {□, □, □}, {□, □, □}}))

Now type in the entries and you are off and running.

T.7) Square matrices as spreadsheets: Gaussian elimination and reduced row echelon form.  How it works and when you might want to use it

T.7.a.i) The basic idea

What is the basic idea behind Gaussian elimination?

You probably met Gaussian elimination way back in junior high school.

The idea is that when you have a system of linear equations, and you add non-zero multiples of one equation to another, then you change the linear system but do not change its solutions.

For instance if your linear system is

(({{1, 1}, {2, 4}})).{x,y} = {2, 10},

Add -2 times the first equation to the second equation and keep the first equation the same to get

(({{1, 1}, {2 -  2, 4 -  2}})).{x,y} = {2, 10 - 4},

This is the same as

(({{1, 1}, {0, 2}})).{x,y} = {2, 6},

And now you add -(1/2)  times the second equation to the first and keep the second equation the same to get

(({{1, 1  -  1}, {0, 2}})).{x,y} = {2 -3, 6},

This is the same as

(({{1, 0}, {0, 2}})).{x,y} = {-1,6},

Divide 2 through the second equation to get

(({{1, 0}, {0, 1}})).{x,y} = {-1 , 3},

At this point you read off the solution {x,y} = {-1,3}.
Try it out

(({{1, 1}, {2, 4}})).{-1, 3} =={2,10}

T.7.a.ii)  The augmented matrix for the linear system 
                 (({{1, 1}, {2, 4}})).{x,y} = {2, 10}  is  (({{1, 1, 2}, {2, 4, 10}})).
Spread sheet operations on the rows the augmented matrix give you the solution via reduced row echelon form

Stay with the same linear system

(({{1, 1}, {2, 4}})).{x,y} = {2, 10},

The augmented matrix for this linear system is :(({{1, 1, 2}, {2, 4, 10}}))

Use spread sheet operations on the rows of the augmented matrix to come up with the solution to the linear system.

You do the same thing you did above, but this time you do it to the rows of the augmented matrix.

Go with the augmented matrix  
                      (({{1, 1, 2}, {2, 4, 10}}))

Add -2 times the first row to the second row and keep the first row the same to get

(({{1, 1, 2}, {2-2, 4-2, 10 -4}})) = (({{1, 1, 2}, {0, 2, 6}}))

And now you add -(1/2)  times the second row to the first and keep the second row the same to get

(({{1-0, 1 -1, 2-3}, {0, 2, 6}})) = (({{1, 0, -1}, {0, 2, 6}}))

And then you divide the second row by 2 to get

(({{1, 0, -1}, {0, (2/2), (6/2)}})) = (({{1, 0, -1}, {0, 1, 3}}))

And read off the solution {x,y} = {-1,3}.

The beauty of it is that Mathematica will do this for you at a flip of the finger:

MatrixForm[RowReduce[ (({{1, 1, 2}, {2, 4, 10}}))]]

T.7.a.iii ) A 3D example powered by Mathematica

Here's a linear system:

A =(({{2.1, 2.3, 4.2}, {4.5, -0.6, 3.1}, {2.1, 1.9, -2.8}}));
dim = 3;
Clear[x, j];
X = 
        Table[x[j], {j, 1, dim}];
Y =  {1.2, 2.5, 
          3.0};
linearsystem = 
      ColumnForm[Thread[A . X == Y]]

Use the reduced row echelon form of the augmented matrix to solve this system.

The linear system is:

linearsystem

The augmented matrix is

augmentedA  = AppendRows[A, Transpose[{Y}]];
MatrixForm[augmentedA]

And look at this:

rowechelon  =  RowReduce[augmentedA];
MatrixForm[rowechelon]

Pick out the fourth column and put:

X =Table[rowechelon[[j,4]],{j,1,3}]

Check:

A . X

Y

Got it.

Go with:

augmented  = 
        AppendRows[A, Transpose[{Y}]];
MatrixForm[augmented]

Subtract appropriate multiples of the first row from the other rows:

pivot1 =  {augmented[[1]] , 
          augmented[[2]]  - 
            augmented[[2, 1]]/augmented[[1, 1]] augmented[[1]], 
          augmented[[3]]  - 
            augmented[[3, 1]]/augmented[[1, 1]] augmented[[1]]};

MatrixForm[pivot1]

Subtract  the  appropriate multiple of the second row from the third row:

pivot2  =  {pivot1[[1]] , pivot1[[2]] , 
          pivot1[[3]]  - 
            pivot1[[3, 2]]/pivot1[[2, 2]] pivot1[[2]]};
MatrixForm[pivot2]

Subtract the appropriate multiples of the third from the first and second rows:

pivot3 =  {pivot2[[1]]  -  
            pivot2[[1, 3]]/pivot2[[3, 3]]  pivot2[[3]], 
          pivot2[[2]]  - 
            pivot2[[2, 3]]/pivot2[[3, 3]]  pivot2[[3]], 
          pivot2[[3]]};
MatrixForm[pivot3]

Now subtract the appropriate multiple of the second from the first row:

pivot4 =  {pivot3[[1]]  -  
            pivot3[[1, 2]]/pivot3[[2, 2]]  pivot3[[2]], 
          pivot3[[2]], pivot3[[3]]};
MatrixForm[pivot4]

Clean it up by dividing each row by the diagonal terms:

reducedrowecheon =  {pivot4[[1]]/pivot4[[1, 1]], 
          pivot4[[2]]/pivot4[[2, 2]], 
          pivot4[[3]]/pivot4[[3, 3]]};
MatrixForm[reducedrowecheon]

Compare:

MatrixForm[RowReduce[augmented]]

On the money.

T.7.a.iii) A 6D example powered by Mathematica

Here's a linear system:

A =(({{3.4, 0.0, 0.0, -3.3, 3.9, -3.5}, {-1.2, 3.0, 1.8, -1.3, -3.9, 2.6}, {0.9, 3.7, 2.4, -1.8, -1.8, 0.7}, {-1.8, -2.4, -2.6, -1.7, -3.0, -1.0}, {4.1, -1.6, 0.0, 2.0, 0.0, -3.7}, {2.6, 3.8, -2.5, 2.7, -3.1, 0.9}}));
dim = 6;
Clear[x, j];
X = 
        Table[x[j], {j, 1, dim}];
Y =  {1.9, 0, 0, -2.6, 1.8, 
          0.8};
linearsystem = 
      ColumnForm[Thread[A . X == Y]]

Use the reduced row echelon form of the augmented matrix to solve this system.

The linear system is:

linearsystem

The augmented matrix is

augmentedA  = AppendRows[A, Transpose[{Y}]];
MatrixForm[augmentedA]

And look at this:

rowechelon  =  RowReduce[augmentedA];
MatrixForm[rowechelon]

Pick out the seventh column and put

X=Table[rowechelon[[j,7]],{j,1,6}]

Check:

A . X

Y

This is getting to be old hat.

T.7.a.iv) Failure in 4D

Here's a linear system:

A = (({{1.2, -0.8, 0.9, 2.3}, {0.5, -2.4, 0.6, 3.1}, {2.2, -5.6, 2.1, 8.5}, {-3.1, 0, -2.1, -3.8}}));
\(\(Y = {\(-1.0\), 0.4, \(-1.4\), 0.7};\)\(
\)\)
Clear[x, k];
X = Table[x[k], {k, 1, 4}];

linearsystem = 
      ColumnForm[Thread[A . X == Y]]

Use the reduced row echelon form of the augmented matrix to try to solve this system.

The linear system is:

linearsystem

The augmented matrix is

augmentedA  = AppendRows[A, Transpose[{Y}]];
MatrixForm[augmentedA]

And look at this:

rowechelon  =  RowReduce[augmentedA];
MatrixForm[rowechelon]

Garbage. 
Reason: A is not invertible.

Inverse[A]

T.7.b.i) When you want to use Gaussian elimination (reduced row echelon form)

If you're everything by hand calculations, you really have no choice but to use Gaussian elimination.
And that's why most printed matrix texts emphasize Gaussian elimination so heavily.

But when you work with larger matrices with decimal entries, hand calculations can get out of control very,very  quickly.  That's why most matrix printed texts usually limit themselves to artificial coefficient matrices with integer entries. 

On the other hand, when you deal with huge matrices, calculation of SVD and inverses can get out of hand. In such situations, you have little choice but to go with Gaussian elimination.

But when you do, use professionally written software.

Naively written routines often break down when an appropriate multiple of a row is very small.

If Mathematica bails out on your large system, then MatLab or any other product using Linpack make good alternative choices.

T.7.b.ii) Pseudoinverse

Is Gaussian elimination (reduced row echelon form)useful for getting at the pseudo inverse?

No.