How to define a set of two random orthogonal unit vectors?

Question

I need to define a random set of two orthogonal unit vectors. The code below doesn't give proper unit and orthogonal vectors :

Clear["Global`*"]

n[k_] := Normalize[{RandomReal[{-1, 1}], RandomReal[{-1, 1}], RandomReal[{-1, 1}]}]

u[k_] := Normalize[{RandomReal[{-1, 1}], RandomReal[{-1, 1}], RandomReal[{-1, 1}]}]

a[k_] := Normalize[u[k] - (u[k].n[k])n[k]]

(* unitary and orthogonality test : *)
Table[n[k].n[k], {k, 1, 5}]  (* this should output five 1 *)
Table[a[k].a[k], {k, 1, 5}]  (* this should output five 1 *)
Table[n[k].a[k], {k, 1, 5}]  (* this should output five 0 *)

Specifically, for any value of $k$ (which should be any real number) : \begin{align} \vec{\boldsymbol{n}}(k) \cdot \vec{\boldsymbol{n}}(k) &= 1, \\ \vec{\boldsymbol{a}}(k) \cdot \vec{\boldsymbol{a}}(k) &= 1. \\ \vec{\boldsymbol{n}}(k) \cdot \vec{\boldsymbol{a}}(k) &= 0. \end{align}

So what is wrong with the code above ?

n[k_]:= defines a function of k. If you call it 3 times with the same argument, it gets evaluated 3 times and you get 3 different random vectors — Niki Estner, Mar 23 '16 at 15:30

Martin Ender · Accepted Answer · 2016-03-23T16:05:35.343

As nikie pointed out, the issue is that each time you call one of your functions (regardless of k) new random vectors will be calculated. A quick fix for this is to memoise the functions. That means they'll be assigned a random value the first time they are called, but will then remember that value for the given k and always return the same. Memoisation is trivial in Mathematica:

n[k_] := n[k] = Normalize[{RandomReal[{-1, 1}], RandomReal[{-1, 1}], RandomReal[{-1, 1}]}]
u[k_] := u[k] = Normalize[{RandomReal[{-1, 1}], RandomReal[{-1, 1}], RandomReal[{-1, 1}]}]
a[k_] := a[k] = Normalize[u[k] - (u[k].n[k]) n[k]]

Table[n[k].n[k], {k, 1, 5}]
Table[a[k].a[k], {k, 1, 5}]
Table[n[k].a[k], {k, 1, 5}]

(* {1., 1., 1., 1., 1.} *)
(* {1., 1., 1., 1., 1.} *)
(* {9.4369*10^-16, 2.77556*10^-17, 5.55112*10^-17, 0., 5.55112*10^-17} *)

Note that the slight discrepancies from zero are due to floating point inaccuracies. These are unavoidable due to the finite representation of "real" numbers in computers. I recommend that you give the floating point guide a thorough read. If you're feeling brave you can also read What Every Computer Scientist Should Know About Floating-Point Arithmetic, but it's a lot more technical.

When working with floats, you should always check that your results fall into a narrow band of tolerance (e.g. $10^{-10}$) of the desired value and never use exact equality.

As for your actual implementation, n and u can be simplified by generating a random vector right away:

n[k_] := n[k] = Normalize[RandomReal[{-1, 1}, 3]]
u[k_] := u[k] = Normalize[RandomReal[{-1, 1}, 3]]

Alternatively, you can use Orthogonalize, where now a and b are your orthonormal vectors, but note that the memoisation gets a bit trickier, since you're generating two corresponding vectors when either one is requested:

a[k_] := First[{a[k], b[k]} = Orthogonalize[RandomReal[{-1, 1}, {2, 3}]]]
b[k_] := Last[{a[k], b[k]} = Orthogonalize[RandomReal[{-1, 1}, {2, 3}]]]

Table[a[k].a[k], {k, 1, 5}]
Table[b[k].b[k], {k, 1, 5}]
Table[a[k].b[k], {k, 1, 5}]
(* {1., 1., 1., 1., 1.} *)
(* {1., 1., 1., 1., 1.} *)
(* {0., 0., 3.46945*10^-17, 5.55112*10^-17, -5.55112*10^-17} *)

