# 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}$

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


• Hello Damien. my method is good for you? Elegant? Commented Apr 21, 2020 at 16:03
• It is indeed. I will try to verify it with the above example before accepting. Commented Apr 21, 2020 at 16:37
• Of course, verification is very important. I hope I didn't made an error in calculations. Commented Apr 21, 2020 at 16:41

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.

• Verified as being a correct solution. Thank you. Commented Apr 21, 2020 at 16:46
• @Damien Thank you for the verification. Thank you for the posting this interesting problem which extended also my knowledge about rotations. Commented Apr 21, 2020 at 17:47
• Interesting is to make a graph ( for example in desmos.com) for $\cos\ y\ =\left(\sin\left(x\right)\right)^{2}+\cos\left(\ t\right)\ \left(\cos\left(x\right)\right)^{2}$ where t is a parameter for the angle of rotation . Such graph shows for a fixed $t$ a relationship between the angle of vector with the rotation plane and obtained angle between vector and its image. Interesting here is , for example , transition of the graph form near value $t=\pi$. Commented Apr 22, 2020 at 5:35
• Second comment: Since the highlighted equation is a square, then there is both + and - cos_gamma solutions and thus there are two axes of rotations for a given a, b and theta Commented Apr 23, 2020 at 22:49
• Third comment: The simultaneous equations $u \dot a = \cos \gamma$ and $u \dot b = \cos \gamma$ is solvable via an under-determined set of linear equations, subject to a unit norm constraint. Easiest via an SVD to get the minimum norm solution (which if > 1 means there is no solution), then using the nullspace (right vector of v) and finding the coefficient that makes the norm(x) == 1. See math.stackexchange.com/questions/2253443/…. There will usually be two solutions that fit the unit norm constraint, but only one will work. Commented Apr 23, 2020 at 22:53