MGM.06 Beyond 3D 
 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].";

Needs["Statistics`ContinuousDistributions`"];

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

B.1) Hitting and hanging in higher dimensions

B.1.a.i) Hitting on 3D and hanging in 2D

Here's a random matrix which hits on 3D points and hangs them in 2D:

hitdim = 3;
hangdim = 2;

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 are points scattered on the surface of the unit sphere in 3D:

hitdim = 3;
\(\(points = 
        Table[Random[NormalDistribution[0, 1]], {k, 1, 
            200}, {hitdim}];\)
\)
\(\(spherepoints = 
        Table[points[[k]]\/\@\(points[[k]] . points[[k]]\), {k, 1, Length[points]}];\)
\)
Clear[pointcolor, k];
\(\(pointcolor[k_] := 
        RGBColor[
          0.5\ \((Sin[
                  N[2\ π]\ spherepoints[[k, 
                      1]]] + 1)\), 
          0.5\ \((Sin[
                  N[π]\ spherepoints[[k, 
                      2]]] + 1)\), 
          0.5\ \((Sin[
                  N[π]\ spherepoints[[k, 
                      3]]] + 1)\)];\)
\)
spherepointplot = 
        Show[Table[
            Graphics3D[{PointSize[0.03], pointcolor[k], 
                Point[spherepoints[[
                    k]]]}], {k, 1, Length[spherepoints]}], 
          Axes -> True, Boxed -> False, 
          AxesLabel -> {"x", "y", "z"}, 
          ViewPoint -> CMView];

Here's what happens to these points after the hit with A:

after = 
        Show[Table[
            Graphics[{PointSize[0.03], pointcolor[k], 
                Point[A . 
                    spherepoints[[
                      k]]]}], {k, 1, 
              Length[spherepoints]}], Axes -> True, 
          AxesLabel -> {"x", "y"}, 
          PlotLabel -> "After"];

Identify and plot the ellipse you see.

Look at the stretch factors:

stretches=SingularValues[A][[2]]

This tells you that a hit with A hangs the results on hangerframe[1] and hangerframe[2].

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

frameup = 
        Show[after, 
          Arrow[stretch[1] hangerframe[1], Tail -> {0, 0}, 
            VectorColor -> Black], 
          Arrow[stretch[2] hangerframe[2], Tail -> {0, 0}, 
            VectorColor -> Black]];

Seeing this, you can easily plot the ellipse you want:

Clear[ellipseplotter, t];
ellipseplotter[t_]=Cos[t] stretch[1] hangerframe[1]+Sin[t] stretch[2] hangerframe[2];

ellipseplot = 
          ParametricPlot[ellipseplotter[t], {t, 0, 2 π}, 
            PlotStyle -> {{Thickness[0.01]}}, 
            DisplayFunction -> Identity];

Show[frameup, ellipseplot];

There you go.

B.1.a.ii) Hitting on 4D and hanging in 2D

Here's a random matrix which hits on 4D  points and hangs them in 2D:

hitdim = 4;
hangdim = 2;

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 is a table of points scattered on the surface of the unit sphere in 4D shown with the 4D coordinates of the first two:

\(\(points = 
        Table[Random[NormalDistribution[0, 1]], {k, 1, 
            200}, {hitdim}];\)
\)
\(\(spherepoints = 
        Table[points[[k]]\/\@\(points[[k]] . points[[k]]\), {k, 1, Length[points]}];\)
\)
spherepoints[[1]]
spherepoints[[2]]

Because these points are in 4D, they cannot be plotted. 
But because A hangs in 2D, you can plot what happens when you hit these points with A:

Clear[pointcolor, k];
\(\(pointcolor[k_] := 
        RGBColor[
          0.5\ \((Sin[
                  N[2\ π]\ spherepoints[[k, 
                      1]]] + 1)\), 
          0.5\ \((Sin[
                  N[π]\ spherepoints[[k, 
                      2]]] + 1)\), 
          0.5\ \((Sin[
                  N[π]\ spherepoints[[k, 
                      3]]] + 1)\)];\)
\)
after = 
        Show[Table[
            Graphics[{PointSize[0.025], pointcolor[k], 
                Point[A . 
                    spherepoints[[
                      k]]]}], {k, 1, 
              Length[spherepoints]}], Axes -> True, 
          AxesLabel -> {"x", "y"}, 
          PlotLabel -> "After"];

Identify and plot the ellipse you see.

Look at the stretch factors:

stretches=SingularValues[A][[2]]

This tells you that a hit with A hangs the results on hangerframe[1] and hangerframe[2].

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

frameup = 
        Show[after, 
          Arrow[stretch[1] hangerframe[1], Tail -> {0, 0}, 
            VectorColor -> Black], 
          Arrow[stretch[2] hangerframe[2], Tail -> {0, 0}, 
            VectorColor -> Black]];

Seeing this, you can easily plot the ellipse you want:

Clear[ellipseplotter, t];
ellipseplotter[t_]=Cos[t] stretch[1] hangerframe[1]+Sin[t] stretch[2] hangerframe[2];

ellipseplot = 
          ParametricPlot[ellipseplotter[t], {t, 0, 2 π}, 
            PlotStyle -> {{Thickness[0.01]}}, 
            DisplayFunction -> Identity];

Show[frameup, ellipseplot];

No problem-o.

B.1.b.i) Hitting on 5D and hanging in 3D

Here's a matrix A which hits on 5D points and hangs them in 3D:

A = {{1.257, -1.807, 1.071, 0.985, 
            1.049}, {-0.380, -1.946, -1.143, -1.106, 
            1.562}, {1.535, 1.813, 
            1.266, -0.280, -0.695}};

MatrixForm[A]

Here is a table of points scattered on the surface of the unit sphere in 5D shown with the 5D coordinates of the first two:

hitdim = 5;
hangdim = 3;
points = 
          Table[Random[NormalDistribution[0, 1]], {k, 1, 300}, {hitdim}];

spherepoints = 
          Table[points[[k]]/(Sqrt[points[[k]] . points[[k]]]), {k, 1, Length[points]}];

spherepoints[[1]]
spherepoints[[2]]

Because these points are in 5D, they can't be plotted, but because A hangs in 3D, you can plot what happens when you hit these points with A:

