6
$\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$
  • 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$ – rchilton1980 Sep 12 '17 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$ – Tolga Birdal Sep 12 '17 at 20:38
11
$\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$
  • $\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$ – Tolga Birdal Sep 12 '17 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 Sep 12 '17 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$ – Tolga Birdal Sep 12 '17 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 Sep 12 '17 at 23:51
  • $\begingroup$ Glad to hear that. Your reply was a great mind opener and of huge help. I appreciate a lot. $\endgroup$ – Tolga Birdal Sep 12 '17 at 23:59
11
$\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$
  • 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$ – Wolfgang Bangerth Sep 12 '17 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 Sep 12 '17 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$ – Tolga Birdal Sep 12 '17 at 12:12
  • $\begingroup$ I've put my reply into a separate answer to the original question. $\endgroup$ – Wolfgang Bangerth Sep 12 '17 at 13:07
  • $\begingroup$ Probably it is also necessary to multiply the unit vectors by $max_i|v_i|$ in order to obtain a well conditioned matrix $A$. In that case it might be possible to derive an upper bound on the condition number of $A$ (?, I have to think about this). For well conditioned matrices even a slightly altered version of classical Gram-Schmidt (according to Åke Björck's Numerical Methods in Matrix Computations) is suitable. $\endgroup$ – wim Sep 12 '17 at 14:41
3
$\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$
  • $\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$ – Tolga Birdal Sep 12 '17 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$
  • $\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$ – Tolga Birdal Sep 12 '17 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$ – Wolfgang Bangerth Sep 12 '17 at 21:44

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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