8
$\begingroup$

I am looking for an alternative and robust alternative to Gram-Schmidt orthogonalization, but with one constraint:

I have a unit vector $\mathbf{v}_1 \in \mathbb{R}^d \,\text{s.t.} \,\|\mathbf{v}_1^T\mathbf{v}_1\|=1$ and I want to find a set of $d-1$ vectors $\{\mathbf{v}_2\cdots \mathbf{v}_{d}\} \perp \mathbf{v}$, orthogonal to $\mathbf{v}$.

Therefore, my goal is not to generate a set of arbitrary vectors, but to find $d-1$ other vectors which form a $d \times d$ square matrix $\mathbf{V}$ with orthonormal columns: $$ \mathbf{V}= \begin{bmatrix} \mathbf{v}_1 & \mathbf{v}_2 & \cdots & \mathbf{v}_d \end{bmatrix} $$

Gram-Schmidt and modified Gram-Schmidt algorithms can accomplish that, but it is known that even modified Gram-Schmidt suffers from numerical inaccuracy. There are methods for doing this in 3D such as this one. However, as cross-product does not generalize to n-dimensions we cannot extend the simple 3D methods.

I am also aware of the similar posts such as this one but my question is different as the first vector is specified.

Ideally, I am looking for a computationally efficient and accurate method to perform this task. Simplicity is important because I will also like to compute the gradient of this operation.

$\endgroup$
2
  • 1
    $\begingroup$ Is MGS really not good enough? The problem statement almost insists upon this algorithm (seeking explicit tabulation of Q and want to include some specific v in Q). Have you tried taking a random DxD matrix, overwriting the first column with v, and MGS'ing it? How did it fail? Alternatively, you could make a random Dx(D-1) matrix, project v off from each column, then use some other QR algorithm (householder, perhaps column pivoting, etc). You would reconstruct V/Q explicitly in a second pass (dorgqr-like routine, or apply Q to identity using dormqr-like routine) $\endgroup$ Commented Sep 12, 2017 at 20:31
  • $\begingroup$ This is similar to @Will P.'s comment. I do not suppose that it is good enough. That's the reason why I ask further opinions. MGS didn't fail yet though. $\endgroup$ Commented Sep 12, 2017 at 20:38

4 Answers 4

14
$\begingroup$

There is a known mathematical question here: you are given a unit vector, lying on the $(n-1)$-sphere in $\mathbb{R}^n$, $v\in S^{n-1}$, and you would like to associate with each such vector a frame in the tangent bundle of $S^{n-1}$ at $v$. This is a map $S^{n-1}\to \mathrm{F}S^{n-1}$ from the sphere to its frame bundle, also known as a global section of the frame bundle of $S^{n-1}$.