Clear[pointcolor, k];
\(\(pointcolor[k_] := 
        RGBColor[
          0.5\ \((Sin[
                  N[2\ π]\ spherepoints[[k, 
                      1]]] + 1)\), 
          0.5\ \((Sin[
                  N[π]\ spherepoints[[k, 
                      2]]] + 1)\), 
          0.5\ \((Sin[
                  N[π]\ spherepoints[[k, 
                      3]]] + 1)\)];\)
\)
\(\(stretches = \(SingularValues[
            A]\)[[2]];\)
\)
ranger = Max[stretches];
after = 
        Show[Table[
            Graphics3D[{PointSize[0.025], pointcolor[k], 
                Point[A . 
                    spherepoints[[
                      k]]]}], {k, 1, 
              Length[spherepoints]}], Axes3D[ranger], Axes -> True, 
          AxesLabel -> {"x", "y", "z"}, Boxed -> False, 
          ViewPoint -> CMView, 
          PlotRange -> {{-ranger, ranger}, {-ranger, 
                ranger}, {-ranger, ranger}}, 
          PlotLabel -> "After"];

Identify and plot the ellipsoid you see.

Look at the stretch factors:

stretches=SingularValues[A][[2]]

This tells you that a hit with A hangs the results on hangerframe[1], hangerframe[2] and hangerframe[3].

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

frameup = 
        Show[after, 
          Table[Arrow[stretch[k] hangerframe[k], Tail -> {0, 0, 0}, 
              VectorColor -> Black, ShaftWidth -> 0.01], {k, 1, 
              3}]];

Seeing this, you can easily plot the ellipsoid:

Clear[ellipsiodplotter, t];
ellipsoidplotter[s_, t_]=Sin[s] Cos[t] stretch[1] hangerframe[1]+Sin[s] Sin[t] stretch[2] hangerframe[2]+Cos[s] stretch[3] hangerframe[3];

{slow, shigh} = {0, π};
{tlow, thigh} = {0, 2 π};

ellipsoidplot=ParametricPlot3D[ellipsoidplotter[s, t],{s, slow, shigh},{t, tlow, thigh},Axes -> True,AxesLabel -> {"x", "y", "z"},Boxed -> False,ViewPoint -> CMView,PlotLabel -> "After",PlotRange -> {{-ranger, ranger}, {-ranger, 
                    ranger}, {-ranger, ranger}},DisplayFunction -> Identity];

Show[ellipsoidplot, 
        Axes3D[ranger], DisplayFunction -> $DisplayFunction];

Rock on.

B.1.b.ii) Shadows on planes determined by the hanger frames

Stay with the same matrix and the same points as in part i) immediately above and look at these plots:

Clear[hang, m, n];
\(\(hang[m_, n_] := \(hang[m, n] = 
          Show[Table[
              Graphics[{PointSize[0.03], pointcolor[k], 
                  Point[{A . 
                        spherepoints[[
                          k]] . hangerframe[m], 
                      A . spherepoints[[
                          k]] . hangerframe[n]}]}], {k, 1, 
                Length[spherepoints]}], Axes -> True, 
            PlotLabel -> {m, n}, 
            PlotRange -> {{\(-ranger\), ranger}, {\(-ranger\), ranger}}, 
            Epilog -> {Thickness[0.01], 
                Circle[{0, 0}, {stretch[m], stretch[n]}]}, 
            DisplayFunction -> Identity]\);\)
\)
Show[hang[1, 2], 
        DisplayFunction -> $DisplayFunction];

Show[hang[1, 3], 
        DisplayFunction -> $DisplayFunction];

Show[hang[2, 3], 
        DisplayFunction -> $DisplayFunction];

What are these plots and how are they related to the ellipsoid and the points plotted in part i)?

▪ The component of A.spherepoints[[k]] in the direction of hangerframe[m] is   
              (A.spherepoints[[k]].hangerframe[m] ) hangerframe[m].
▪ The component of A.spherepoints[[k]] in the direction of hangerframe[n] is   
              (A.spherepoints[[k]].hangerframe[n] ) hangerframe[n].

So the projection of  A.spherepoints[[k]] onto the plane determined by hangerframe[m] and hangerframe[n] is
  (A.spherepoints[[k]]. hangerframe[m]) hangerframe[m] + (A.spherepoints[[k]]. hangerframe[m]) hangerframe[m].

The plots above use hangerframe[m] and hangerframe[n] as coordinate axes and plot the perpendicular frame coordinates
          {(A.spherepoints[[k]]. hangerframe[m]), ((A.spherepoints[[k]]. hangerframe[n]))}

To see what they mean, look at this:

Clear[section];
\(\(section[m_, n_] := \(section[m, n] = 
          ParametricPlot3D[
            r\ \((stretch[m]\ Cos[t]\ hangerframe[m] + 
                  stretch[n]\ Sin[t]\ hangerframe[n])\), {r, 0, 1}, {t, 0, 
              2\ π}, PlotPoints -> {2, Automatic}, 
            PlotLabel -> {m, n}, Axes -> True, 
            AxesLabel -> {"x", "y", "z"}, 
            Boxed -> False, ViewPoint -> CMView, 
            PlotRange -> All, 
            DisplayFunction -> Identity]\);\)
\)
Show[section[1, 2], frameup, 
        DisplayFunction -> $DisplayFunction];

hang[1,2] plots the shadow (projection) of these points on the ellipse you see above:

Show[hang[1, 2], 
        DisplayFunction -> $DisplayFunction];

To see what hang[1,3] plots, look at this:

Show[section[1, 3], frameup, 
        DisplayFunction -> $DisplayFunction];

hang[1,3] plots the shadow (projection) of these points on the ellipse you see above:

Show[hang[1, 3], 
        DisplayFunction -> $DisplayFunction];

To see what hang[2,3] plots, look at this:

Show[section[2, 3], frameup, 
        DisplayFunction -> $DisplayFunction];

hang[2,3] plots the shadow (projection) of these points on the ellipse you see above:

Show[hang[2, 3], 
        DisplayFunction -> $DisplayFunction];

Here they all are:

Show[
        GraphicsArray[{hang[1, 2], hang[1, 3], hang[2, 3]}]];

The three main cross sections of the points and the ellipsoid plotted above.

B.1.c.i) Hitting on 5D and hanging in 4D

Here's a random matrix which hits on 5D points and hangs them in 4D:

hitdim = 5;
hangdim = 4;

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 is a table of points scattered on the surface of the unit sphere in 5D shown with the 5D coordinates of the first two:

\(\(points = 
        Table[Random[NormalDistribution[0, 1]], {k, 1, 
            200}, {hitdim}];\)
\)
\(\(spherepoints = 
        Table[points[[k]]\/\@\(points[[k]] . points[[k]]\), {k, 1, Length[points]}];\)
\)
spherepoints[[1]]
spherepoints[[2]]

This matrix hits on 5D points and hangs them in 4D.
A sample:

A . spherepoints[[1]]

Because these hit points are in 4D, there is no hope in plotting these points, but you can be sure that they are all packed in an ellipsoid in 4D.
How can you try to look at them?

Go to the well again:

stretches=SingularValues[A][[2]]

