Assuming we already have a sequence of rotations $\{R_i\}_{i=1}^N$ that when chained together is close to identity rotation E: $$ \prod_{i=1}^NR_i\approx E, $$ then the problem is defined to be finding a small rotation perturbation $\Delta R$ which will make the above equation exactly identity: $$ \prod_{i=1}^NR_i\Delta R=E. $$
The following notations are used for the convenience of expression: $$ R_{i,j}:=\prod_{k=i}^jR_i,\\ [\mathbf {w} ]^{\wedge}:={\begin{bmatrix}\,\,0&\!-w_{3}&\,\,\,w_{2}\\\,\,\,w_{3}&0&\!-w_{1}\\\!-w_{2}&\,\,w_{1}&\,\,0\end{bmatrix}},\\ \{[\mathbf {w}]^{\wedge}\}^{\vee}:=\mathbf{w},\\ R^{\alpha}:=\exp(\alpha\ln(R)). $$
In the 2D case the solution can be easiy expressed as $\Delta R=R_{1,N}^{-\frac{1}{N}}$ because rotation composition is commutative. But in the 3D case obtaining the exact solution involves solving a complex transcendental equation. In practice I only need a simple iterative method which promises to converge to the real solution.
I tried the heuristic solution as $$ \Delta R \gets R_{1,N}^{-\frac{1}{N}},\\ R_{i}\gets R_i\Delta R, $$ and repeat the above iteration. The basic idea is that rotation composition is nearly commutative near identity $E$. However, this method seems to converge to a strange rotation that still has some gap between $E$, i.e., $R_{1,N}\to R_0 \neq E$.
I then derived another solution as follows: Let $\exp([\mathbf{w}]^\wedge)= \Delta R$, then assuming $\Delta R$ and $R_{1,N}$ are close to $E$, we have $\Delta R\approx E+[\mathbf{w}]^{\wedge}+O(||\mathbf{w}||^2)$ and $R_{1,N}\approx E+[\mathbf{w}_c]^\wedge$. Then the equation can be expressed as: \begin{align} E&=\prod_{i=1}^NR_i(E+[\mathbf{w}]^\wedge+O(||\mathbf{w}||^2)\\ &=R_{1,N}+\sum_{i=1}^NR_{1,i}[\mathbf{w}]^\wedge R_{i+1,N}+O(||\mathbf{w}||^2) \end{align} Right multipliing $R_{1,N}^T\approx E-[\mathbf{w}_c]^\wedge$ on both side and ignoring the second order infinitesimal, we have: \begin{align} E-[\mathbf{w}_c]^\wedge=E+\sum_{i=1}^NR_{1,i}[\mathbf{w}]^\wedge R_{1,i}^T. \end{align} Because $^\vee$ is a linear operator on skew-symmetric matrices and $\{R[\mathbf{w}]^\wedge R^T\}^\vee=R\mathbf{w}$, we obtain the following equation: $$ (\sum_{i=1}^NR_{1,i})\mathbf{w}=-\mathbf{w}_c. $$ Similarily if $\Delta R$ is multiplied on the left ($\Delta RR_i$), we have: $$ (\sum_{i=1}^NR_{i,N})^T\mathbf{w}=-\mathbf{w}_c. $$ This method can be viewed as mathcing the first-order gradient of the rotation vector and can be implemented iteratively: $$ \Delta R\gets\exp(\{-(\sum_{i=1}^NR_{1,i})^{-1}[\ln(R_{1,N})]^\vee\}^\wedge),\\ R_i\gets R_i\Delta R. $$ However, this solution seems to be oscilliting and I have no idea why it would behave like this.
Anyone help is greatly appreciated!
Note: Different kinds of perturbation is also allowed like a small perturbation in the Lie algebra or in the rotation exponent if it can result in simple iterative method when solving this rotation perturbation problem.