This notebook contains additional material to the Mathematica GuideBooks  (Springer-Verlag, New York, 2004) http://www.MathematicaGuideBooks.org . It uses the stylesheet of the Mathematica GuideBooks, which is embedded.

Numerics Volume, Chapter 1

We start with two function definitions for normalizing a vector and calculating the distance between two points.

𝓃 = #/Sqrt[#.#]&; ℓ = Sqrt[(#2 - #1).(#2 - #1)]&;

The vertices of a hypercube in dim dimensions.

(* the vertices of a dimD unit cube *)
vertices[dim_] := vertices[dim] = 
Flatten[Outer[List, Sequence @@ Table[{-1, 1}/2, {dim}]], dim - 1]

The (2^dum) edges of a hypercube in dim dimensions. We will use them later for visualizations.

(* the edge vertices of a dimD unit cube *)
edges[dim_] := edges[dim] = 
(* all n-tuples with one differing element *)
Map[Position[vertices[dim], #][[1, 1]]&, Union[Sort /@ Flatten[
        Table[Table[{vertices[dim][[j]], 
                     MapAt[-#&, vertices[dim][[j]], k]}, {k, dim}],
                      {j, Length[vertices[dim]]}], 1]], {2}]

And these are the dim (dim-1)(2^(dim-3))  2D faces of a hypercube in dim dimensions. We will also use them for visualizations and for minimization.

edgeToFaces[edge_] := 
(* sweep out a face by moving an edge *)
Module[{posis = Flatten[Position[Equal @@@ Transpose[edge], True]]},
       Table[Join[edge, MapAt[-#&, #, posis[[k]]]& /@ 
                  Reverse[edge]], {k, Length[posis]}]]

(* the face vertices of a dimD unit cube *)
faces[dim_] := faces[dim] = 
(* all n-tuples with two differing elements *)
Map[Position[vertices[dim], #][[1, 1]]&, #[[1, -1]]& /@ 
(* sweep out all faces and keep one of each *)
Split[Sort[{Sort[#], #}& /@ Flatten[edgeToFaces /@ 
               Map[vertices[dim][[#]]&, edges[dim], {-1}], 1]],
               #1[[1]] === #2[[1]]&], {2}]

Next we carry out the straightforward implementation of the objective function representing the cumulative relative deviation of all pairwise vertex distances. The two arguments of the function relativeVertexPositionErrorSum are two vectors in (ℝ^dim) that span a 2D hyperplane.

relativeVertexPositionErrorSum[{d1_?(VectorQ[#, NumberQ]&), d2_}] :=
Module[{dim = Length[d1], 𝒹1, 𝒹2, verticesProjected},
 (* two orthogonal vectors *)
 𝒹1 = 𝓃[d1]; 𝒹2 = 𝓃[d2 - 𝒹1.d2 𝒹1]; 
 (* the projected vertices *)
 verticesProjected = {𝒹1.#, 𝒹2.#}& /@ vertices[dim]; 
 (* sum of squares of edge length errors *)
 Sum[Abs[ℓ[verticesProjected[[i]], verticesProjected[[j]]]/
         ℓ[vertices[dim][[i]], vertices[dim][[j]]] - 1]^2, 
           {i, 2^dim}, {j, i - 1}]]

The function projectedEdgesGraphics shows the projected edges.

projectedEdgesGraphics[{d1_?VectorQ, d2_?VectorQ}, opts___] :=
Module[{dim = Length[d1], 𝒹1, 𝒹2, verticesProjected},
 (* two orthogonal vectors *)
 𝒹1 = 𝓃[d1]; 𝒹2 = 𝓃[d2 - 𝒹1.d2 𝒹1];
 (* the projected vertices *)
 verticesProjected = {𝒹1.#, 𝒹2.#}& /@ vertices[dim]; 
 (* show projected vertices and edges *)
 Graphics[{{Hue[0], Thickness[0.001], 
            Line /@ Map[verticesProjected[[#]]&, edges[dim], {-1}]},
           {GrayLevel[0], PointSize[0.02], 
            Point /@ verticesProjected}},
      opts, PlotRange -> All, AspectRatio -> Automatic]]
      
projectedEdgesGraphics[{min_?NumericQ, vals_}] := 
projectedEdgesGraphics[Partition[Last /@ vals, 
                                 Length[vals]/2], PlotLabel -> min]

Here are the results, the number of function evaluations, and the timings for the four methods "SimulatedAnnealing", "RandomSearch", "DifferentialEvolution", and "NelderMead" shown. Modulo rotations, the four projections are identical. And the minima all agree to many digits. To see the progress in finding the minimum, we define an auxiliary function relativeVertexPositionErrorSumL that keeps a list of the function values obtained.

relativeVertexPositionErrorSumL[
             {d1_?(VectorQ[#, NumberQ]&), d2_}, L_] :=
With[{res = relativeVertexPositionErrorSum[{d1, d2}]},
     L = {L, res}; res]

(* show graphic and return timings and function evaluations *)
(Show[GraphicsArray[#[[2, 2]]& /@ #]]; 
 {#[[2, 1]], #[[2, 3]], #[[1]]}& /@ #)&[
Timing[{#, ℒ[#] = {}; projectedEdgesGraphics[
 (* find optimal projection hyperplane parameters *)
 NMinimize[relativeVertexPositionErrorSumL[
                    {{d1x, d1y, d1z}, {d2x, d2y, d2z}}, 
           (* bag for values *) Unevaluated[ℒ[#]]],
                     {d1x, d1y, d1z, d2x, d2y, d2z},
           PrecisionGoal -> 2, Method -> #]], 
           (* function evaluations *) Length[Flatten[ℒ[#]]]}]& /@ 
           (* the four methods *)
            {"SimulatedAnnealing", "RandomSearch",
             "DifferentialEvolution", "NelderMead"}]

Displaying the calculated the objective function values depending on the iteration number shows that "SimulatedAnnealing" and "NelderMead" find a good minimum after less than 200 search steps. The "RandomSearch" plot shows many locations for a good projection.

convergenceListPlot[l_, opts___] := 
Module[{ℓ = Flatten[l], λ = Last[l]}, 
    Block[{$DisplayFunction = Identity}, 
       (* logarithm of absolute error *)
       ListPlot[DeleteCases[MapIndexed[{#2[[1]], #1}&, 
                 Log[10, Abs[(# - λ)& /@ ℓ]]],
                 {_, Indeterminate}], PlotRange -> All]]];

Show[GraphicsArray[#]]& /@ 
Partition[convergenceListPlot[ℒ[#], PlotLabel -> #]& /@ 
           {"SimulatedAnnealing", "RandomSearch",
            "DifferentialEvolution", "NelderMead"}, 2];

The last result suggests the use of "SimulatedAnnealing" and/or "NelderMead" for higher dimensions in the following. To speed up the calculations, we implement a compiled version of relativeVertexPositionErrorSum. We also allow for projections into 3D, 4D, … hyperspaces.

(* the main compiled function *)
cf[dim_] := cf[dim] = 
Compile[{{ds, _Real, 2}},
Module[{𝒹𝓈 = ds, λ = Length[ds], verticesProjected, 𝒱s = #},
 (* make orthogonal vectors *)
 Do[𝒹𝓈[[k]] = #/Sqrt[#.#]&[ds[[k]] - Sum[𝒹𝓈[[j]] 𝒹𝓈[[j]].ds[[k]], 
    {j, λ - 1}]], {k, λ}];
 (* the projected vertices *)
 verticesProjected = Table[𝒱s[[k]].#& /@ 𝒹𝓈, {k, Length[𝒱s]}];
 (* sum of squares of edge length errors *)
 Sum[Sum[Abs[Sqrt[#.#]&[verticesProjected[[i]] - verticesProjected[[j]]]/
             Sqrt[#.#]&[𝒱s[[i]] - 𝒱s[[j]]] - 1]^2,
        {j, i - 1}], {i, Length[𝒱s]}]]]&[vertices[dim]];
        
(* interface function to the compiled function *)
relativeVertexPositionErrorSumCF[
        {d1_?(VectorQ[#, NumberQ]&), d2__}] := cf[Length[d1]][{d1, d2}];

Here are these two methods applied to hypercubes of dimensions four, five, and six. Again, the two methods agree quite well with one another.

Do[Show[GraphicsArray[#]]&[projectedEdgesGraphics[
     NMinimize[{relativeVertexPositionErrorSumCF[
                        Table[d[j][k], {j, 2}, {k, dim}]]},
               Flatten[Table[d[j][k], {j, 2}, {k, dim}]],
               PrecisionGoal -> 2, Method -> #]]& /@ 
               {"SimulatedAnnealing", "NelderMead"}], {dim, 4, 6}]

For comparison we show projections resulring from determining principal components.

(* make a square k×k array of graphics *)
toSquareArray[l_, k_] := 
Module[{mod = Mod[Length[l], k]}, 
If[mod === 0, Partition[l, k],
   Append[Partition[l, k], Join[Take[l, -mod], 
          Table[Graphics[{}], {k - mod}]]]]]

Function[d, (* print dimension *)
 StylePrint["dim = " <> ToString[d], "PrintText"];
 (* project onto principal components *)
 Module[{svd = SingularValueDecomposition[Transpose@N[vertices[d]]][[1]]},
  Show[GraphicsArray[#]]& /@
  toSquareArray[Flatten[Table[projectedEdgesGraphics[{svd[[i]], svd[[j]]}],
                       {i, 2, d}, {j, i + 1, d}]], 6]]] /@ {6, 7}

Next we project into 3D and show the 2D faces of the resulting 3D objects as holed planes. Here are some function to show the projected 2D faces of the hypercubes.

(* quadrilateral spanned by p1, p2, p3, p4 *)
makeHoled3DPolygon[{p1_, p2_, p3_, p4_}, pp_, f_] := 
Module[{points = Table[p1 + i j/pp^2 (p1 + p3 - p2 - p4) +
                       i/pp (p2 - p1) + j/pp (p4 - p1), 
                {i, 0, pp}, {j, 0, pp}], polys},
 polys = Table[Polygon[{points[[i, j]], points[[i + 1, j]], 
                        points[[i + 1, j + 1]], points[[i, j + 1]]}], 
               {i, pp}, {j, pp}];
 makeHole[#, f]& /@ Flatten[polys]]

makeHole[Polygon[l_], f_] := (* cut a hole in a polygon *)
Module[{mp = Plus @@ l/Length[l], L}, L = (mp + f(# - mp))& /@ l;
       {MapThread[Polygon[Join[#1, Reverse[#2]]]&, 
                  Partition[Append[#, First[#]], 2, 1]& /@ {l, L} ]}]

projectedFacesGraphics3D[{d1_, d2_, d3_}, opts___] :=
Module[{dim = Length[d1], 𝒹1, 𝒹2, 𝒹3, verticesProjected},
 (* three orthogonal vectors *)
 𝒹1 = 𝓃[d1]; 𝒹2 = 𝓃[d2 - 𝒹1.d2 𝒹1]; 𝒹3 = 𝓃[d3 - 𝒹1.d3 𝒹1 - 𝒹2.d3 𝒹2]; 
 (* the projected vertices *)
 verticesProjected = {𝒹1.#, 𝒹2.#, 𝒹3.#}& /@ vertices[dim]; 
Graphics3D[{{SurfaceColor[Hue[Random[]], Hue[Random[]], 3 Random[]], 
             EdgeForm[], makeHoled3DPolygon[#, 6, 0.82]}& /@ 
             Map[verticesProjected[[#]]&, faces[dim], {-1}],
            {GrayLevel[0], PointSize[0.02], Point /@ verticesProjected}},
      opts, PlotRange -> All, BoxRatios -> Automatic]]
      
projectedFacesGraphics3D[{min_?NumericQ, vals_}] := 
projectedFacesGraphics3D[Partition[Last /@ vals, Length[vals]/3]]

These are the resulting 3D projections. We use only the "NelderMead" method.

Show[GraphicsArray[
    Table[projectedFacesGraphics3D[
      NMinimize[relativeVertexPositionErrorSumCF[
                        Table[d[j][k], {j, 3}, {k, dim}]],
             Flatten[Table[d[j][k], {j, 3}, {k, dim}]],
             PrecisionGoal -> 2, Method -> "NelderMead"]], {dim, 4, 5}]]]

We end with minimizing the area differences of all 2D faces in 2D projections. The function relativeErrorSumArea calculates the cumulative relative face area deviation.

(* the main compiled function *)
cfArea[dim_] := cfArea[dim] = 
Compile[{{d1, _Real, 1}, {d2, _Real, 1}},
Module[{𝒹1, 𝒹2, verticesProjected, p1, p2, p3, p4,
        a, b, c, 𝒜1, 𝒜2, area = 0., 𝒱s = #1, ℱs = #2},
 (* make orthogonal vectors *)
 𝒹1 = #/Sqrt[#.#]&[d1]; 𝒹2 = #/Sqrt[#.#]&[d2 - 𝒹1.d2 𝒹1]; 
 (* the projected vertices *)
 verticesProjected = {𝒹1.#, 𝒹2.#}& /@ 𝒱s; 
 (* sum of squares of face area errors *)
 Sum[{p1, p2, p3, p4} = Map[verticesProjected[[#]]&, ℱs[[i]], {-1}];
     (* split quadrilateral into two triangles *)
     {a, b, c} = Sqrt[#.#& /@ {p1 - p2, p2 - p3, p3 - p1}]; 
     𝒜1 = Sqrt[(-a + b + c) (a + b - c) (a - b + c) (a + b + c)];
     {a, b, c} = Sqrt[#.#& /@ {p3 - p4, p4 - p1, p1 - p3}];
     𝒜2 = Sqrt[(-a + b + c) (a + b - c) (a - b + c) (a + b + c)]; 
     area = (𝒜1 + 𝒜2)/4;                         
     (1/area - 1)^2, {i, Length[ℱs]}]]]&[vertices[dim], faces[dim]]
     
(* interface function to the compiled function *)
relative2DFaceAreaErrorSumCF[
        {d1_?(VectorQ[#, NumberQ]&), d2__}] := cfArea[Length[d1]][d1, d2];

This time, we obtain very symmetric projections. Again, the minima and the resulting projections agree for the two methods.

Do[Show[GraphicsArray[projectedEdgesGraphics[
 NMinimize[{relative2DFaceAreaErrorSumCF[Table[d[j][k], {j, 2}, {k, dim}]]},
           Flatten[Table[d[j][k], {j, 2}, {k, dim}]],
           PrecisionGoal -> 3, Method -> #]]& /@ 
                        {"SimulatedAnnealing", "NelderMead"}]], 
    {dim, 3, 6}]