This tells you that a hit with A hangs points on this 4D perpendicular frame:

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

Table[
      hangerframe[k], {k, 1, 4}]

You cannot plot this frame, but you can do the next best thing. 

You can plot the shadow (projection) of the hit points on the two dimensional slicer of 4D defined by hangerframe[m] and hangerframe[n]. 

Here you go:

\(\(ranger = Max[stretches];\)
\)
Clear[pointcolor, k];
\(\(pointcolor[k_] := 
        RGBColor[
          0.5\ \((Sin[
                  N[2\ π]\ spherepoints[[k, 
                      1]]] + 1)\), 
          0.5\ \((Sin[
                  N[π]\ spherepoints[[k, 
                      2]]] + 1)\), 
          0.5\ \((Sin[
                  N[π]\ spherepoints[[k, 
                      3]]] + 1)\)];\)
\)
Clear[hang, m, n];
\(\(hang[m_, n_] := \(hang[m, n] = 
          Show[Table[
              Graphics[{PointSize[0.03], pointcolor[k], 
                  Point[{A . 
                        spherepoints[[
                          k]] . hangerframe[m], 
                      A . spherepoints[[
                          k]] . hangerframe[n]}]}], {k, 1, 
                Length[spherepoints]}], Axes -> True, 
            PlotLabel -> {m, n}, 
            PlotRange -> {{\(-ranger\), ranger}, {\(-ranger\), ranger}}, 
            Epilog -> {Thickness[0.01], 
                Circle[{0, 0}, {stretch[m], stretch[n]}]}, 
            DisplayFunction -> Identity]\);\)
\)
Show[hang[1, 2], 
        DisplayFunction -> $DisplayFunction];

Think of it. 
Few have ever been able to peer into 4D with the hit of a key the way you are doing right now. 
To see all the 2D cross sections of the 4D ellipsoid you get by hitting the 5D unit sphere with A, just activate:

Show[
        GraphicsArray[{{hang[1, 2], hang[1, 3], hang[1, 4]}, {hang[2, 3], 
              hang[2, 4], hang[3, 4]}}]];

B.1.c.ii) Another high D sample

Does the same thing work for all higher dimensions?

Sure does.
Try another:

hitdim = 5;
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 is a table of points scattered on the surface of the unit sphere in 5D shown with the 5D coordinates of the first two:

\(\(points = 
        Table[Random[NormalDistribution[0, 1]], {k, 1, 
            200}, {hitdim}];\)
\)
\(\(spherepoints = 
        Table[points[[k]]\/\@\(points[[k]] . points[[k]]\), {k, 1, Length[points]}];\)
\)
spherepoints[[1]]
spherepoints[[2]]

This matrix hits on 5D points and hangs them in 7D.
A sample:

A . spherepoints[[1]]

Because these hit points are in 7D, there is no hope in plotting these points, but you can be sure that they are all packed in an ellipsoid in 7D.
To get a look at them, go to the well again:

stretches = SingularValues[A][[2]]

There are five nonzero stretch factors here.

This tells you that A hangs its hits on this five dimensional perpendicular frame within 7D.

Clear[alignerframe, stretch, hangerframe, k];
alignerframe[
          k_] := Eigenvectors[Transpose[A] . A][[
          k]];
stretch[
          k_] := (Sqrt[A . alignerframe[k] . A . alignerframe[k]]);
\(\(hangerframe[k_] := A . alignerframe[k]\/stretch[k];\)
\)
rank = 5;
Table[hangerframe[k], {k, 1, rank}]

You cannot plot this frame, but you can do the next best thing. 

You can plot the shadow (projection) of the hit points on the two dimensional slicer of 7D defined by hangerframe[m] and hangerframe[n]. 
Here you go:

\(\(ranger = Max[Table[stretch[k], {k, 1, hitdim}]];\)
\)
Clear[pointcolor, k];
\(\(pointcolor[k_] := 
        RGBColor[
          0.5\ \((Sin[
                  N[2\ π]\ spherepoints[[k, 
                      1]]] + 1)\), 
          0.5\ \((Sin[
                  N[π]\ spherepoints[[k, 
                      2]]] + 1)\), 
          0.5\ \((Sin[
                  N[π]\ spherepoints[[k, 
                      3]]] + 1)\)];\)
\)
Clear[hang, m, n];
\(\(hang[m_, n_] := \(hang[m, n] = 
          Show[Table[
              Graphics[{PointSize[0.03], pointcolor[k], 
                  Point[{A . 
                        spherepoints[[
                          k]] . hangerframe[m], 
                      A . spherepoints[[
                          k]] . hangerframe[n]}]}], {k, 1, 
                Length[spherepoints]}], Axes -> True, 
            PlotLabel -> {m, n}, 
            PlotRange -> {{\(-ranger\), ranger}, {\(-ranger\), ranger}}, 
            Epilog -> {Thickness[0.01], 
                Circle[{0, 0}, {stretch[m], stretch[n]}]}, 
            DisplayFunction -> Identity]\);\)
\)
Show[hang[1, 2], 
        DisplayFunction -> $DisplayFunction];

To see all of these 5D ellipses within 7D you get by hitting the 5D unit sphere with A, just activate:

Show[
        GraphicsArray[{{hang[1, 2], hang[1, 3], hang[1, 4], 
              hang[1, 5]}, {hang[2, 3], hang[2, 4], hang[2, 5]}, {hang[3, 4], 
              hang[3, 5], hang[4, 5]}}]];

Life in higher dimensions.

B.1.c.iii) Subspaces

Wait a minute. The last matrix hit 5D into 7D but its hangerframe contained only 5 mutually perpendicular vectors:

rank = 5;
Table[hangerframe[k], {k, 1, rank}]

What happened to the other two dimensions?

The hit by A on 5D had no chance of filling up all of 7D. When you hit A on all of 5D, you get what folks like to call a five dimensional subspace of 7D.
This subspace consists of all the 7D vectors Y of the form

Y = (∑_(k = 1)^rank)(Y.hangerframe[k]) hangerframe[k]

More on this later.

B.2) Rank of a matrix. The advantage you get from full rank.

B.2.a.i) Rank  is the number of non-zero stretch factors.
What this means

Here's a matrix A which hits on 6D points and hangs them in 5D:

\(\(A = {{1.257, \(-1.807\), 1.071, 0.985, 
            1.049, \(-3.247\)}, {\(-0.380\), \(-1.946\), \(-1.143\), \(-1.106\), 1.562, 2.150}, 
\ \ \ \ \ \ \ {1.535, 1.813, 
            1.266, \(-0.280\), \(-0.695\), \(-0.945\)}, {\(-1.915\), \(-3.759\), \(-2.409\), \(-0.826\), 2.257, 3.095}, 
\ \ \ \ \ \ \ \ {3.172, 1.952, 
            3.480, 1.811, \(-1.208\), \(-6.342\)}};\)\)
