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


 A: 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.
