Let's say we have two sets of vectors $A$, $B$ and $A'$, $B'$. These vectors are connected such that
$A' = R. A$, and
$B'=R.B$. Where R is a $3D$ rotation matrix.
Now let's assume we know $A$, $B$ and $A'$, $B'$. Then is there any way to know $R$? If yes, then what should be the approach to find it.
On
First, let's assume $$ \vec{A} = \left [ \begin{matrix} A_x \\ A_y \\ A_z \end{matrix} \right ] \quad \text{and} \quad \vec{B} = \left [ \begin{matrix} B_x \\ B_y \\ B_z \end{matrix} \right ] $$ are nonzero ($\lVert\vec{A}\rVert\ne 0$, $\lVert\vec{B}\rVert\ne 0$) and linearly independent ($\lVert\vec{A}\cdot\vec{B}\rVert \lt \lVert\vec{A}\rVert \lVert\vec{B}\rVert$, or in other words, there does not exist a real $\lambda$ for which $\lVert\vec{A} - \lambda \vec{B}\rVert = 0$).
Let's also assume that there exists a rotation matrix $\mathbf{B}$, $$\mathbf{R} = \left [ \begin{matrix} R_{11} & R_{12} & R_{13} \\ R_{21} & R_{22} & R_{23} \\ R_{31} & R_{32} & R_{33} \end{matrix} \right ]$$ so that $$\mathbf{R}\vec{A} = \vec{a} = \left [ \begin{matrix} a_x \\ a_y \\ a_z \end{matrix} \right ] \quad \text{and} \quad \mathbf{R}\vec{B} = \vec{b} = \left [ \begin{matrix} b_x \\ b_y \\ b_z \end{matrix} \right ]$$
If $\mathbf{R}$ is a rotation matrix, then also $$\mathbf{R}(\vec{A} \times \vec{B}) = \vec{a} \times \vec{b} = \left [ \begin{matrix} a_y b_z - a_z b_y \\ a_z b_x - a_x b_z \\ a_x b_y - a_y b_x \end{matrix} \right ]$$ These three matrix-vector products give you a system of nine equations: $$\left\lbrace\begin{aligned} R_{11} A_x + R_{12} A_y + R_{13} A_z &= a_x \\ R_{21} A_x + R_{22} A_y + R_{23} A_z &= a_y \\ R_{31} A_x + R_{32} A_y + R_{33} A_z &= a_z \\ R_{11} B_x + R_{12} B_y + R_{13} B_z &= b_x \\ R_{21} B_x + R_{22} B_y + R_{23} B_z &= b_y \\ R_{31} B_x + R_{32} B_y + R_{33} B_z &= b_z \\ R_{11} (A_y B_z - A_z B_y) + R_{12} (A_z B_x - A_x B_z) + R_{13} (A_x B_y - A_y B_x ) &= a_y b_z - a_z b_y \\ R_{21} (A_y B_z - A_z B_y) + R_{22} (A_z B_x - A_x B_z) + R_{23} (A_x B_y - A_y B_x ) &= a_z b_x - a_x b_z \\ R_{31} (A_y B_z - A_z B_y) + R_{32} (A_z B_x - A_x B_z) + R_{33} (A_x B_y - A_y B_x ) &= a_x b_y - a_y b_x \\ \end{aligned}\right.$$ which you can solve for $R_{11}, \dots, R_{33}$.
I checked the solution using Maple. There are quite a few terms, so I'll just tabulate them here. All $R_{11}, \dots, R_{33}$ have a common divisor, d = $(\vec{A} \times \vec{B})\cdot(\vec{A} \times \vec{B})$, where $\times$ is vector cross product operation.
d = (A_x * B_y - A_y * B_x)*(A_x * B_y - A_y * B_x)
+ (A_x * B_z - A_z * B_x)*(A_x * B_z - A_z * B_x)
+ (A_y * B_z - A_z * B_y)*(A_y * B_z - A_z * B_y)
R_11 = ( A_y * A_y * B_x * b_x + A_z * A_z * B_x * b_x
- A_x * A_y * B_y * b_x - A_x * A_z * B_z * b_x
+ A_x * a_x * B_y * B_y + A_x * a_x * B_z * B_z
- A_y * a_x * B_x * B_y - A_z * a_x * B_x * B_z
+ A_y * a_y * B_z * b_z - A_z * a_y * B_y * b_z
- A_y * a_z * B_z * b_y + A_z * a_z * B_y * b_y ) / d
R_12 = ( A_x * A_x * B_y * b_x - A_x * A_y * B_x * b_x
- A_x * B_x * B_y * a_x - A_x * B_z * a_y * b_z
+ A_x * B_z * a_z * b_y - A_y * A_z * B_z * b_x
+ A_y * B_x * B_x * a_x + A_y * B_z * B_z * a_x
+ A_z * A_z * B_y * b_x + A_z * B_x * a_y * b_z
- A_z * B_x * a_z * b_y - A_z * B_y * B_z * a_x) / d
R_13 = ( A_x * A_x * B_z * b_x - A_x * A_z * B_x * b_x
- A_x * B_x * B_z * a_x + A_x * B_y * a_y * b_z
- A_x * B_y * a_z * b_y + A_y * A_y * B_z * b_x
- A_y * A_z * B_y * b_x - A_y * B_x * a_y * b_z
+ A_y * B_x * a_z * b_y - A_y * B_y * B_z * a_x
+ A_z * B_x * B_x * a_x + A_z * B_y * B_y * a_x ) / d
R_21 = ( A_x * B_y * B_y * a_y - A_x * A_y * B_y * b_y
- A_x * A_z * B_z * b_y + A_x * B_z * B_z * a_y
+ A_y * A_y * B_x * b_y - A_y * B_x * B_y * a_y
- A_y * B_z * a_x * b_z + A_y * B_z * a_z * b_x
+ A_z * A_z * B_x * b_y - A_z * B_x * B_z * a_y
+ A_z * B_y * a_x * b_z - A_z * B_y * a_z * b_x ) / d
R_22 = ( A_x * A_x * B_y * b_y - A_x * A_y * B_x * b_y
- A_x * B_x * B_y * a_y + A_x * B_z * a_x * b_z
- A_x * B_z * a_z * b_x - A_y * A_z * B_z * b_y
+ A_y * B_x * B_x * a_y + A_y * B_z * B_z * a_y
+ A_z * A_z * B_y * b_y - A_z * B_x * a_x * b_z
+ A_z * B_x * a_z * b_x - A_z * B_y * B_z * a_y ) / d
R_23 = ( A_x * A_x * B_z * b_y - A_x * A_z * B_x * b_y
- A_x * B_x * B_z * a_y - A_x * B_y * a_x * b_z
+ A_x * B_y * a_z * b_x + A_y * A_y * B_z * b_y
- A_y * A_z * B_y * b_y + A_y * B_x * a_x * b_z
- A_y * B_x * a_z * b_x - A_y * B_y * B_z * a_y
+ A_z * B_x * B_x * a_y + A_z * B_y * B_y * a_y ) / d
R_31 = ( A_x * B_y * B_y * a_z - A_x * A_y * B_y * b_z
- A_x * A_z * B_z * b_z + A_x * B_z * B_z * a_z
+ A_y * A_y * B_x * b_z - A_y * B_x * B_y * a_z
+ A_y * B_z * a_x * b_y - A_y * B_z * a_y * b_x
+ A_z * A_z * B_x * b_z - A_z * B_x * B_z * a_z
- A_z * B_y * a_x * b_y + A_z * B_y * a_y * b_x ) / d
R_32 = ( A_x * A_x * B_y * b_z - A_x * A_y * B_x * b_z
- A_x * B_x * B_y * a_z - A_x * B_z * a_x * b_y
+ A_x * B_z * a_y * b_x - A_y * A_z * B_z * b_z
+ A_y * B_x * B_x * a_z + A_y * B_z * B_z * a_z
+ A_z * A_z * B_y * b_z + A_z * B_x * a_x * b_y
- A_z * B_x * a_y * b_x - A_z * B_y * B_z * a_z ) / d
R_33 = ( A_x * A_x * B_z * b_z - A_x * A_z * B_x * b_z
- A_x * B_x * B_z * a_z + A_x * B_y * a_x * b_y
- A_x * B_y * a_y * b_x + A_y * A_y * B_z * b_z
- A_y * A_z * B_y * b_z - A_y * B_x * a_x * b_y
+ A_y * B_x * a_y * b_x - A_y * B_y * B_z * a_z
+ A_z * B_x * B_x * a_z + A_z * B_y * B_y * a_z ) / d
In a computer program, the above calculations use floating-point arithmetic, which means rounding errors will creep in. It is therefore often necessary to orthonormalize the matrix (i.e., make sure all row vectors are unit vectors and perpendicular to each other, and all column vectors are unit vectors and perpendicular to each other). This is nontrivial, and is a completely separate issue in general.
One possible approach is to use an unit quaternion (or versor) $Q$, to describe the matrix $\mathbf{R}$. For each pair of vectors $\vec{A}_i, \vec{a}_i$ you calculate the versor $Q_i = (w_i, x_i, y_i, z_i)$: $$\begin{aligned} w &= \frac{\sqrt{1 + R_{11} + R_{22} + R_{33}}}{2} \\ x &= \frac{R_{32} - R_{23}}{4 w} \\ y &= \frac{R_{13} - R_{31}}{4 w} \\ z &= \frac{R_{21} - R_{12}}{4 w} \end{aligned}$$ You can then calculate the sum of the versors, and normalize it to unit length: $$\begin{aligned} d &= \sum_i w_i^2 + x_i^2 + y_i^2 + z_i^2 \\ w &= \frac{1}{d}\sum_i w_i \\ x &= \frac{1}{d}\sum_i x_i \\ y &= \frac{1}{d}\sum_i y_i \\ z &= \frac{1}{d}\sum_i z_i \end{aligned}$$ Note that each versor describes the orientation in axis-angle format. If $(a_x , a_y , a_z)$ is the axis unit vector, and $\theta$ is the rotation angle around that axis, then $$\begin{aligned} \theta &= 2 \operatorname{atan2}(\sqrt{x^2 + y^2 + z^2}, \; w) \\ a_x &= \frac{x}{\sqrt{x^2 + y^2 + z^2}} \\ a_y &= \frac{y}{\sqrt{x^2 + y^2 + z^2}} \\ a_z &= \frac{z}{\sqrt{x^2 + y^2 + z^2}} \\ \end{aligned}$$ where $\operatorname{atan2}(y, x) = \arctan(y/x)$, except that it takes the quadrant into account too (and therefore can report all angles -180° to +180° or 0° to 360°). You can also reconstruct the rotation matrix $\mathbf{R}$ from the normalized sum of unit quaternions, using $$\begin{aligned} R_{11} &= 1 - 2 ( y^2 + z^2) \\ R_{21} &= 2 ( x y + z w ) \\ R_{31} &= 2 ( x z - y w ) \\ \end{aligned} \quad \begin{aligned} R_{12} &= 2 ( x y - z w ) \\ R_{22} &= 1 - 2 ( x^2 + z^2 ) \\ R_{32} &= 2 ( y z + x w ) \\ \end{aligned} \quad \begin{aligned} R_{13} &= 2 ( x z + y w ) \\ R_{23} &= 2 ( y z - x w ) \\ R_{33} &= 1 - 2 ( x^2 + y^2 ) \\ \end{aligned}$$
In some molecular dynamics simulations, temperature controls can cause particle clusters to rotate (making the simulation not very useful). There, the actual problem is that one has a large number of vector pairs (vector triplets!) that follow $$\mathbf{R} (\vec{A}_i - \vec{o}) = \vec{o} + \vec{a}_i + \vec{\epsilon}_i$$ where $\vec{o}$ is a point on the rotation axis (usually the center of mass for the particle cluster), $\vec{A}$ and $\vec{a}$ the locations of the particle at two different time steps, and $\vec{\epsilon}_i$ corresponds to the relative movement of particle $i$ with respect to the cluster or particle clump it is in, between those two time steps.
Although this seems to be the same problem, it is actually quite different, because to cancel the rotation, one wants to find the $\mathbf{R}$ and $\vec{o}$ that minimize $\sum_i^N\lVert\vec{\epsilon}_i\rVert^\gamma$, where $\gamma$ is usually 2 for least squares fitting. Because we don't know what the $\vec{\epsilon}_i$ actally were, even this only yields an approximate, "best guess" rotation matrix. (This is also why the rotation is usually cancelled over a number of time steps, until the average rotation angle per time step is under some specific limit.)
As long as $A$ and $B$ are linearly independent, then $R$ is uniquely determined. For $R(A\times B)=A'\times B'$, and as $A$, $B$, $A\times B$ form a basis for $\Bbb R^n$ then $R$ is determined by its action on these three vectors.