MatrixForm[A]

Here are the SVD stretch factors for A:

stretches = SingularValues[A][[2]]

Three positive stretch factors.  
This tells you that the rank of A is:

rank  =  Length[stretches]

What does this mean?

Exactly three non-zero stretch factors signals that with A are framed up by a perpendicular frame consisting of 3 5D vectors .
They are:

Clear[hangerframe];
{hangerframe[1], hangerframe[2], hangerframe[3]} = SingularValues[A][[1]]

Watch it happen for a random hit with A:

Clear[x, i];
\(	\(x[i_]\  := \ Random[Real, {\(-10\), 10}];\)
		\)
randomhit =  A . {x[1], x[2], x[3], x[4], x[5], x[6]}

To see whether this randomhit is framed by the three hangerframe vectors check whether
              randomhit=(∑_(k = 1)^rank)(randomhit.hangerframe[k]) hangerframe[k]

Clear[k];
(∑_(k = 1)^rank)(randomhit.hangerframe[k]) hangerframe[k]

They agree. 
And it will happen this way every time because this is what rank means.

B.2.a.ii) Another example

Here's a matrix A which hits on 5D points and hangs them in 7D:

A = {{1.257, -1.807, 1.071, 0.985, 
            1.049}, 
{-0.380, -1.946, -1.143, -1.106, 1.562}, 
{1.535, 1.813, 
            1.266, -0.280, -0.695}, 
{-1.915, -3.759, -2.409, -0.826, 2.257}, 
{3.172, 1.952, 
            3.480, 1.811, -1.208}, 
{2.345, 4.169, 
            3.052, -4.064, 1.094}, 
{0.827, -2.217, 
            0.428, 5.875, -2.302}};
MatrixForm[A]

Here are the SVD stretch factors for A:

stretches = SingularValues[A][[2]]

Four positive stretch factors.  The rank of A is:

rank  =  Length[stretches]

What does this mean?

Exactly four non-zero stretch factors signals that hits with A are framed up by a perpendicular frame consisting of 4 7D vectors .
They are:

Clear[hangerframe];
{hangerframe[1], hangerframe[2], hangerframe[3], 
        hangerframe[4]} = SingularValues[A][[1]]

Watch it happen for a random hit with A:

Clear[x, i];
\(	\(x[i_]\  := \ Random[Real, {\(-10\), 10}];\)
		\)
randomhit =  A . {x[1], x[2], x[3], x[4], x[5]}

To see whether this randomhit is framed by the three hangerframe vectors check whether
              randomhit=(∑_(k = 1)^rank)(randomhit.hangerframe[k]) hangerframe[k]

Clear[k];
(∑_(k = 1)^rank)(randomhit.hangerframe[k]) hangerframe[k]

They agree. 
And it will happen this way every time because that's what rank means.

B.2.b.i) Full rank

Given a matrix A that hits on hitdimD and hangs in hangdimD, what does it mean to say that A is of full rank?
Why do most folks regard full rank as a big bonus?

Saying that A is of full rank is the same as saying the rank of A is the same as hitdim. 
This is the same as saying that none of the SVD stretch factors for A is 0.
This is the same as saying hits with A squash NO non-zero vectors onto  {0,0,0,...,0}.

Full rank is a big bonus because if A is of full rank, then you know that corresponding linear systems A.X = Y cannot be under determined (have more than one solution). More on this below.

B.3) Linear algebra:  
Given a linear system A.X = Y,  put 
            Ytest = (∑_(k = 1)^rank)(Y . 
                  hangerframe[k]) hangerframe[k]
If Ytest ≠ Y, then the linear system A.X = Y has no solution X.

If Ytest = Y, then the linear system A.X = Y has a solution X given by
          Xsol=(∑_(k = 1)^rank)(((Y . hangerframe[k])/(stretch[k]))) alignerframe[k].
If  A is of full rank, then Xsol is the only solution of A . X = Y (exactly determined).
If  A is not of full rank, then Xsol is one of many solutions of A . X = Y (under determined)

B.3.a.i) Given a linear system A.X = Y,  put 
            Ytest=(∑_(k = 1)^rank)(Y.hangerframe[k]) hangerframe[k].
If Ytest ≠ Y, then the linear system A.X = Y has no solution for X.
If Ytest = Y, then the linear system A.X = Y has a solution for X

Given a linear system A.X = Y,  do an SVD analysis of the coefficient matrix A and put 
            Ytest=(∑_(k = 1)^rank)(Y.hangerframe[k]) hangerframe[k]

Explain these two facts:
▪Saying that the linear system A . X = Y has at least one solution for X is the same as saying that  Ytest = Y.

▪ Saying that the linear system A . X = Y has no solution for X is the same as saying that Ytest ≠ Y

In the problem on rank above, you saw that if Y is a hit with A, then

Y =(∑_(k = 1)^rank)(Y.hangerframe[k]) hangerframe[k]= Ytest.

So if Y ≠ Ytest, there is no hope of solving A.X = Y for X. 

If Y = Ytest, then you will see in part ii) below why

(∑_(k = 1)^rank)(((Y . hangerframe[k])/(stretch[k]))) alignerframe[k]

is a guaranteed solution of A.X = Y.

B.3.a.ii) Given a linear system A.X = Y,  put     
          Ytest=(∑_(k = 1)^rank)(Y.hangerframe[k]) hangerframe[k].
 If Ytest = Y, then a guaranteed solution is
           Xsol=(∑_(k = 1)^rank)(((Y . hangerframe[k])/(stretch[k]))) alignerframe[k]

Explain this:
Given a linear system A.X = Y,  put

Ytest=(∑_(k = 1)^rank)(Y.hangerframe[k]) hangerframe[k].

If Ytest = Y, then a guaranteed solution is

Xsol=(∑_(k = 1)^rank)(((Y . hangerframe[k])/(stretch[k]))) alignerframe[k].

In this setup, you are given that given that

Y = Ytest  =(∑_(k = 1)^rank)(Y.hangerframe[k]) hangerframe[k].

To do the explanation (proof), just  remember that
                A.alignerframe[k]=stretch[k] hangerframe[k], 
so that
                hangerframe[k]=((A.alignerframe[k])/(stretch[k])).
This gives

Y=(∑_(k = 1)^rank)(Y.hangerframe[k]) ((A . alignerframe[k])/(stretch[k])).

Thanks to linearity of matrix hits, this is the same as

Y= A.(∑_(k = 1)^rank)(((Y.hangerframe[k])/(stretch[k]))) alignerframe[k].