Also, you should note that this method (as well as your original method) will not create vectors which are uniformly distributed over the sphere. The vectors you get from Normalize[RandomReal[{-1,1},3]] will be less commonly aligned with the axes and more commonly point along the diagonals of the space. If you want to ensure a uniform distribution of these vectors, you'll have to do something more clever (rejection sampling, or using a weighted distribution in spherical polar coordinates maybe). I am also not sure whether a uniform distribution of input vectors will result in a uniform distribution of planes spanned by the two orthogonal vectors you end up with, although it seems likely on grounds of symmetry.

A simple question : Is this code the best (proper or simplest) way to define random orthogonal unit vectors ? Or is there a better way of doing this with Mathematica ? — Cham, Mar 23 '16 at 15:48
@Cham You can definitely simplify n and u with RandomReal[{-1, 1}, 3] to get a random vector. I'll edit the answer in a bit. — Martin Ender, Mar 23 '16 at 15:49
About the random unit vector uniformly distributed on a sphere, I would simply use n = {Sqrt[1 - u^2]Cos[phi], Sqrt[1 - u^2]Sin[phi], u}, with u = RandomReal[{-1, 1}] and phi = RandomReal[{0, 2Pi}]. This is automatically uniformly distributed on a sphere. So no complicated rejection code. — Cham, Mar 23 '16 at 19:41
In the case of the vector $\vec{a}$, I would use some easy rotation matrices. We have $\vec{n} = R_z(\phi) R_y(\theta) \hat{z}$. The orthogonal vector is then defined as this : $\vec{a}(\alpha) = R_z(\phi) R_y(\theta) R_z(\alpha) \hat{x}$, where $\alpha$ is a random angle uniformly distributed on $(0, , 2 \pi)$. It is easy to show that $\vec{a} \cdot \vec{n} \equiv 0$, for all angles $\alpha$. I'll then use the random variable $u = \cos{\theta}$ and $\sqrt{1 - u^2} = \sin{\theta}$. — Cham, Mar 23 '16 at 19:55
This question addresses methods of choosing vectors on a sphere with a uniform distribution. — David Z, Mar 23 '16 at 21:21

unlikely · Answer 2 · 2016-03-23T18:11:56.797

