Is there an elegant way to solve $b = R\left(\theta \hat{u} \right) a$ where $a, b \in \mathbb{R}^3; R \in SO(3)$ for $\hat{u}$

Question

I have a problem where:

$ b = R \left(\theta \hat{u} \right) a $

where $a, b \in \mathbb{R}^3$ are known, and $R\left(\theta \hat{u} \right) \in SO(3)$ is a rotation matrix constructed from angle $\theta$ about axis $\hat{u}$ of unit length. In this particular problem, $\theta$ is known and I wish to solve for $\hat{u}$.

I can brute force this from the Rodrigues rotation formula and solve for the elements of $\hat{u}$ (which turns out to be a three quadric intersection problem) but is there a more elegant way to solve this problem?

Clarifying notes:

Angle $\theta$ is fixed in this problem and is key. I am speculating that this should reduce the family of rotations that would otherwise satisfy the first equation to a unique (or at least finite) number of solutions.

---

EDIT: Counter-example to the cross-product answer.

Consider the rotation matrix constructed as follows:

uhat = [4; 6; 8] / norm([4; 6; 8])

uhat =

   0.37139
   0.55709
   0.74278

theta = 0.5;
R = RodriguesRotation(uhat, theta)

R =

   0.89447  -0.33078   0.30085
   0.38144   0.91557  -0.12740
  -0.23331   0.22871   0.94512

With vectors $a$ and $b$ being:

a = [1; 2; 3]
b = R * a

b =

   1.1355
   1.8304
   3.0595

Then computing $\hat{u} = \frac{a \times b}{|a\times b|}$:

axb = cross(a, b);
unit_axb = axb / norm(axb)

unit_axb =

   0.74582
   0.41213
  -0.52336

Which is clearly does not recover the original rotation axis.

---

Edit 2: Widawensen's solution can be verified:

unit_a = a / norm(a);
unit_b = b / norm(b);
cos_theta = cos(theta);

cos_beta = sqrt((cos_alpha - 1) / (cos_theta - 1));
beta = acos(cos_beta_plus);
gamma = pi/2 - beta;

Verifying the dot products:

dot(unit_a, uhat) - cos(gamma)
ans =   -8.8818e-16

dot(unit_b, uhat) - cos(gamma)
ans =   -8.8818e-16

Answer 1

Without loss of generality we can assume that given vectors $a,b$ are of unit length (in other case we can normalize them to unit length). Also axis vector $u$ is of unit length. This gives more compact way of describing the dot product for these vector as they can be expressed as cosines of angles between unit vectors.

As we know skew-symetric matrix $S(u)$, made from components of vector $u$, is a second-rank matrix so the image of the space in transformation $S(u)$ is a plane and $u$ is perpendicular to this plane.

We know in the problem the angle of the rotation $\theta$, what means that the angle between orthogonal projections of $a$ and $b$ on plane $S(u)$ is just $\theta$ (rotation is operation acting on vectors lying in the plane $S(u)$, perpendicular to the axis of rotation).

If we denote as $\beta$ the angle between vector $a$ and its projection $a_S=(I-uu^T)a $ on plane $S(u)$ then the length for projected vector $a_S$ is $\cos\beta$.
The same can be said about projection $b_S $ of vector $b$.
Notice that at the same time angle between $a$ (also $b$) and $u$ equals to $90-\beta$.

We can write down the appropriate equation for dot product of $a_S$ and $b_S$ in two ways (left hand side and right hand side of the equation below, transformed below is LHS)

$((I-uu^T)a)^T(I-uu^T)b=\cos\theta(\cos\beta)^2 $

$(a-uu^Ta)^T(b-uu^Tb) =\cos\theta(\cos\beta)^2$

$(a^T-a^Tuu^T)(b-uu^Tb)=\cos\theta(\cos\beta)^2$

$a^Tb-(a^Tu)(u^Tb)=\cos\theta(\cos\beta)^2$

$a\circ b - (a\circ u)( b \circ u)=cos\theta(\cos\beta)^2$

Denote the angle between $a$ and $b$ as $\alpha$.

Then we have

$\cos\alpha - \cos(90-\beta)\cos(90-\beta)= \cos\theta(\cos\beta)^2$