This signals that

Xsol = (∑_(k = 1)^rank)(((Y.hangerframe[k])/(stretch[k]))) alignerframe[k]

is a guaranteed solution of the linear system A.X = Y.

B.3.a.iii) If Ytest = Y, and A is of full rank, then Xsol is the only solution of A.X = Y (exactly determined).
 If Ytest = Y, and A is of not of full rank, then Xsol is not the only solution of A.X = Y (under determined).

Given a linear system A.X = Y with Ytest = Y; so that A.X = Y has at least one solution.
Explain these facts:
▪ If  A is of full rank, then Xsol is the only solution of A . X = Y (exactly determined).
▪ If  A is not of full rank, then Xsol is not the only solution of A . X = Y (under determined).

▪ Given that A is of full rank,  you already know that A.Xsol = Y.
The goal  is to explain why Xsol is the only solution.
To this end, suppose you've got your hands on another solution Xother. This gives you
           A . Xother = Y
and     
            A . Xsol = Y.
Subtract to get
           A . Xother - 
          A . Xsol = Y - Y = {0, 0, . . ., 
            0}.
This gives you
          A . (Xother - Xaol) = {0, 0, . . ., 
          0}.
Because A is of full rank, hits with A squash NO non-zero vectors onto  {0,0,...,0}. 
This tells you that
           Xother - Xsol = {0, 0,  ... , 0}.
And this tells you that
          Xother = Xsol.
The upshot:

Xsol=(∑_(k = 1)^rank)(((Y . hangerframe[k])/(stretch[k]))) alignerframe[k].

▪ Given that A is of not of full rank,  you already know that A.Xsol = Y.
The goal  is to explain why Xsol is the not the only solution.

Because A is not of full rank, you are guaranteed that there is an Xzero with 
           Xzero ≠ {0, 0, . . ., 0}  
 that is squashed onto {0, 0, . . ., 0}  with a hit with A.

In other words,
           A . Xzero = {0, 0, . . ., 0}.
This also gives
           A . (t Xzero) =  t A . Xzero  =  {0, 
            0, . . ., 0}
for any number t.
Add this to
           A . Xsol = Y 
to get
           A . Xsol +  
          A . (t Xzero) =  Y  + {0, 
              0, . . ., 0}  =  Y.
for any number t.
This gives
            A . (Xsol +  t Xzero) =  A . Xsol  +  
            A . (t Xzero) =   Y
 for any number t.
 As t varies from -infinity to infinity, this gives you infinitely many solutions of 
           A . X = Y.

B.3.b.i) Sample case: Exactly one solution

Here's a matrix which hits on 3D and hangs in 5D:

\(\(A\  = \ {{0.722, \(-2.855\), 
            2.841}, 
\ \ \ \ \ \ \ \ \ \ \ \ {1.315, \(-0.229\), \(-1.267\)}, 
\ \ \ \ \ \ \ \ \ \ \ \ {2.884, \(-1.006\), 
            0.789}, 
\ \ \ \ \ \ \ \ \ \ \ \ {0.192, 
            2.081, \(-2.131\)}, 
\ \ \ \ \ \ \ \ \ \ \ \ {1.929, \(-1.171\), 
            2.328}};\)
\)
MatrixForm[A]

This matrix hits on 3D and hangs its hits in 5D.

Clear[x, y, z];
A . {x, y, z}

Use this matrix to make a system of linear equations
          A . X = Y:

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

Given a specific 
         Y = {y[1], y[2], y[3], y[4], y[5]}
saying that the corresponding linear system has a solution is the same as saying that you can find specific numbers x[1], x[2] and x[3] so that when you put 
          X = {x[1], x[2], x[3]}, 
then you get 
           A . X = Y.

Stay with the same linear system as above but go with this specific    
       Y = {y[1], y[2], y[3], y[4], y[5]}:

Y = {24.1787, -6.56968, 11.4069, -16.8513, 
          20.8961};
hitdim = 3; hangdim = 5;
Clear[x, k];
\(\(X = Table[x[k], {k, 1, hitdim}];\)
\)
linearsystem = ColumnForm[Thread[A . X == Y]]

If this system has a solution, then come up with one. 
Say whether it has additional solutions.

Do your SVD analysis of the stretch factors of the coefficient matrix A:

stretches = SingularValues[A][[2]]

This determines the rank of A:

rank = Length[stretches]

Dial up the hangerframe vectors corresponding to the non-zero stretch factors and form Ytest:

Clear[hangerframe];
hangerframe[k_] := SingularValues[A][[1]][[k]];
 Ytest = (∑_(k = 1)^rank)(Y . 
                hangerframe[k]) hangerframe[k]

Compare:

Y

Bingo.
               Y  = Ytest.
This tells you that there is a solution. One guaranteed solution is:

       \(Xsol = ∑\+\(k = 1\)\%rank\((Y . hangerframe[k]\/stretch[k])\)\ alignerframe[k]\).:

Clear[alignerframe, stretch, hangerframe, k];
hangerframe[
        k_] := SingularValues[
            A][[1]][[k]]
stretch[
        k_] := SingularValues[
            A][[2]][[k]]
\(\(alignerframe[
          k_] := \(\(SingularValues[
              A]\)[[3]]\)[[k]];\)
\)
\(\(rank = 
        Length[\(SingularValues[
              A]\)[[2]]];\)
\)
Xsol = 
      Expand[(∑_(k = 1)^rank)(( Y . hangerframe[k] alignerframe[k])/stretch)[k]]

Check:

A . Xsol
Y

Good.

To see whether there are more solutions, compare hitdim and the rank:

hitdim

rank

Because hitdim = rank, the coefficient matrix A is of full rank. 
This tells you the Xsol is the only solution.

B.3.b.ii) Sample case: No solutions

Stay with the same linear system as in part i) but go with this specific    
          Y = {y[1], y[2], y[3], y[4], y[5]}:

Y = {-9.62776, 6.26124, 0.436794, 
          8.02976, -6.23616};
hitdim = 3;
\(\(hangdim = 5;\)
\)
Clear[x, k];
\(\(X = Table[x[k], {k, 1, hitdim}];\)
\)
linearsystem = ColumnForm[Thread[A . X == Y]]

If this system has a solution, then come up with one. 
Say whether it has additional solutions.

Do your SVD analysis of the stretch factors of the coefficient matrix A:

stretches = SingularValues[A][[2]]

This determines the rank of A:

rank = Length[stretches]

Dial up the hangerframe vectors corresponding to the non-zero stretch factors and form Ytest:

Clear[hangerframe];
hangerframe[k_] := SingularValues[A][[1]][[k]];
 Ytest = (∑_(k = 1)^rank)(Y . 
                hangerframe[k]) hangerframe[k]

