Numerical derivative of a matrix depending on a Rotation matrix I want to derivate a matrix $A$ w.r.t. a rotation matrix numerically. 
The rotation matrix I am talking about is $R_i \in SO(3)$ and it is coming from a choice of Euler angles sequence $(1,2,3)$. The angles involved are the so called roll-pitch-yaw angles $\phi_i,\theta_i,\psi_i$  (or Cardan Angles, or Tait-Bryan angles).
Therefore:
$R_i(\phi_i,\theta_i,\psi_i)=R_i=R_x(\phi_i)R_y(\theta_i)R_z(\psi_i) \;\;\;\;\;\;(1)$
If I would need to derivate the matrix $A$ numerically w.r.t. a scalar $\alpha$ I would do: 
$\frac{A(\alpha)-A(\alpha+\delta)}{\delta}\;\;\;\;\;\;(2)$
Instead in my case I was thinking to do, since I can write a rotation matrix as $R(\omega,\beta)=exp(\hat{\omega}\beta)$ (where $\hat{\omega}$ is the skew-symmetric matrix associated to the vector $\omega \in \mathbb{R}^3$): 
$\frac{A(R_i)-A(R_i\cdot exp(e_i\delta))}{\delta}\;\;\;\;\;\; \forall i \in\{1,2,3\}\;\;\;\;\;\;(3)$
with 
$
e_1=
\left(
\begin{array}{c}
1\\0\\0
\end{array}
\right)
\;\;
,
e_2=
\left(
\begin{array}{c}
0\\1\\0
\end{array}
\right)
,
e_3=
\left(
\begin{array}{c}
0\\0\\1
\end{array}
\right)
$
I am not sure the $(3)$ is correct and I am not sure its denominator is correct. 
Could you please help me? 
Thanks a lot.
 A: If you have a sequence of rotations about the local axes $\mathbf{z}_1$, $\mathbf{z}_2$ and $\mathbf{z}_3$ with angles $q_1$, $q_2$ and $q_3$ then you have the following properties
$$ \begin{align} \mathtt{R} &= \mathrm{Rot}(\mathbf{z}_1,q_1)\mathrm{Rot}(\mathbf{z}_2,q_2) \mathrm{Rot}(\mathbf{z}_3,q_3) \\
\dot{\mathtt{R}} & = [\boldsymbol{\omega} \times] \mathtt{R} \\
\boldsymbol{\omega} & = \mathbf{z}_1 \dot{q}_1 + \mathrm{Rot}(\mathbf{z}_1,q_1) \left( \mathbf{z}_2 \dot{q}_2 + \mathrm{Rot}(\mathbf{z}_2,q_2) \mathbf{z}_3 \dot{q}_3 \right)
\end{align} $$
NOTE: The $[\omega \times]$ notation is the 3x3 anti-symmetric cross product operator matrix. $$ [ \begin{pmatrix} x\\y\\z \end{pmatrix} \times ] = \begin{bmatrix} 0 & -z & y \\ z & 0 & -x \\ -y& x & 0 \end{bmatrix} $$ 
You can easily prove the above with two rotations if you accept the rules on how to differentiate vectors that ride on rotating frames: (search rotating frame derivative). Expanding it to three rotations is more tedious, but following the same logic. $$ \dot{\mathtt{R}} = \dot{\rm R}_1 {\rm R}_2 {\rm R}_3 +{\rm R}_1 \dot{\rm R}_2 {\rm R}_3+{\rm R}_1 {\rm R}_2 \dot{\rm R}_3 $$
A: The gradient you seek is a 4th order tensor$,\,{\mathcal F},\,$ which satisfies
$$\eqalign{
dA_{ij} &= {\mathcal F}_{ijkl}\,dR_{kl}
\cr\cr
}$$
Since $dR$ has 9 independent elements, you'll to calculate need 9 numerical derivatives
$$\eqalign{
\frac{A(R+hE_{kl})-A(R)}{h}
}$$
where $E_{kl}$ is the matrix which has the $(k,l)$ element equal to unity, and all other elements equal to zero. The $h$ parameter should be chosen to balance round-off error against accuracy. For IEEE double precision its around $10^{-8}$
Update
Here's an explicit version of the proposed technique 
$$\eqalign{
{\mathcal F}_{ijkl} &= \frac{A_{ij}(R+hE_{kl})-A_{ij}(R)}{h}
}$$