$\cos\alpha - (\sin\beta)^2 = \cos\theta(\cos\beta)^2$

$ \bbox[yellow,5px,border:2px solid red] { \cos\alpha = (\sin\beta)^2+ \cos\theta(\cos\beta)^2} $

This nice looking equation allows easy solution for $(\cos\beta)^2$ and then for $\cos\beta$ and $\sin\beta$

We can denote $90-\beta =\gamma$, where $\gamma$ is the angle between vector $a$ and $u$.

This value of angle $\gamma$ is the same for vectors $b$ and $u$.

Knowing the angle $\gamma$ you have two equations from dot product:

$u \circ a= \cos\gamma$ and $u \circ b= \cos\gamma$.

Additionally, these dot products are equal to the same value ($\cos\gamma$) what means that $a-b$ is perpendicular to the axis vector $u$.

Here we have really two independent unknowns in $u$ as the third one is linked with two others by unit length of vector $u$ (the third equation).

Three mentioned equations lead to the desired solution.

Answer 2

The axis of rotation $\hat{u}$ for a rotation that sends $a$ to $b$, lies in the plane of symmetry between $a$ and $b$. This plane is spanned by two vectors $v_1$ and $v_2$ where

$\hat{v}_1 = \dfrac{a+b}{\| a + b \|} $

$ \hat{v}_2 = \dfrac{ a \times b }{\| a \times b \| } $

Note that $v_1$ and $v_2$ are unit vectors and are orthogonal to each other.

Therefore, the unit vector $\hat{u}$ can be expressed as follows:

$ \hat{u} = \cos \psi \ \hat{v}_1 + \sin \psi \ \hat{v}_2 $

Where $\psi$ is to be found. It is straight forward to see that if the rotation angle $\theta$ is positive, then we must have $ \sin \psi \ge 0 $ and vice versa.

Now from the Rodrigues rotation formula, we have

$ b = { \hat{u} \hat{u} }^T a + \cos \theta (I - { \hat{u} \hat{u} }^T ) a + \sin \theta S_{\hat{u}} a \tag{1}$

We now have that

$ {\hat{u} \hat{u} }^T = \cos^2 \psi \ { \hat{v}_1 \hat{v}_1 }^T + \sin^2 \psi \ { \hat{v}_2 \hat{v}_2 }^T + \cos \psi \sin \psi \ ( { \hat{v}_1 \hat{v}_2}^T + { \hat{v}_2 \hat{v}_1}^T ) $

And

$ S_{\hat{u}} = \cos \psi S_{v_1} + \sin \psi S_{v_2} $

(Remember the matrix $S_w a = w \times a $ )

Dotting both sides of equation $(1)$ with the vector $a$ gives

$ a \cdot b = \cos \theta (a \cdot a) + (1 - \cos \theta) \bigg( \cos^2 \psi (a \cdot \hat{v}_1)^2 + \sin^2 \psi \ (a \cdot \hat{v}_2)^2 + 2 \cos \psi \sin \psi ( a \cdot \hat{v}_1 )(a \cdot {v}_2 ) \bigg)$

which is a simple trigonometric equation in double the angle $\psi$. It can be simplified to

$ a \cdot b = \cos \theta (a \cdot a) + \dfrac{1}{2}(1 - \cos \theta ) \bigg( (a \cdot \hat{v}_1)^2 + (a \cdot \hat{v}_2)^2 + \bigg( (a \cdot \hat{v}_1)^2 - (a \cdot \hat{v}_2)^2 \bigg) \cos(2 \psi) \\ + 2 \sin(2 \psi) (a \cdot \hat{v}_1)(a \cdot \hat{v}_2 ) \bigg) $

From here, it is straight forward to find the four possible solutions for this trigonometric equation, then choose the two valid solutions, which are the one the same sign of $\sin \psi$ as $\theta$.

Once the two possible values of $\psi$ are determined, then the two possible axes of rotation are determined from the formula above.

To verify the validity of the above procedure, I let vector $a$ be

$ a = \begin{bmatrix} 1 \\ 2 \\ 3 \end{bmatrix}$