Compare:

Y

Close, but no cigar.
Y is not the same as Ytest.
The call: This linear system has no solution.

B.3.b.iii) Sample case: More than one solution

Here's a new matrix A:

\(\(A = {{0.6146, 0.5210, \(-0.1295\), 1.9475, 
            1.9705, \(-0.4572\)}, 
\ \ \ \ \ 		{0.8440, \(-0.2634\), \(-0.1994\), \(-0.7530\), \(-0.6166\), 7.0777}, 
\ \ \ \ \ 		{\(-0.8955\), 
            1.2243, \(-0.4330\), \(-1.4189\), 0.2747, 
            4.2659}, 
	\ \ \ \ \ 	{1.2645, \(-2.7612\), 3.9456, 1.4984, 
            0.1693, \(-4.5239\)}};\)
\)
MatrixForm[A]

Here is a 4D 
            Y = {y[1], y[2], y[3], y[4]}
together with the corresponding linear system:

Y = {10.2452, 2.8353, -8.44979, -5.24377};
hitdim = 6;
\(\(hangdim = 4;\)
\)
Clear[x, k];
\(\(X = Table[x[k], {k, 1, hitdim}];\)
\)
linearsystem = ColumnForm[Thread[A . X == Y]]

If this system has a solution, then come up with one. Say whether it has additional solutions.

Do your SVD analysis of the stretch factors of the coefficient matrix A:

stretches = SingularValues[A][[2]]

This determines the rank of A:

rank = Length[stretches]

Dial up the hangerframe vectors corresponding to the non-zero stretch factors and form Ytest:

Clear[hangerframe];
hangerframe[k_] := SingularValues[A][[1]][[k]];
 Ytest = (∑_(k = 1)^rank)(Y . 
                hangerframe[k]) hangerframe[k]

Compare:

Y

Bingo.
               Y  = Ytest.
This tells you that there is a solution. One guaranteed solution is:

       \(Xsol = ∑\+\(k = 1\)\%rank\((Y . hangerframe[k]\/stretch[k])\)\ alignerframe[k]\).:

Clear[alignerframe, stretch, hangerframe, k];
hangerframe[
        k_] := SingularValues[
            A][[1]][[k]]
stretch[
        k_] := SingularValues[
            A][[2]][[k]]
\(\(alignerframe[
          k_] := \(\(SingularValues[
              A]\)[[3]]\)[[k]];\)
\)
\(\(rank = 
        Length[\(SingularValues[
              A]\)[[2]]];\)
\)
Xsol = 
      Expand[(∑_(k = 1)^rank)(( Y . hangerframe[k] alignerframe[k])/stretch)[k]]

Check:

A . Xsol
Y

Good.

To see whether there are more solutions, compare hitdim and the rank:

hitdim

rank

Because hitdim >  rank, the coefficient matrix A is of not of full rank. 
This tells you the Xsol is the not the only solution.

B.3.b.iv) Getting all the solutions

Stay with the same linear system as in part i) immediately above.
Come up with all solutions.

You already know Xsol is a solution. Check:

A . Xsol
Y

This gives one solution namely:

Xsol

To get the others, look at the rank of A and the hit dim of A

rank
hitdim

The fact that rank < hitdim signals that there are more solutions.
Reason: Hits with A use only the four (= rank of A) aligner frames calculated above and zeroes out 6-4 = 2 others. They are:

Clear[zeroedalignerframe];
zeroedalignerframe[k_] := NullSpace[A][[
          k]];
Table[zeroedalignerframe[k], {k, 1, hitdim - rank}]

Check:

Table[A . zeroedalignerframe[k], {k, 1, hitdim - rank}]

All the solutions are parameterized by:

Clear[solution, s, t];
solution[s_, t_] = 
      Xsol + s zeroedalignerframe[1] + t zeroedalignerframe[2]

Try it:

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

Perfecto.

B.4) If the linear system A.X = Y has no solutions at all (overdetermined), then  
         Xapprox=(∑_(k = 1)^rank)(((Y . hangerframe[k])/(stretch[k]))) alignerframe[k]
is the best approximate solution you can get. 
In this case, folks call Xapprox by the name "best least squares solution"

B.4.a.i) Best approximate solution is  Xapprox=(∑_(k = 1)^rank)(((Y . hangerframe[k])/(stretch[k]))) alignerframe[k]

Here's a linear system:

A=(({{-0.91, 0.83}, {-0.62, -0.28}, {-0.30, 1.20}}));

\(\(Y = {\(-1.05\), 0.35, \(-1.41\)};\)\(
\)\)
hitdim = 2;
\(\(hangdim = 3;\)\(
\)\)
Clear[x, k];
X = Table[x[k], {k, 1, hitdim}];

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

In search of a solution X = {x[1], x[2]},you do your SVD analysis of the coefficient matrix A:

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

Dial up the hangerframe vectors corresponding to the non-zero stretch factors and form Ytest:

Clear[hangerframe];
hangerframe[k_] := SingularValues[A][[1]][[k]];
 Ytest = (∑_(k = 1)^rank)(Y . 
                hangerframe[k]) hangerframe[k]

Compare:

Y

Y ≠ Ytest; so there is no exact solution. 
What do you do to come up with the best approximate solution you can get?

Look at Y and Ytest again:

Y
Ytest

Ytest is not all that different from Y and you can find Xapprox with A.Xapprox = Ytest. 
You can get it from the same formula use to get solutions. Just put

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

with full assurance that A.Xpprox = Ytest

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

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

Xapprox=Expand[(∑_(k = 1)^rank)(((Y . hangerframe[k])/(stretch[k]))) alignerframe[k]]

And check:

A . Xapprox
Y

Xapprox is the best approximate solution you can get.
If you want to see why Xapprox is the best go to the next part.

B4.a.ii) Why that worked

Stay with the same setup as in part i) and explain why Xapprox  is the best approximate solution of the linear system A.X = Y.is so close to Y.

All possible hits with A are hung on the frame consisting of:

Table[hangerframe[k], {k, 1, rank}]

This plots out as a plane in 3D.
Here is a plot of that plane in the vicinity of A.Xapprox shown together with a plot of Y:

