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

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?

-
Why don't you expect vectors u2 away from the equatorial plane? Suppose you start with v1={1,0,0} and v2={Sqrt[1-eps^2],0,eps} where eps is small, then u2 will be {0,0,1}. – Heike Feb 3 '12 at 21:18

You could also construct a basis by using Orthogonalize[{v1,v2, Cross[v1,v2]}] – Heike Feb 3 '12 at 21:11