To create a rotated vector $b$, I let the rotation axis be along $N$

$ N = \begin{bmatrix} 2 \\ 3 \\ 4 \end{bmatrix} $

So that the unit vector along the axis of rotation is

$ u = \dfrac{N}{\| N \|} = \begin{bmatrix} 0.371391 \\ 0.557086 \\ 0.742781 \end{bmatrix}$

From here, and using an angle of rotation equal to $\theta = 0.5 \text{radian} $, and using the Rodrigues rotation formula, I found that the rotated vector is

$ b = \begin{bmatrix} 1.135461 \\ 1.830388 \\ 3.059478 \end{bmatrix} $

That is all the necessary input information.

Using the above procedure, I solve for $\psi$ and chose the valid values (the ones that have $\sin(\psi) \gt 0 $ because $\theta \gt 0 $).

The valid solutions are

$ \psi_1 = 3.023478 $ and $ \psi_2 = 0.118115 $

Corresponding to $\psi_1$ is the axis

$ \hat{u} = \begin{bmatrix} -0.19562 \\ -0.45996 \\ -0.86613 \end{bmatrix} $

And corresponding to $\psi_2$ is the axis

$ \hat{u} = \begin{bmatrix} 0.371391 \\ 0.557086 \\ 0.742781 \end{bmatrix}$

Which is the axis we set up, along with the rotation angle $\theta$, to find vector $b$.

I then verified that indeed the other axis gives the same rotated vector $b$, with the given angle of rotation.

This validates the above procedure.

Answer 3

$ \def\lR#1{\Big(#1\Big)} \def\LR#1{\left(#1\right)} \def\op{\operatorname} \def\q{\quad} \def\qq{\qquad} \def\c{\gamma} \def\s{\sigma} \def\t{\theta} \def\l{\lambda} $The problem can be approached from a purely numerical perspective, without requiring any kind of deep geometric insight.

Use the known $\t$ to define the scalar constants $$\eqalign{ \c &= \cos\t,\qq \s=\sin\t,\qq \tau=(1-\c) \\ }$$ Although $u$ is unknown, we can use it as a parameter to define the matrices $$\eqalign{ U &= \op{skew}(u)=(u\times I),\qq R=\c I+\s U + \tau uu^T \\ }$$ The equation which must be satisfied is $$\eqalign{ b &= Ra \\ &= \c a +\s Ua +\tau(u^Ta)u \\ &= \c a +\s(u\times a) +\tau(u\cdot a)u \\ }$$ Rearrange this to isolate $u$ on the LHS $$\eqalign{ u &= \frac{b-\c a -\s(u\times a)}{\tau(u\cdot a)} \\ }$$ While this is not a solution, it does suggest an iterative scheme: $$\eqalign{ u \;=& random\ \ unit\ \ vector \\ {\sf do}\q \\ &c = \frac{b-\c a -\s(u\times a)}{\tau(u\cdot a)} \\ &u = u + \tfrac18c \\ &u = \frac{u}{\|u\|} \\ {\sf until} &{\sf converged} \\ }$$ A dampening factor of $\frac18$ (or smaller) seems to be required for convergence. If you need to solve this sort of problem often, you might want to invest some time in finding the optimal dampening factor.

Update

Here is an iteration which converges much more rapidly.

Define $$\eqalign{ s &= (b-\c a),\qq A=\op{skew}(a) \\ }$$ isolate $u$ as before $$\eqalign{ u &= \frac{b-\c a -\s(u\times a)}{\tau(u\cdot a)} \;\equiv\; \frac{s+\s Au}{\tau a\cdot u} \\ }$$ then rearrange it into a matrix equation $$\eqalign{ \lR{\LR{\tau a\cdot u}I - \s A}\,u = s \\ }$$ The following scheme converges without any dampening factors $$\eqalign{ u \;=& random\ \ unit\ \ vector \\ {\sf do}\q \\ &u = \lR{\LR{\tau a\cdot u}I - \s A}^{-1}s \\ &u = \frac{u}{\|u\|} \\ {\sf until} &{\sf converged} \\ }$$ $\sf NB\!:$ Random starting points converge to different solutions (see Hosam's excellent post)