h = 0.2;
k = 0.2;
\(\(planeplot = 
        ParametricPlot3D[
          A . \((Xapprox + s\ alignerframe[1] + 
                t\ alignerframe[2])\), {s, \(-h\), h}, {t, \(-k\), k}, 
          PlotPoints -> {2, 2}, 
          DisplayFunction -> Identity];\)
\)
Ypointplot = Graphics3D[{PointSize[0.03], Point[Y]}];
\(\(Ylabel = 
        Graphics3D[Text["Y", Y + {h\/3, k\/3, k\/3}]];\)
\)
\(\(AXapproxpointplot = 
        Graphics3D[{Red, PointSize[0.03], Point[A . Xapprox]}];\)
\)
Ylabel = Graphics3D[Text["Y", Y + {0, 0, (k/5)}]];
\(\(AXapproxlabel = 
        Graphics3D[
          Text["A.Xapprox", A . Xapprox - {0, 0, k\/3}]];\)
\)
setup  = 
        Show[planeplot, Ypointplot, Ylabel, AXapproxpointplot, AXapproxlabel, 
          Axes -> True, AxesLabel -> {"x", "y", "z"}, 
          Boxed -> False, ViewPoint -> CMView, 
          PlotRange -> All, 
          DisplayFunction -> $DisplayFunction];

The vector running from A.Xapprox to Y seems to be perpendicular to the plane on which A hangs its hits. Confirm:

(Y - A . Xapprox) . hangerframe[1]

(Y - A . Xapprox) . hangerframe[2]

This tells you that Y - A.Xapprox is perpendicular to both the vectors that frame the plane on which A hangs its hits. 
This is enough to guarantee that Y - A.Xapprox is perpendicular to the plane on which A hangs its hits.

And because A . Xapprox is on the plane on which A hangs its hits and because the perpendicular distance is the shortest distance, this tells you that

In short: Xsol is as close to a solution of:

linearsystem

as you or anyone else will ever get.
That's why lost of folks call X the best approximate solution for this linear system.

B.4.a.iii) Other linear systems without exact solutions

Does the same thing work for other linear systems without solutions?

You bet your sweet mouse it does!
As a matter of fact, it works for the same reasons examined above.
To whit:
When you analyze the coefficient matrix A, you get its rank and then you know that A hangs its hits on the part of hangdimD framed by
           {hangerframe[1],hangerframe[2],,. . .,hangerframe[rank]}.
You put

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

When you hit Xapprox with A you get

A. Xapprox=(∑_(k = 1)^rank)(((Y . hangerframe[k])/(stretch[k]))) A.alignerframe[k]

=(∑_(k = 1)^rank)(((Y . hangerframe[k])/(stretch[k])))  stretch[k] hangerframe[k]

=(∑_(k = 1)^rank)(Y . 
                hangerframe[k])   hangerframe[k]

So 
        Y - A . Xapprox = 
        Y - ((∑_(k = 1)^rank)(Y . hangerframe[k]) hangerframe[
                k]).
  
Consequently for each p with 1 ≤  p ≤  rank, you get 
           (Y - A . Xapprox).hangerframe[p]  
        =    (Y.hangerframe[p])-
	  ((∑_(k = 1)^rank)(Y . 
                    hangerframe[k])( hangerframe[k].)hangerframe[p])
        =    (Y.hangerframe[p])-
	       (Y . hangerframe[p]) (hangerframe[p].hangerframe[p])

=   (Y.hangerframe[p])-(Y . hangerframe[p])

= 0.
The upshot:
           (Y - A . Xapprox).hangerframe[p] =0 
for all the individual vectors, hangerframe[p], in the hangerframe for A. 

And this tells you that Y - A . Xapprox is perpendicular to the part of hangdimD on which A hangs its hits.

And because in all dimensions, the perpendicular distance is the shortest distance, this tells you that A.Xapprox is closer to Y than any other possible hit with A.

B.5)  Inverses for invertible matrices and pseudoinverses for non-invertible matrices
If  A=hanger.stretcher.aligner and A is invertible, then (A^(-1))=(aligner^t)unstretcher.(hanger^t) . 
If  A=hanger.stretcher.aligner, then a hit on Y with PseudoInverse[A]  =(aligner^t)unstretcher.(hanger^t)  gives you Xapprox.

B.5.a.i)  If  A=hanger.stretcher.aligner, then a hit on Y with PseudoInverse[A]  =(aligner^t)unstretcher.(hanger^t)  gives you Xapprox

Explain this:
If  A is any matrix and 
               A = hanger.stretcher. aligner, 
then when you hit Y with

PseudoInverse[A]  =(aligner^t)unstretcher.(hanger^t)

you get

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

In other words, A.[PseudoInverse[A].Y] is as close to Y as any other hit with A.

When you hit (aligner^t)unstretcher.(hanger^t) on Y, you get

(aligner^t)unstretcher.(hanger^t).Y

=  (aligner^t)unstretcher.{Y . hangerframe[1] , 
              Y . hangerframe[2], . . ., 
              Y . hangerframe[rank]}

=  (aligner^t).{ ((Y . hangerframe[1])/(stretch[1])),((Y . hangerframe[2])/(stretch[2])),. . .,((Y . hangerframe[rank])/(stretch[rank]))}

=   (∑_(k = 1)^rank)(((Y . hangerframe[k])/(stretch[k]))) alignerframe[k]

= Xapprox =  PseudoInverse[A].Y

Try it out on a sample random matrix A

