Consider an equation $S(t)b(t) = a$, where $a, b(t) \in S^{n-1}$ are given and the vector $b(t)$ is continuous, i.e. its endpoint traces a continuous curve on the unit sphere. The task is to find continuous solution $S(t) \in SO(n)$ (numerically). Consider a time net $t_0 < ... < t_n$.

  1. In 2D case we can count the angle of $b(t)$ on the unit circle and then $S(t)$ will be the rotation matrix by this angle (inversed) corrected by the angle of $a$.
  2. In 3D case we can use the Euler theorem and correct the matrix $S(t_k)$ on each step of our iterational method by multiplying it by the rotation matrix $S'$ that sends $b(t_k)$ to $b(t_{k+1})$.
  3. In general case assume that $a=(1,0,...,0)^{T}$. I tried to construct some orthonormed basis ${b(t_0),e_2,...,e_n}$ at first step and to let $S(t_0) = [b(t_0),e_2,...,e_n]^{-1}$ using Gram-Schmidt process. Then at each step of our method I constructed the new one using Gram-Schmidt process and the previous basis as initial condition. It works, but it is time-consuming enough.

My question is how to obtain such $S(t)$ in general case more quickly?

  • $\begingroup$ Could you describe what are the spaces $S^{n-1}$ and $SO(n)$? $\endgroup$
    – Paul
    Mar 21, 2012 at 14:00
  • $\begingroup$ $SO(n)$ is the special orthogonal group of dimension $n$ (I'm assuming over the reals), which is the set of matrices $Q \in \mathbb{R}^{n \times n}$ such that $QQ^{T} = Q^{T}Q = I$, and $\det(Q) = 1$. $S^{n-1}$ is the $(n-1)$-dimensional unit sphere. $\endgroup$ Mar 21, 2012 at 14:32
  • $\begingroup$ I'm curious about the description of method 3. When you say that you construct $S(t_{k+1})$ from $S(t_{k})$, are you applying Gram-Schmidt to $[b(t_{k+1}), e_{2}, \ldots, e_{n}]^{-1}$ again? What do you mean by "[using] the previous basis as initial condition"? $\endgroup$ Mar 21, 2012 at 14:36
  • $\begingroup$ What's preventing you from using the infinitesimal generator of rotation in the general case? $\endgroup$ Mar 21, 2012 at 15:23
  • $\begingroup$ @GeoffOxberry, I mean that I apply the Gram-Schmidt process to $b(t_{k+1}), e_2(t_k),...,e_n(t_{k})$ and I receive the basis $b(t_{k+1}), e_2(t_{k+1}),...,e_n(t_{k+1})$. $\endgroup$
    – Appliqué
    Mar 21, 2012 at 15:53

2 Answers 2


Let $A(t_{k}) = [b(t_{k}), e_{2}, \ldots, e_{n}]$. If $b(t_{k+1}) = b(t_{k})$, nothing changes, so without loss of generality, assume that $b(t_{k+1}) - b(t_{k}) \neq 0$, and also assume that $A(t_{k})$ and $A(t_{k+1})$ are both invertible.

You can get $A(t_{k+1})^{-1}$ from $A(t_{k})^{-1}$ by noting that $A(t_{k+1})$ is a rank-one update of $A(t_{k})$, and using the Sherman-Morrison formula, with $u = b(t_{k+1}) - b(t_{k})$, and $v = e_{1}$; this rank-one update would be cheaper than naïvely inverting $A(t_{k+1})$ using an LU decomposition.

Then, since $A(t_{k+1})^{-1}$ is a rank-one update of $A(t_{k})^{-1}$, you can update your QR factorization (at least, that's how I would do Gram-Schmidt numerically) also using a rank-one updating scheme. One algorithm for accomplishing a rank-one update of a QR factorization can be found in Section 12.5.1 of the third edition of Matrix Computations by Golub and van Loan.

Both of these updates should reduce the complexity of $\mathcal{O}(n^{3})$ operations to $\mathcal{O}(n^{2})$ operations; you can even calculate $A(t_{0})^{-1}$ and its QR factorization in this fashion because $A(t_{0})$ is a rank-one update of the identity matrix.

  • $\begingroup$ Hm, $S(t_{k}) = [b(t_{k}),e_1(t_{k}),...,e_n(t_k)]^{T}$. Then I have to see only for update of decomposition $[b(t_{k+1}),e_1(t_{k}),...,e_n(t_k)] = QR$ and take $S(t_{k+1}) = Q^{T}$. Is it true? $\endgroup$
    – Appliqué
    Mar 21, 2012 at 17:47
  • $\begingroup$ That's another way you could do it. The first vector in your basis should just be scaled in a QR decomposition, so once you find the first QR decomposition (i.e., the identity matrix), you can repeatedly rank-one update it. $\endgroup$ Mar 21, 2012 at 18:11
  • $\begingroup$ So if we receive $S(t_{k+1})$ using orthoghonalisation of $[b(t_{k+1}),e_2(t_k),...,e_n(t_k)]$ by one hand and updating a QR decomposition by the other we will get two different resulting matrices? Why the last will sent $b(t_{k+1})$ to $[1,0,...,0]^{T}$ and will be close to $S(t_k)$? I think that update algorithm doesn't garantee this $\endgroup$
    – Appliqué
    Mar 21, 2012 at 18:29
  • $\begingroup$ I don't know that the update algorithm will guarantee that the first column of the updated QR decomposition will be a scaled version of $b(t_{k+1})$, but the first update to the QR decomposition would be correct. $\endgroup$ Mar 21, 2012 at 18:36
  • $\begingroup$ I tried the updating in MATLAB; it won't guarantee that the first column of the updated QR decomposition is a scaled version of $b(t_{k+1})$. $\endgroup$ Mar 21, 2012 at 19:04

Rewrite your equation $$ b(t)=B(t)a $$

$B(t+\Delta t)=B^\prime B(t)$ where $$B^\prime=\left(\begin{matrix} b/\|b\| & \dot b/\|\dot b\| \end{matrix}\right)\left(\begin{matrix} \cos \theta & \sin \theta \\ -\sin \theta & \cos \theta \end{matrix}\right)\left(\begin{matrix} b/\|b\| & \dot b/\|\dot b\| \end{matrix}\right)^\ast +\\ I-\left(\begin{matrix} b/\|b\| & \dot b/\|\dot b\| \end{matrix}\right)\left(\begin{matrix} b/\|b\|& \dot b/\|\dot b\| \end{matrix}\right)^\ast$$ and $\theta$ such that $B^\prime$ takes $b(t)$ to $b(t+\Delta t)$. Now $S(t)=B(t)^\ast$.

B was constructed such that your coordinate frame is defined by $b$ and $\dot b$. Since all other directions are unaffected, we project onto the complement of $\mbox{span }\{b,\dot b\}$ (second term on the right) and then rotate within $b$ and $\dot b$ (first term on right).

You may also want to consider that your problem can be understood as

$$ S\dot b +\dot S b = 0$$ such that $$ S\in SO(n).$$

  • $\begingroup$ It's important to note that the differential condition is necessary, but not sufficient. $\endgroup$ Mar 21, 2012 at 19:41
  • $\begingroup$ @GeoffOxberry: Yes, I've omitted the initial conditions as well, and of course any rotation on the complement of span $\{b,\dot b\}$ would still be a solution. Obviously $B^\prime$ is not unqiue. It would work without the projector as well. I just prefer having invariant spaces over null spaces :-) $\endgroup$ Mar 21, 2012 at 19:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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