For the 3D rotation operation $R^{-1}(R(\omega_0)*R(\omega))$, how can we compute the derivative wrt $\omega$? How can we compute the Jacobian derivative of the function:
$$f(\omega) = R^{-1}(R(\omega_0) R(w))$$
with respect to $\omega$, where $\omega_0 \in \mathbb{R}^3$ is some fixed/constant vector, the function:
$$R(\omega) \triangleq \exp \big( [\omega]_{\times} \big) \triangleq \exp 
\left( 
\begin{bmatrix}
 0    & -w_z  &  w_y \\
 w_z  &  0    & -w_x \\
-w_y  &  w_x  &  0 \\
\end{bmatrix}
\right),
$$
is the Rodrigues-vector-to-rotation mapping, and $R^{-1}(\cdot)$ is the corresponding inverse function?
For the sake of this question, we are primarily interested in taking the derivative about the point $\omega=0$, since other values of $\omega$ can be absorbed into the constant $\omega_0$ with a little extra effort. So everything should be expressible in terms of $\omega_0$.

Edit: Using some dirty empirical methods, I was able to get the first few coefficients of the Taylor expansion as:
$$
\frac{d f}{d \omega} = Z^0 + \frac{1}{2!} Z^1+ \frac{1}{2 (3!)}Z^2 - \frac{1}{6 (5!)}Z^4 + \frac{1}{6 (7!)} Z^6 - (...)
$$
where $Z=[\omega_0]_{\times}$. (Evidently, there no higher order odd terms beyond the first.) So this appears to have a similar form to Rodrigues, as:
$$\frac{d f}{d \omega} = I + \frac{1}{2} Z + \beta(\omega_0) Z^2$$
but getting a functional form for $\beta$ seems nontrivial. (I'd venture a guess than $\beta$ is purely a function of $\|\omega_0\|$.)
 A: $
\def\a{\alpha}
\def\b{\beta}
\def\g{\gamma}
\def\m{\lambda}
\def\e{\varepsilon}
\def\t{\theta}
\def\s{\sigma}
\def\w{\omega}
\def\o{{\tt1}}
\def\c#1{\color{red}{#1}}
\def\l{\left}\def
\r{\right}\def\lr#1{\l(#1\r)}
\def\LR#1{\bigg(#1\bigg)}
\def\fracLR#1#2{\lr{\frac{#1}{#2}}}
\def\p{\partial}
\def\grad#1#2{\frac{\p #1}{\p #2}}
\def\gradLR#1#2{\fracLR{\p #1}{\p #2}}
\def\cay#1{\operatorname{cay}\lr{#1}}
\def\skew#1{\operatorname{skew}\lr{#1}}
\def\pp#1{\operatorname{psi}\lr{#1}}
\def\ndot{\mathop{\bullet}\limits}
$First, a sketch of the Rodrigues and Gallego & Yezzi formulas.
Let $\e$ denote the Levi-Civita tensor, then given a vector $(a),\,$ its associated magnitude $(\a)$, skew matrix $(A)$, and rotation matrix $(R)$ can be calculated as
$$\eqalign{
\a &= \|a\|,\;\qquad A = -\e\cdot a \\
R &= e^A = I + \rho_1 A + \rho_2 A^2 \\
\rho_1 &= \frac{\sin(\a)}{\a} \quad\quad \rho_2=\frac{\o-\cos(\a)}{\a^2} \\
}$$
The skew matrix has an interesting periodic property for any $n\ge\o$
$$A^{n+2} = -\a^2A^n$$
allowing any analytic function of $A$ to be reduced to a quadratic polynomial.

Gallego & Yezzi utilize one such function
$$\eqalign{
G &= \pp{A} = I + \g_1 A + \g_2 A^2 \\
\g_1 &= \frac{\cos(\a)-\o}{\a^2} = -\rho_2 \\
\g_2 &= \frac{\a-\sin(\a)}{\a^3} = \frac{\o-\rho_1}{\a^2}  \\
}$$
where $\phi(x) = \fracLR{e^x-1}{x}\;$ is a function commonly seen in exponential integrators
$\;$ and $\;\pp{x} = \phi(x)\,e^{-x}$
They actually employ a different (but equivalent) expression for $G$
$$G = \fracLR{aa^T+(R^T-I)A}{\a^2}$$
and derive an amazing formula for the gradient of $Ru$
(where $u$ is an arbitrary vector)
$$\eqalign{
\grad{(Ru)}{a} &= R\cdot\lr{\e\cdot u}\cdot G \\
}$$
Multiplying each side by $e_k\,$ yields their so-called
compact formula for the gradient of $R$
$$\eqalign{
&\grad{(Ru)}{a_k} &= +R\cdot\lr{\e\cdot u}\cdot(Ge_k) \\
&\grad{R}{a_k}\cdot u &= -R\cdot\lr{\e\cdot g_k}\cdot u \\
&\grad{R}{a_k} &= -R\cdot\e\cdot g_k \\
&\grad{R}{a} &= -R\cdot\e\cdot G \\
\\
}$$

In the current problem, instead of a single vector $\{a\}$, there are three vectors $\{f,w_0,w\}$ whose respective rotation matrices satisfy the relationship
$$R_f = R_0\,R_w $$
The differential of this relationship leads directly to
the desired gradient
$$\eqalign{
dR_f &= R_0\cdot dR_w \\
\gradLR{R_f}{f}\cdot df &= R_0\cdot\gradLR{R_w}{w}\cdot dw \\
\c{R_f\cdot\e}\cdot G_f\cdot df &= \c{R_0\cdot R_w\cdot\e}\cdot G_w\cdot dw \\
G_f\cdot df &= G_w\cdot dw \\
\grad{f}{w} &= G_f^{-1}G_w \\\\
}$$

Update
There was a question in one of the comments on how to invert $G$.
Since the inverse function will also be a quadratic polynomial
one can immediately write
$$\eqalign{
B &= G^{-1} = I + \b_1 A + \b_2 A^2 \\ 
}$$
The equation $BG=I$ generates a $2\times 2$ system for the
$\beta$-coefficients in terms of the $\g$-coefficients, which can be solved using Cramer's Rule. Then the $\{\g_k\}$ are substituted by the
$\{\rho_k\}$ to obtain the result in terms of the Rodrigues coefficients.
$$\eqalign{
\b_1 &= \frac{1}{2}, \quad
\b_2 = \frac{2\rho_2-\rho_1}{2\a^2\rho_2} \\
}$$
The various coefficients are trigonometric functions of $\a$ which have finite limits as $\a\to 0$
$$\eqalign{
\rho_1 &= 1, \quad &\rho_2 &= \frac{1}{2}, \quad
\g_1 &= \frac{-1}{2}, \quad &\g_2 &= \frac{1}{6}, \quad
\b_1 &= \frac{1}{2}, \quad &\b_2 &= \frac{1}{12} \\
}$$
A: To begin with, let's get rid of the logarithm map, and hit both sides by an arbitrary test vector $\mathbf{v}$:
$$R_{f(\omega)}\mathbf{v} = R_{\omega_0} R_{\omega}\mathbf{v}.$$
Now we can differentiate both sides, making use of the formulas (derived using geometric arguments) in Gallego and Yezzi, A Compact Formula for the Derivative of a 3-D Rotation in Exponential Coordinates, Journal of Mathematical Imaging and Vision, Volume 51 Issue 3, March 2015, 378--384:
$$\left[\left(d R_{\omega}\right) \delta \omega\right]\mathbf{v} = -R_{\omega} [\mathbf{v}]_\times T(\omega)\delta \omega,$$
for variation $\delta\omega$ and 
$$T(\omega) = \begin{cases}I, & \|\omega\| = 0\\
\frac{\omega\omega^T + (R_{-\omega}-I)[\omega]_\times}{\|\omega\|^2}, & \|\omega\| > 0.\end{cases}$$
We have
$$\left[dR_{f(\omega)} df(\omega)\delta\omega\right] \mathbf{v} = -R_{\omega_0} R_{\omega} [\mathbf{v}]_{\times} T(\omega) \delta \omega$$
$$-R_{f(\omega)} \left[\mathbf{v}\right]_{\times} T[f(\omega)] df(\omega) \delta \omega =-R_{\omega_0} R_{\omega} [\mathbf{v}]_{\times} T(\omega) \delta \omega.$$
The rotations at the beginning cancel, and since the variation $\delta\omega$ and test vector $\mathbf{v}$ are arbitrary, we must have equality of matrices
$$df(\omega) = T[f(\omega)]^{-1}T(\omega),$$
with special case
$$df(0) = T(\omega_0)^{-1}.$$