hitdim = Random[Integer, {3, 5}];
\(\(hangdim = Random[Integer, {4, 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 is (aligner^t)unstretcher.(hanger^t):

rank=Length[SingularValues[A][[2]]];
Clear[alignerframe, k];
alignerframe[k_]:=SingularValues[A][[3,k]];
aligner = Table[alignerframe[k], {k, 1, rank}];
Clear[stretch];
stretch[k_]:=SingularValues[A][[2,k]];
stretcher = 
            DiagonalMatrix[Table[stretch[k], {k, 1, rank}]];

Clear[hangerframe, k];
hangerframe[k_]:=SingularValues[A][[1,k]];
hanger = Transpose[Table[hangerframe[k], {k, 1, rank}]];

unstretcher=DiagonalMatrix[Table[(1/(stretch[k])),{k, 1, rank}]];
MatrixForm[Transpose[aligner].unstretcher.Transpose[hanger]]

Compare with Mathematica's built in instruction:

MatrixForm[PseudoInverse[A]]

B.5.a.ii)  If  A=hanger.stretcher.aligner and A is invertible, then (A^(-1))= PseudoInverse[A]  =(aligner^t)unstretcher.(hanger^t)

Explain this:
If  A is any invertible matrix and 
               A = hanger.stretcher. aligner, 
then (A^(-1))

(A^(-1)) = PseudoInverse[A]  =(aligner^t)unstretcher.(hanger^t).

Go with a kD invertible matrix A and any vector Y in kD.
The solution of the linear system A.X = Y is
            Xsol =  (A^(-1)).Y .
On the other hand,
            Xappox = PseudoInverse[A].Y 
is as close as it is possible to get to  solution.  
And because Xsol is the exact solution, 
            Xapprox = Xsol.
So 
            PseudoInverse[A].Y  =  (A^(-1)).Y
for all Y's in kD.
This tells you that  PseudoInverse[A] and (A^(-1)) are the same matrix.

Try it out on this sample invertible matrix:

hitdim = 5;
\(\(hangdim = 5;\)\(
\)\)
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 calculation of (A^(-1))

MatrixForm[Inverse[A]]

And here's Mathematica's calculation of PseudoInverse[A]

MatrixForm[PseudoInverse[A]]

B.5.a.iii) Comparing inverse and pseudoinverse

When and how do you use the pseudoinverse instead of the inverse?

You can use the inverse on matrices that are invertible.
You want to use the pseudoinverse on matrices that are not invertible.

Here's a sample of what you can do with the pseudo inverse:

A =  (({{0.45, 0.43, -1.06, -1.65}, {-1.44, -1.28, 0.31, 1.60}, {0.079, -1.96, 0.86, 1.73}, {-0.57, 0.58, 1.82, -1.46}, {1.39, 0.32, 0.44, -0.78}, {1.89, -0.49, -1.12, 1.11}, {-1.21, 0.55, 1.50, -0.64}}));
MatrixForm[
      A]

This matrix is not invertible:

Inverse[A]

Yet there is no problem with the pseudoinverse of A:

MatrixForm[PseudoInverse[A]]

Example: Nailing a solution using the pseudoinverse
You can use PseudoInverse[A] to check whether the linear system A.X = Y has a solution for this 7D Y:

Y  =  {-6.4855, 5.2153, 
          8.63207, -4.7117, -0.9053, 
          5.6479, -2.9046};
Xapprox = PseudoInverse[A].Y

Compare:

A . Xapprox

Y

The call: 
Nailed an solution.

Example: Nailing a reasonably good approximate solution using the pseudoinverse
You can use PseudoInverse[A] to check whether the linear system A.X = Y has a solution for this new 7D Y:

Y  =  {-5.28, 2.44, 3.09, -5.24, -1.13, 
          5.35, -2.71};
Xapprox = PseudoInverse[A].Y

Compare:

A . Xapprox

Y

The call:
Strictly speaking, there is no solution, but the approximate solution Xapprox is not all that bad.

Example: Using the pseudoinverse to see a linear system that is not even close to having a solution
You can use PseudoInverse[A] to check whether the linear system A.X = Y has a solution for this new 7D Y:

Y  =  {-4.28, 0.44, 
          6.09, -3.24, -0.13, -4.35, -0.71};
Xapprox = PseudoInverse[A].Y

Compare:

A . Xapprox

Y

The call: Xapprox is the best that can be done. But the approximation sucks.
For this Y, the linear system A.X = Y is not even close to having a solution. 
It's a sorry excuse for a linear system.

B.5.a.iii)  An invertible matrix is a square matrix of full rank. Other matrices are not invertible

Explain these shots from the hip:
▪ If A is a matrix with hitdim < hangdim, then A cannot be invertible.
▪ If A is a matrix with hitdim > hangdim, then A cannot be invertible.
▪ If A is a matrix with hitdim = hangdim (i.e. A is a square matrix) , then saying that A is invertible is the same as saying that A is of full rank.

▪ If A is a matrix with hitdim < hangdim, then A cannot be invertible.
Explanation: Hits with A can fill up at most hitdim dimensions of hangdimD. Because hitdim < hangdim, hits with A cannot fill all of hangdim D. There are at least hangdim-hitdim dimensions left out.

The upshot:  There are many vectors Y in hangdim D that are not hits with A. You cannot undo A on these vectors. So A is not invertible.

▪ If A is a matrix with hitdim > hangdim, then A cannot be invertible.
Explanation: In this case, hits with A squash dimensions. There must be a vector Xsquash with ||Xsquash|| ≠ 0 in hitdimD that a hit with A squashes to {0,0,...,0}. 
But A.{0,0,...,0} = {0,0,...,0}.
So         
         A.Xsquash = A.{0,0,...,0} = {0,0,...,0}
If A is invertible, then 
             Xsquash =  (A^(-1))A.Xsquash = (A^(-1))A.{0,0,...,0} = {0,0,,...,0}.
This cannot happen because ||Xsquash|| ≠ 0.

▪ If A is a matrix with hitdim = hangdim (i.e. A is a square matrix) , then saying that A is invertible is the same as saying that A is of full rank..
Explanation: If A is a matrix with hitdim = hangdim (i.e. A is a square matrix) , then saying that A is invertible is the same as saying that given a Y in handdimD, the corresponding linear system
        AX = Y
is exactly determined. 
This is the same as saying that A is of full rank.

B.6) Transposes in higher dimensions:
The rank of A is the same as the rank of (A^t).
The nonzero stretch factors for A are the same as the nonzero stretchfactors for (A^t)
The non-zero aligner frame for A is the non-zero hangerframe for (A^t). 
The non-zero hanger frame for A is the non-zero aligner frame for (A^t).

B.6.a.i) A test case

Here's a matrix A:

A=(({{0.43, -1.32, -2.26, 1.47, 4.13, 2.32}, {-0.34, 0.34, 0.00, 7.89, -2.96, 0.89}, {2.58, 1.34, 3.82, 2.84, -1.44, 0.78}, {0.61, 1.35, 1.91, 17.2, -6.64, 2.17}, {1.72, -0.65, -0.35, 2.89, 3.41, 2.71}}));
MatrixForm[A]

Calculate the rank of A:

rankA = 
      Length[SingularValues[
            A][[2]]]

Here's the transpose (A^t) of A:

MatrixForm[Transpose[A]]

Calculate the rank of (A^t):

rankAtrans = 
      Length[SingularValues[
            Transpose[
              A]][[2]]]

The rank of (A^t)is the same of A.
Is this an just an accident?

Uh-oh.

If you want to see why, then go on.

B.6.a.ii) Why rank of A = rank of (A^t).

Explain this:
No matter what matrix A you go with, you are guaranteed that
         rank of A = rank of (A^t).

Remember
        A=hanger.stretcher.aligner.
So
        (A^t)=(aligner^t).(stretcher^t).(hanger^t)
             = (aligner^t).stretcher.(hanger^t) 
 because, as a diagonal matrix, stretcher = (stretcher^t). 
           
This that you that A and (A^t) both use the same perpendicular frames for their hits. 
Except the roles are reversed. 
The perpendicular frame A uses for non-zero aligning is used by (A^t)  for non-zero hanging.
The perpendicular frame A uses for non-zero hanging is used by (A^t)  for non-zero aligning.

And they both have the same non-zero stretch factors.
That's why matter what matrix A you go with, you are guaranteed that
         rank of A = rank of (A^t).