Wahba's Problem seeks to find the $3\times3$ rotation matrix which minimises:
$$ J(\boldsymbol{\mathrm{R}}) = \frac{1}{2} \sum_{k}|| \boldsymbol{{w}}_k - \boldsymbol{\mathrm{R}}\boldsymbol{v}_k||^2 $$
I have a similar but potentially more difficult problem, where my transformation between two vectors is given by
$$ \boldsymbol{w}_k \approx \boldsymbol{\mathrm{R}}_1\, \boldsymbol{\mathrm{R}}_{2,k}\, \boldsymbol{\mathrm{R}}_3\, \boldsymbol{v}_k $$ where $\boldsymbol{\mathrm{R}}_{2,k}$ is a known Euler rotation about two axes, and $k=1, 2, \cdots, N$.
I wish to determine the best fit rotations $\boldsymbol{\mathrm{R}}_1$, and $\boldsymbol{\mathrm{R}}_3$ which describe the two sets of vectors.
With the exception of using a large set of data in a nonlinear optimisation, is there any way to estimate these two rotation matrices?
Any help is appreciated - Thank you!
$ \def\p{\partial} \def\grad#1#2{\frac{\p #1}{\p #2}} \def\op#1{\operatorname{#1}} \def\LR#1{\left(#1\right)} \def\fracLR#1#2{\LR{\frac{#1}{#2}}} \def\vecc#1{\op{vec}\LR{#1}} \def\diag#1{\op{diag}\LR{#1}} \def\Diag#1{\op{Diag}\LR{#1}} \def\trace#1{\op{Tr}\LR{#1}} \def\qiq{\quad\implies\quad} \def\m#1{\left[\begin{array}{r|r}#1\end{array}\right]} \def\mc#1{\left[\begin{array}{r}#1\end{array}\right]} \def\skew#1{\,\left[\,#1\,\right]_\times} \def\zz{0\;\,} \def\bx{\boxtimes} $The first step is to collect the $v_k\LR{{\rm and}\;w_k}$ vectors into the columns of a matrix
$$\eqalign{ V &= \m{\:v_1 & \:v_2 & \ldots & \:v_n} &\;=\; \sum_k v_ke_k^T \\ W &= \m{w_1 & w_2 & \ldots & w_n} &\;=\; \sum_k w_ke_k^T \\ }$$ where $e_k$ denote the canonical basis vectors.
The Frobenius product can be defined in terms of the trace $${ A:B = \trace{A^TB} }$$ while the Khatri-Rao $(\bx)$ product is the column-wise Kronecker $(\otimes)$ product $$\eqalign{ U &= \m{u_1\,&u_2&\ldots&\,u_n} \\ U &= {V\bx W} \quad\iff\; u_k = {v_k\otimes w_k} \\ }$$ and has the following occasionally useful property $$\eqalign{ \vecc{W\,\Diag{x}\:V^T} = \LR{V\bx W} x \;=\; Ux \\ }$$
Define the matrix variables $$\eqalign{ \def\R{R_b^T\otimes R_a} \def\Rab{R_p} Q_k &= R_a\,R_k\,R_b &\qiq Q_k^TQ_k = I \\ s_k &= \vecc{R_k} &\qiq S = \sum_k s_ke_k^T \\ M &= -2\,US^T \\ \Rab &= \R &\qiq \Rab^T\Rab = I \\ }$$
Use the above ideas to rewrite the cost function $$\eqalign{ \def\J{J_0} J &= \sum_k \LR{w_k - Q_kv_k}^T\LR{w_k - Q_kv_k} \\ &= \sum_k \LR{w_k^Tw_k - 2w_k^TQ_kv_k + v_k^Tv_k} \\ \J &= \LR{J-\sum_k \LR{w_k^Tw_k + v_k^Tv_k}} \\ &= -2\sum_k w_k:Q_kv_k \\ &= -2\sum_k We_k:Q_kVe_k \\ &= -2\sum_k W\LR{e_ke_k^T}V^T:Q_k \\ &= -2\sum_k W\Diag{e_k}V^T : R_aR_kR_b \\ &= -2\sum_k \LR{V\bx W}e_k : \LR{R_b^T\otimes R_a}s_k \\ &= -2\LR{V\bx W}: \LR{R_b^T\otimes R_a}\sum_k s_ke_k^T \\ &= -2U:\Rab\, S \\ &= M:\Rab \\ }$$ Now proceed as per the usual Orthogonal Procrustes Problem, i.e.
compute the SVD of $M$, then $\Rab$ equals the product of the orthogonal factors.
The only twist is that you must decompose $\Rab$ back into $R_a\;{\rm and}\;R_b$
Update
I don't know what your data looks like, however it seems unlikely that $\Rab$ will factor neatly into two orthogonal matrices.
Instead, you can partition $M$ into $3\times 3$ blocks, distribute the Frobenius product, and calculate the differential of the cost function like so $$\eqalign{ \def\Sj{\sum_j} \def\Sk{\sum_k} \def\Mjk{{M}_{jk}} \def\Ejk{E_{jk}} \def\Ekj{E_{kj}} \Ejk &= e_je_k^T \;=\; \Ekj^T \\ M &= \Sj\Sk \Ejk\otimes\Mjk \\ \\ \J &= \LR{\Sj\Sk \Ejk\otimes\Mjk}:\LR{R_b^T\otimes R_a} \\ &= \Sj\Sk\LR{\Ekj:R_b}\LR{\Mjk:R_a} \\ \\ d\J &= \LR{\Sj\Sk\LR{\Ekj:R_b}\Mjk}:dR_a \;+\; \LR{\Sj\Sk\LR{\Mjk:R_a}\Ekj}:dR_b \\ &= G_a:dR_a \;+\; G_b:dR_b \\ }$$ The problem is that $R_a\;{\rm and}\;R_b$ are orthogonal matrices, which constrains the directions of $dR_a\;{\rm and}\;dR_b$ in complicated ways.
However, a rotation matrix can be constructed from an unconstrained vector
via the Rodrigues formula $$\eqalign{ \def\E{\color{red}{\cal E}} \skew{b} &= -\E\cdot b = \mc{ \zz & -b_3 & b_2 \\ b_3 & \zz & -b_1 \\ -b_2 & b_1 & \zz \\ } \\ R_b &= \exp(\skew{b})\\ }$$ where $\E$ is the Levi-Civita tensor.
A paper by Yezzi-Gallego contains an amazing formula for differentiating $R_b$ $$\eqalign{ \def\g{{\large g}} Y_b &= \frac{bb^T + \LR{R_b^T-I}\skew b}{b^Tb} \\ dR_b &= -\LR{R_b\cdot\E\cdot Y_b}\cdot db \\ }$$ Substituting this into the differential of the cost function yields $$\eqalign{ d\J &= &-\LR{G_a:\LR{R_a\cdot\E\cdot Y_a}}\cdot da \\ & &-\LR{G_b:\LR{R_b\cdot\E\cdot Y_b}}\cdot db \\ \grad{\J}{a} &= &\g_a \;=\; -{G_a:\LR{R_a\cdot\E\cdot Y_a}} \\ \grad{\J}{b} &= &\g_b \;=\; -{G_b:\LR{R_b\cdot\E\cdot Y_b}} \\ }$$ These unconstrained vector gradients can be used in an Alternating Least Squares iteration scheme to find optimal $\LR{a,b}$ vectors such that $\,\g_a=\g_b=0,\:$ and the corresponding rotation matrices can be constructed from the optimal vectors.
Note that dot products and double-dot (Frobenius) products extend to tensors $$\eqalign{ \def\A{{\cal A}} \def\B{{\cal B}} \def\C{{\cal G}} \def\c#1{\color{red}{#1}} \def\F{{\cal F}} \F &= \A\cdot\B &\qiq \F_{a\ldots j\,\ell\ldots z} = \Sk\A_{a\ldots j\c{k}}\,\B_{\c{k}\ell\ldots z} \\ \C &= \A:\B &\qiq \C_{a\ldots i\,\ell\ldots z} = \Sj\Sk\A_{a\ldots i\c{jk}}\,\B_{\c{jk}\ell\ldots z} \\ }$$