An easily-to-implement solution is to formulate a linear equation system for the coefficients of $\mathbf{R}_{ab}$ (solvable in MatLab, Python, ...), although these nine coefficients are not independent of each other (constraints of SO3-group), anyway it works fairly well
\begin{equation}
\left[
\begin{array}{cccccc}
f_{ax}& f_{ay} & f_{az} & 0 & 0 &0 & 0 & 0 &0 \\
0 & 0 &0 & f_{ax}& f_{ay} & f_{az} & 0 & 0 &0 \\
0 & 0 &0 & 0 & 0 &0 & f_{ax}& f_{ay} & f_{az}\\
g_{ax}& g_{ay} & g_{az} & 0 & 0 &0 & 0 & 0 &0 \\
0 & 0 &0 & g_{ax}& g_{ay} & g_{az} & 0 & 0 &0 \\
0 & 0 &0 & 0 & 0 &0 & g_{ax}& g_{ay} & g_{az}\\
h_{ax}& h_{ay} & h_{az} & 0 & 0 &0 & 0 & 0 &0 \\
0 & 0 &0 & h_{ax}& h_{ay} & h_{az} & 0 & 0 &0 \\
0 & 0 &0 & 0 & 0 &0 & h_{ax}& h_{ay} & h_{az}
\end{array}
\right]
\left[
\begin{array}{c}
R_{xx}\\ R_{xy}\\ R_{xz}\\ R_{yx}\\ R_{yy}\\ R_{yz}\\ R_{zx}\\ R_{zy}\\ R_{zz}
\end{array}
\right]
=
\left[
\begin{array}{c}
f_{bx}\\ f_{by}\\ f_{bz}\\
g_{bx}\\ g_{by}\\ g_{bz}\\
h_{bx}\\ h_{by}\\ h_{bz}
\end{array}
\right]
\end{equation}
Alternatively, you may choose parametrizations with less parameters, and consequently solve less, but nonlinear equations (choose the convenient ones in the problem at hand). For example with Euler-parameters ($a$, $b$, $c$ $d$), you would express all components of $\mathbf{R}(a,b,c,d)$ in dependency on four parameters: $R_{xx}=a^2+b^2-c^2-d^2$, $\ $ $R_{xy}=2(bc-ad)$, $\ $ $\dots$
If you were free to choose orthonormal basis vectors in one of the two systems (A or B), then you would directly get the rows in one direction (other direction is the transpose $\mathbf{R}_{ab}^\text{T} = \mathbf{R}_{ba}$, since rotation matrices $\mathbf{R}$ are orthogonal).
scipy.spatial.transformmay provide many operations you need. – Dominik Kern Oct 17 '22 at 14:51