(*Function*)
nal[k_Integer?Positive] := 
 RandomVariate[UniformDistribution[{{0, \[Pi]}, {-\[Pi], \[Pi]}}], 
       2 k] // Map@Prepend[1] // 
     CoordinateTransform["Spherical" -> "Cartesian", #] & // 
    Partition[#, 2] & // Map@Orthogonalize // Transpose

(*Use*)
{nl, al} = nal[5]
Norm /@ nl
Norm /@ al
MapThread[Dot, {nl, al}] // Chop

(*Or*)
Clear[n, a]
Evaluate@{Array[n, 5], Array[a, 5]} = nal[5]
Table[n[k].n[k], {k, 1, 5}]
Table[a[k].a[k], {k, 1, 5}]
Table[n[k].a[k], {k, 1, 5}] // Chop

Results

{{{-0.825575, 0.22288, 0.518411}, {0.128707, 0.483466, 
   0.865849}, {0.959535, -0.178934, 0.21743}, {0.780066, -0.574996, 
   0.24673}, {-0.921722, 0.0682483, 0.3818}}, {{0.48972, -0.173447, 
   0.854453}, {0.0995755, 0.862398, -0.496341}, {-0.220525, 
   0.00265385, 0.975378}, {0.302404, 
   0.00124804, -0.953179}, {-0.28588, 0.545705, -0.787704}}}

{1., 1., 1., 1., 1.}

{1., 1., 1., 1., 1.}

{0, 0, 0, 0, 0}

{{{-0.50804, 0.653068, 0.561602}, {0.625118, -0.315714, 
   0.71383}, {0.748392, -0.629871, 
   0.207781}, {0.382584, -0.0387238, -0.923109}, {0.0835878, 
   0.288501, -0.953824}}, {{-0.370342, -0.754291, 
   0.542118}, {-0.135175, 0.856936, 0.497382}, {0.335259, 
   0.0889442, -0.937918}, {-0.918986, -0.119072, -0.375881}, \
{0.685458, -0.711399, -0.155105}}}

{1., 1., 1., 1., 1.}

{1., 1., 1., 1., 1.}

{0, 0, 0, 0, 0}

score 1 · Answer 3 · answered Mar 23 '16 at 18:23

another simple way to get uniformly distributed vectors, by rejecting random vectors that fall outside a unit sphere:

randv := Module[{v}, 
  Normalize[While@(Norm[v = RandomReal[{-1, 1}, 3]] > 1); v]]
rs = Table[Orthogonalize[{randv, randv}], {500}];
Norm /@ rs[[All, 1]] // Union  (*1.*)
Norm /@ rs[[All, 2]] // Union  (*1.*)
Dot @@ # & /@ rs // Chop // Union (*0.*)

David G. Stork · Answer 4 · 2016-09-25T17:43:04.990

Any of the above "solutions" that is based on uniformly sampling a unit cube and simply projecting a point onto a sphere will not give a uniform random pair of orthogonal vectors. That is because the projection does not uniformly sample the sphere. You can see that by looking at the cluster of points on the unit sphere close to the corners of the cube in this figure:

Graphics3D[{
  {Opacity[0.5], Sphere[]},
  Point[Normalize /@ RandomReal[{-1, 1}, {50000, 3}]]
  }]

Instead, one should sample from a spherical Gaussian distribution, which has no such preferred directions.

This code (truly) randomly selects a direction from the unit sphere, rotates each of the two unit orthogonal basis vector by a random angle $0 \le \theta \leq 2 \pi$ around that direction.

u = {1, 0, 0}; 
v = {0, 1, 0};
orthogonalvectors = Orthogonalize@(RotationMatrix[
       RandomReal[{0, 2 π}], 
       RandomVariate[
        MultinormalDistribution[{0, 0, 0}, 
         IdentityMatrix[3]]]].# & /@ {u, v})

(* {{0.988191, -0.112487, -0.104045}, {-0.0381163, 0.477218, -0.877958}} *)

Confirm by:

ortogonalvectors[[1]].ortogonalvectors[[1]] //Chop

(* 1 *)

ortogonalvectors[[2]].ortogonalvectors[[2]] //Chop

(* 1 *)

ortogonalvectors[[1]].ortogonalvectors[[2]] //Chop

(* 0 *)

Run the code as often as you like to get new pairs of orthogonal vectors.

Compare closely the above sphere figure with this one that properly samples the unit sphere:

Graphics3D[{
  {Opacity[0.5], Sphere[]},
  Red, Point[
   Normalize /@ 
    RandomVariate[
     MultinormalDistribution[{0, 0, 0}, IdentityMatrix[3]], 5000]]
  }]

This solution appears to be a bit similar to the one I've adopted, that I posted above. However, I don't use any gaussian distribution. — Cham, Mar 24 '16 at 02:13
@Cham, the point is that the normal distribution is the right distribution to use for uniformly sampling a sphere. Certainly cheaper than evaluating trig functions. — J. M.'s missing motivation, Mar 24 '16 at 02:31

score 1 · Answer 5 · edited Mar 10 '17 at 09:42

Here's the solution I've adopted :

(* Uniformly distributed random values of Cos[theta] : *)
u[k_] := u[k] = RandomReal[{-1, 1}]

(* Uniformly distributed random value of phi : *)
phi[k_] := phi[k] = RandomReal[{0, 2Pi}]

(* Random unit vector : *)
n[k_] := n[k] = {
        Sqrt[1 - u[k]^2]Cos[phi[k]],
        Sqrt[1 - u[k]^2]Sin[phi[k]],
        u[k]
    }

(* Uniformly distributed random angle AROUND the vector n[k] : *)
alpha[k_] := alpha[k] = RandomReal[{0, 2Pi}]

(* Random unit vector, orthogonal to n[k] by construction : *)
a[k_] := a[k] = {
        u[k]Cos[alpha[k]]Cos[phi[k]] - Sin[alpha[k]]Sin[phi[k]],
        u[k]Cos[alpha[k]]Sin[phi[k]] + Sin[alpha[k]]Cos[phi[k]],
        -Sqrt[1 - u[k]^2]Cos[alpha[k]]
    }

Here's a plot of 5000 points from the n[k] vector (in red), and 5000 points from the a[k]vector (in blue) :

random points on a sphere

Code to plot this picture :

Graphics3D[
    Table[
    {
        {Red,  PointSize -> 0.003, Point[n[k]]},
        {Blue, PointSize -> 0.003, Point[a[k]]}
    }, {k, 1, 5000}], 
SphericalRegion -> True,
Method -> {"RotationControl" -> "Globe"},
ImageSize -> {600, 600}
]

Is there any drawback to this method ?

How to define a set of two random orthogonal unit vectors?

5 Answers5