Here, the frame bundle (https://en.wikipedia.org/wiki/Frame_bundle) is the set of all pairs $(v, (u_1, \ldots, u_{n-1}))$, of a unit vector $v$ and $n-1$ orthonormal vectors all orthogonal to $v$, and a global section means a map $v \mapsto (u_1, \ldots, u_{n-1})$ that is defined for all $v$ (https://en.wikipedia.org/wiki/Section_(fiber_bundle)#Extending_to_global_sections).

Manifolds that can have such a map are known as parallelizable (https://en.wikipedia.org/wiki/Parallelizable_manifold), and as Wikipedia says, among the unit spheres only $S^0$, $S^1$, $S^3$ and $S^7$ are parallelizable. So if you are looking for a way to get the orthogonal complement $v^\perp$ of $v$ that is valid and smooth in $v$ for all $v$, you can only do that for $n=1,2,4,8$.

For $n=3$, in fact, the proof of impossibility is simple: if the two orthogonal vectors such an algorithm would compute are $f(v)$ and $g(v)$, and both depend continuously on $v$, then $f(v)$ defines a tangent vector at $v$ to the unit sphere at the point $v$, and such a tangent vector field would contradict the hairy ball theorem (https://en.wikipedia.org/wiki/Hairy_ball_theorem).

For $n=2$, you have the map $(a,b) \mapsto (-b, a)$.

For $n=4$, you have the map $$ (a, b, c, d) \mapsto \begin{pmatrix} -b & -c & -d\\ a & d & -c\\ -d & a & b\\ c & -b & a \end{pmatrix} . $$ This corresponds to writing down the quaternion $w = a+\mathrm{i}b+\mathrm{j}c+\mathrm{k}d$, where the dot product is computed as $v_1^\top v_2 =\Re(\bar w_1 w_2)$), and the three orthogonal vectors are $\mathrm{i}w$, $\mathrm{j}w$, $\mathrm{k}w$.

For $n=8$, you can read the answer off of the Cayley table for octonions (https://en.wikipedia.org/wiki/Octonion).

For $n=3$ and other $n$, I think the best you can do is to pick a "pole", such as the first basis vector $e_1$, choose the rotation matrix that rotates in the $(e_1, v)$-plane only and maps $Re_1 = v$, and use the $(n-1)$-frame $(Re_2, \ldots, Re_n)$. The matrix $R$ can be computed as a sequence of Givens rotations, which means its directional derivatives should be straightforward to compute. This map has a singularity at $v=-e_1$. This is equivalent to parallel-transporting a tangent frame at $e_1$ from $e_1$ to $v$ along the shortest path (a great circle on $S^{n-1}$) from $e_1$ to $v$, so there is a geometric interpretation to doing this.

More explicitly, write $v = \cos\theta e_1 + \sin\theta e_v$, where $e_v$ is a unit vector, $e_1^\top e_v = 0$. Then the rotation matrix in the $(e_1,e_v)$-basis is $$R_0 = \begin{pmatrix}\cos\theta&-\sin\theta\\\sin\theta&\cos\theta\end{pmatrix},$$ and the change of basis from $(e_1,e_2,\ldots)$ to $(e_1,e_v,\ldots)$ is given by any orthonormal matrix $X$ such that $Xe_1 = e_1$, $Xe_2 = v$. This can be computed by QR decomposition of the $n\times(n+1)$ matrix with columns $(e_1,v,e_2,\ldots,e_n)$, and the answer is the matrix $R_0$ with the change-of-basis formula applied to it, $$ R = X \begin{pmatrix} R_0&0\\0&I \end{pmatrix} X^{-1}. $$ This is independent of which basis $X$ is chosen, because however the QR decomposition completes the basis $R$ is just the identity on it, and it would also only fail on the two special cases, $v=\pm e_1$.

$\endgroup$
5
  • $\begingroup$ This is amazing. I am tending towards accepting this answer due to the simplicity and the mathematical thoroughness. One of my cases was exactly quaternions ($d=4$) and this is properly addressed in your reply. If you like to elaborate more for the properties of the case of $d=4$ that would be very welcome. $\endgroup$ Commented Sep 12, 2017 at 23:38
  • $\begingroup$ @TolgaBirdal Using quaternions to represent 4d, and unit quaternions to represent the unit sphere is a pretty rich topic, in particular it's well-used in computer graphics to represent rotations (because the rotation group in 3d is double-covered by the 3-sphere, which is the group unit quaternions), so anything I can say more is probably already said even on Wikipedia (en.wikipedia.org/wiki/Quaternion). $\endgroup$
    – Kirill
    Commented Sep 12, 2017 at 23:42
  • $\begingroup$ Even though I am kinda aware of the quaternions and parallel transport theorem, such a map wasn't immediate for me. Thanks for pointing that out. $\endgroup$ Commented Sep 12, 2017 at 23:45
  • $\begingroup$ @TolgaBirdal I think your question is very nice: at first I thought the difficulty of computing this was numerical, but it turned out to be purely mathematical instead. $\endgroup$
    – Kirill
    Commented Sep 12, 2017 at 23:51
  • $\begingroup$ Glad to hear that. Your reply was a great mind opener and of huge help. I appreciate a lot. $\endgroup$ Commented Sep 12, 2017 at 23:59
14
$\begingroup$

In order to find such a set of vectors, you can use the Householder QR factorization. Let your unit vector $v$ be given. Define a nonsingular matrix $A$ by $$ A = \left(\begin{array}{cccc} v & a_2 & \cdots & a_d \end{array}\right), $$ where it is not so important how you obtain the columns $a_2, \ldots, a_d$ of $A$, as long as $\{v, a_2, \ldots, a_d\}$ forms a linearly independent set of vectors.

Then, the QR factorization of $A$ results in an orthogonal matrix $Q$ such that $$ A = QR, $$ where $R$ is an upper triangular matrix. We can write $$ R = \left( \begin{array}{cccc} r_{11} & \times & \cdots & \times \\ 0 & \times & \cdots & \times \\ %0 & 0 & \cdots & \times \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \times \end{array}\right). $$

Since $Q$ is an orthogonal matrix, its columns, which we will denote $q_j$, form an orthonormal set of vectors. The relationship $$ Q^T A = R $$ tells us that $q_1^T v = r_{11}$, and $q_j^T v = 0$ for $j > 1$. Therefore, we have found an orthonormal set of vectors $$ \{ v, q_2, \ldots, q_d \}.$$

In fact, since both $q_1$ and $v$ are orthogonal to the remaining $(d-1)$ vectors $q_j$, we have that $q_1$ and $v$ belong to the same one-dimensional subspace and hence are collinear. Since they are both unit vectors, we know $q_1 = \pm v$.

Regarding the numerical stability of this algorithm, if the QR factorization is performed using Householder reflectors, then $Q$ is generated as a product of orthogonal matrices $Q_d Q_{d-1} \cdots Q_2 Q_1$ of the form $I - 2 w_j w_j^T$. The product of orthogonal matrices is numerically stable, and hence we can expect that the rows of $Q$ will indeed be numerically orthogonal (up to machine precision).

$\endgroup$
6
  • 7
    $\begingroup$ The difficult part, however, is to deal with the statement "where it is not so important how you obtain the columns $a_2,…,a_2,…,a_d$" :-) $\endgroup$ Commented Sep 12, 2017 at 2:08
  • $\begingroup$ @WolfgangBangerth Maybe search a $k$ such that $|v_k|=max_i |v_i|$ and then use the standard basis unit vectors $e_1, ..., e_{k-1}, e_{k+1}, ..., e_d$ ? $\endgroup$
    – wim
    Commented Sep 12, 2017 at 10:47
  • $\begingroup$ @wim can you give more details about this? It would be good for completeness, as this is an important stage. I am still willing to hear other opinions though, because QR is rather complex and not really suitable for differentiation for example. How about a geometric approach to this? $\endgroup$ Commented Sep 12, 2017 at 12:12
  • $\begingroup$ I've put my reply into a separate answer to the original question. $\endgroup$ Commented Sep 12, 2017 at 13:07
  • 2
    $\begingroup$ I think you can simplify this procedure if you use the $n\times(n+1)$ matrix $A=(v,e_1,\ldots,e_n)$ instead, with column vectors $v$ and any $n$ basis vectors. It saves you the trouble of figuring out if the $n-1$ basis vectors just happen to span $v$ or not. $\endgroup$
    – Kirill
    Commented Sep 13, 2017 at 0:12
4
$\begingroup$

How about doing an SVD to get the null space of the $1 \times d$ matrix $A = [\bar{v}]$ with $\bar{v} := v^T/\|v\|$? The null function in MATLAB does exactly this:

d = 4; v = rand(d,1); vbar = v.'/norm(v); Q = null(vbar).'

Then your orthonormal vectors are the rows of $ \begin{bmatrix} \bar{v} \\ Q \end{bmatrix}.$

$\endgroup$
1
  • $\begingroup$ That is true. In fact, the MATLAB $\textit{ortho}$ function does it in a similar way through SVD orthogonalization. Again, a complex operation, but worth mentioning. $\endgroup$ Commented Sep 12, 2017 at 20:39
1
$\begingroup$

You can achieve something like what you are looking for by taking the equivalent of the cross product to keep producing vectors that are perpendicular to the ones you already have.

If you're only interested in the 2d or 3d cases, here is a function that does this (including explanation of how it works): https://github.com/geodynamics/aspect/blob/master/source/utilities.cc#L511

The documentation of the function is here: https://github.com/geodynamics/aspect/blob/master/include/aspect/utilities.h#L213

$\endgroup$
2
  • $\begingroup$ Actually I'm interested in N-d case (well at least 4d). Cross product is not defined in higher dimensions properly I believe. Could you point me to an equivalent stable extension of that, if you know? Maybe this one: mathoverflow.net/questions/94312/… $\endgroup$ Commented Sep 12, 2017 at 20:41
  • $\begingroup$ I haven't thought about the higher dimensional case. But just like in the 3d case, where you have to have come up with a vector that "is definitely not collinear with $v$", you'll just have to come up with a total of $n-2$ such vectors. It should be possible to construct them similarly to the approach used there, going down the list of largest entries of $v$ (or making something up for the small entries if it only has one large entry). $\endgroup$ Commented Sep 12, 2017 at 21:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.