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.

Symbolics Volume, Chapter 3

The average length (OverBar[(Subscript[x, d]^(<))])=(min_(1≤j≤d))|(Subscript[x, j])| of the absolute value of the smallest component of a point {(Subscript[x, 1]),(Subscript[x, 2]),…,(Subscript[x, d])} on a dD unit sphere is [✶]

(OverBar[(Subscript[x, d]^(<))])=((Γ(d/2) )/(Γ((d+1)/2)))(∫_0^∞)(((1-erf(y)))^d)dy.

For each 1≤d≤10, calculate (10^6) random points and the resulting average of the absolute value of the smallest component.

A random point (with respect to the area) on a dD unit sphere can be calculated in the following recursive way [✶]. Given a (random) point {(Subscript[x, 1]),(Subscript[x, 2]),…,(Subscript[x, d-1])} on a (d-1)D unit sphere, the point {r (Subscript[x, 1]),r (Subscript[x, 2]),…,r (Subscript[x, d-1]),±(Sqrt[1-(r^2)])} is a random point on a dD unit sphere. Here ± is chosen with equal probability to generate a point on the upper/lower hemisphere of the dD unit sphere. Here r=(Subscript[f, d])(𝕣) where 𝕣 is a uniformly in [0,1] distributed variable and the inverse function (Subscript[f, d]^(-1))(𝕣) is defined through

(Subscript[f, d]^(-1))(ξ)=(((∫_0^ξ)(((ρ^(d-2))(1-(ρ^2)))^(-1/2))dρ)/((∫_0^1)(((ρ^(d-2))(1-(ρ^2)))^(-1/2))dρ)).

Starting with {±1} (± with equal probability), this procedure generates randomly distributed points on a dD unit sphere.

(* avoid messages due to singularity at ρ == 1 *)
Off[NDSolve::ndsz]; Off[InterpolatingFunction::dmval];

Module[{pp = 200, Ε = 10^100, fInv, fInvSol, ρ}, 
Do[(* solve for f_d *)
   fInvSol = NDSolve[{fInv'[ρ] == Which[ρ == 0, 0, ρ == 1, Ε, 
                                        True, ρ^(d - 2)/Sqrt[1 - ρ^2]], 
                      fInv[0] == 0}, 
                  fInv, {ρ, 0, 1}, PrecisionGoal -> 16,
                  AccuracyGoal -> 16, WorkingPrecision -> 30,
                  MaxSteps -> 10^4][[1, 1, 2]];
   FInv1 = fInvSol[1];
   (* interpolate inverted data *)
   f[d] = Interpolation[Table[{If[ρ == 1, 1, fInvSol[ρ]/FInv1], ρ} // N, 
                              {ρ, 0, 1, 1/pp}]],
   {d, 2, 21}]]

Here are the (Subscript[f, d])(ξ) for 2≤d≤21.

Plot[Evaluate[Table[f[d][ρ], {d, 2, 21}]], {ρ, 0, 1}, 
     PlotRange -> All, PlotStyle -> Table[Hue[n/26], {n, 2, 21}]]

Given a vector x on the (d-1)D unit sphere, the function liftStep generates a vector on the dD unit sphere.

liftStep[x_List] := 
With[{r = f[Length[x] + 1][Random[]]}, 
     Join[r x, {(2. Random[Integer] - 1) Re[Sqrt[1 - r^2]]}]]

Nesting the function liftStep yields the function random𝒩SpherePointList. For a given positive integer d, it returns a list of random points on a 1D sphere, 2D sphere, …, dD sphere.

random𝒩SpherePointList[d_] := 
       NestList[liftStep, {2. Random[Integer] - 1.}, d - 1]

The ensemble average of such points is uniformly distributed with respect to the surface area measure of the spheres. The following graphics show 1000 points on a circle and a classical sphere.

With[{pointList = Table[random𝒩SpherePointList[3], {1000}]}, 
Show[GraphicsArray[
     {Graphics[Point[#[[2]]]& /@ pointList,
               AspectRatio -> Automatic,
               PlotRange -> 1.66 {{-1, 1}, {-1, 1}}],
      Graphics3D[Point[#[[3]]]& /@ pointList,
                 PlotRange -> 1.1 {{-1, 1}, {-1, 1}, {-1, 1}}]}]]]

The function randomxMinList returns a list of the minimal length of the components of a random point on a 1D sphere, 2D sphere, …, dD sphere.

randomxMinList[d_] := Min[Abs[#]]& /@ random𝒩SpherePointList[d]

Modelling now (10^6) random points on spheres up to dimension 10 can be carried out in a few minutes.

dMax = 10; o = 10^6; 
SeedRandom[1];
(averages = Sum[randomxMinList[dMax], {o}]/o);  // Timing

The resulting averages agree excellently with the theoretical predictions.

(* theoretical average value *)
xMinTh[d_] := Gamma[d/2]/Gamma[(d + 1)/2]*
              NIntegrate[(1 - Erf[y])^d, {y, 0, Infinity},
                         PrecisionGoal -> 12]

(* list of {dimension, error in percent} *)
MapIndexed[{#2[[1]], #1}&, 
           Round[10000 (Table[xMinTh[d], {d, dMax}]/
                              averages - 1)]/100. "%"]

For a computational more efficient way to produce random points on high-dimensional spheres, see [✶].