Time-stable SO(n) matrix synthesis algorithm

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?

• Could you describe what are the spaces $S^{n-1}$ and $SO(n)$?
– Paul
Mar 21, 2012 at 14:00
• $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. Mar 21, 2012 at 14:32
• 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"? Mar 21, 2012 at 14:36
• What's preventing you from using the infinitesimal generator of rotation in the general case? Mar 21, 2012 at 15:23
• @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})$. Mar 21, 2012 at 15:53

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.

• 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? Mar 21, 2012 at 17:47
• 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. Mar 21, 2012 at 18:11
• 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 Mar 21, 2012 at 18:29
• 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. Mar 21, 2012 at 18:36
• 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})$. 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).$$

• It's important to note that the differential condition is necessary, but not sufficient. Mar 21, 2012 at 19:41
• @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 :-) Mar 21, 2012 at 19:55