How to fix errors in Gram-Schmidt process when using random vectors?

Question

I first make a function to get a random vector on unit sphere in a swath around the equator. That is what the parameter $\gamma$ controls; if $\gamma = 1/2$, the vectors can be chosen anywhere on the sphere; if $\gamma = 0$, the vectors can only be chosen from the equatorial plane etc.

randomVectorOnUnitSphere[\[Gamma]_] := Module[{},
\[Phi] = RandomReal[{0, 2 \[Pi]}];
w = RandomReal[{1/2 - \[Gamma], 1/2 + \[Gamma]}];
nx = Cos[\[Phi]] Sin[ArcCos[2 w - 1]];
ny = Sin[\[Phi]] Sin[ArcCos[2 w - 1]];
nz = Cos[ArcCos[2 w - 1]];
Return[{nx, ny, nz}];
]

These two random vectors span a plane. I would like to find an orthonormal basis for that plane so I use the Gram-Schmidt process as so:

randMod[\[Gamma]_] := Module[{},
v1 = randomVectorOnUnitSphere[\[Gamma]];(*random unit vector 1*)
v2 = randomVectorOnUnitSphere[\[Gamma]];(*random unit vector 2*)
u1 = v1; (* start Gram-Schmidt process *)
proj = u1.v2/u1.u1 u1;
u2temp = v2 - proj;
u2 = u2temp/Sqrt[u2temp.u2temp];
Return[{u1,u2}] (*the output here is 2 orthonormal vectors that span the plane*)
];

I look at 10000 pairs and plot on unit sphere:

randTable1 = ParallelTable[randMod[1/10][[1]], {i, 1, 10000}];
randTable2 = ParallelTable[randMod[1/10][[2]], {i, 1, 10000}];
randPlot1 = 
ListPointPlot3D[{randTable1, {{0, 0, 1}}}, 
PlotStyle -> {Black, {PointSize -> 0.025, Red}}];
randPlot2 = 
ListPointPlot3D[{randTable2, {{0, 0, 1}}}, 
PlotStyle -> {Black, {PointSize -> 0.025, Red}}];
{Show[Graphics3D[Sphere[{0, 0, 0}, 1]], randPlot1], 
Show[Graphics3D[Sphere[{0, 0, 0}, 1]], randPlot2]} 

and get the output for all the $u_1$ (top) and $u_2$ (bottom):

The $u_2$ vectors are concentrated at the equator but there is a sizable number that are spreadout over the sphere. I would expect to get something that looked just like $u_1$ instead.

The initial vectors $v_1$ and $v_2$ lie in the same 'swath' as $u_1$ so I know that it is the G-S process that is failing.

I have looked at some of the culprit points and have noticed that their initial starting vectors were near parallel. I think this is the cause of the problem. I tried increasing WorkingPrecision when I call RandomReal but that did not seem to have any noticeable effect. Is there a way I can fix this problem?

Answer 1

I'm not sure why you'd expect the resulting vectors to remain in the original swath because you're forcing the angle between the two resulting vectors to be 90 deg. Pairs of vectors with cross products far from the poles of the sphere will necessarily result in one of the vectors "popping out" of the swath.

I'd also like to point out that you don't need Gram-Schmidt to create an orthonormal basis for the two vectors. Instead, you can just do the following: Given two vectors v1 and v2, define u1=v1 then construct u2=Normalize[v1 x v2]. Finally, u3 = u2 x u1. Now u1 will be aligned with v1 and u3 will be in the plane defined by v1 and v2 but it will be orthogonal to u1 and in the general direction of v2.