Your understanding is correct. Given two unit vectors $a$ and $b$, their geometric product $R = a b$ is a rotor. The rotor is:
$R = a b = a \cdot b + a \wedge b$
$R = \cos \theta + I \sin \theta$
where $I = \frac{a \wedge b}{\|a \wedge b\|}$. When $R$ is applied to a blade $X$ using the versor product $R X \tilde R$, the total angle of rotation is $2 \theta$. For convenience the rotor is defined to rotate $\theta/2$ so that when the versor product is applied the total rotation angle amounts to $\theta$.
It can be seen algebraically developing the versor product:
$R X \tilde R = (\cos \theta + I \sin \theta) X (\cos \theta - I \sin \theta)$
$R X \tilde R = \cos^2\theta \ X + (\cos \theta \sin \theta)( −X I + I X ) −
\sin^2\theta \ I X I $
Using the fact that $- X I + I X = -2 X \cdot I$
$R X \tilde R = \cos^2\theta \ X - 2(\cos\theta \sin\theta) X \cdot I − \sin^2\theta \ I X I$
It is convinient for our purpose to express the reflection $I X I$ of the blade $X$ with respect to the dual of the bivector $I^* = I e_{321}$ instead of $I$, so we get $I X I = (-1)^k I^* X I^* = -2 I^* (I^* \cdot X) + X$. where $k$ is the grade of the blade $X$. So we get:
$R X \tilde R = (\cos^2\theta - \sin^2\theta) \ X - 2(\cos\theta \sin\theta) \ X \cdot I + 2 \sin^2\theta \ I^* (I^* \cdot X)$
Applying the following trigonometric identities:
$\cos 2\theta = \cos^2\theta − \sin^2\theta$
$\sin 2\theta = 2 \cos\theta \sin\theta$
$1−\cos 2\theta = 2 \sin^2\theta$
we finally get:
$R X \tilde R = \cos 2\theta \ X - \sin 2\theta \ X \cdot I + (1 − \cos 2\theta) \ I^* (I^* \cdot X)$
This is the Geometric Algebra version of Rodrigues' formula. As can be seen, the versor product produces a rotation with total angle of rotation of $2 \theta$. That is why the convention is to take angle $\theta/2$.
I think of quaternions as actions, not states. If I want to represent the state of an object's rotation by a quaternion, this must always be in reference to some unrotated state.
So when I see that the rotation of some object is represented by the quaternion $q_{current}$, I interpret it as "the object has been rotated away from its initial state by $q_{current}$."
If I want to reconfigure the object so it has rotation $q_{target}$, what I mean is, "I want to apply some rotation $q_{delta}$ to the object so that the result has been rotated from its initial state by $q_{new}$." When interpreted this way, we see that first we need to rotate back to the initial state, then to its new state:
$$ q_{delta} = q_{target}q_{current}^{-1} $$
Remember that functions act on the left! The quaternions are acting on the object, so the actions should be read right-to-left: first we do $q_{current}^{-1}$, then we do $q_{target}$. As you say in your comment above, this why matrix libraries often define combination order from right to left --- they're also often interpreted as acting on vectors.
Does this make sense? Let's check:
$$ q_{new} = q_{delta}q_{current} = q_{target}q_{current}^{-1}q_{current} = q_{target} $$
just as desired.
Note: Your (3) is not the same as the previous equation! In fact, quaternions do not commute, so the equation $q_{new} = q_{current}^{-1}q_{target}q_{current}^{-1}$ all but guarantees that $q_{new}$ will not be equal to $q_{target}$.
Best Answer
First, you must understand basic linear algebra and matrix-vector multiplication. For example, $$\left [ \begin{matrix} m_{11} & m_{12} & m_{13} \\ m_{21} & m_{22} & m_{23} \\ m_{31} & m_{32} & m_{33} \end{matrix} \right ] \left [ \begin{matrix} x_1 \\ x_2 \\ x_3 \end{matrix} \right ] = \left [ \begin{matrix} y_1 \\ y_2 \\ y_3 \end{matrix} \right ]$$ is the same as the system of linear equations $$\left\lbrace ~ \begin{aligned} m_{11} x_1 + m_{12} x_2 + m_{13} x_3 &= y_1 \\ m_{21} x_1 + m_{22} x_2 + m_{23} x_3 &= y_2 \\ m_{31} x_1 + m_{32} x_2 + m_{33} x_3 &= y_3 \\ \end{aligned} \right .$$ except that "Linear Algebra" not only makes it easier to write, but also contains a lot of tools how to efficiently manipulate such systems, and to solve them. (Among other things, of course.)
In general, when we have some matrix $\mathbf{M}$, it has some eigenvectors $\vec{v}_k$ and corresponding eigenvalues $\lambda_k$, such that $$\mathbf{M} \vec{v}_k = \lambda_k \vec{v}_k$$ In other words, when multiplied by the matrix, eigenvectors are those vectors that get only scaled by the corresponding eigenvalue without any change in "direction".
So, if you happened to have a pure rotation/reflection matrix $\mathbf{R}$, such that $$\vec{p}_i = \mathbf{R} \vec{q}_i$$ the axis of rotation of $\mathbf{R}$ is the eigenvector corresponding to the eigenvalue closest to $1$. (The other eigenvalues are often complex; your library may or may not provide those too.)
Singular values are the absolute values of the eigenvalues, $\lvert\lambda_k\rvert$.
Singular value decomposition "decomposes" the matrix $\mathbf{M}$ with all real components (so this will not apply if you use complex numbers!) into three parts, two unitary matrices $\mathbf{U}$ and $\mathbf{V}^T$ (where ${}^T$ denotes transpose, i.e. replacing columns with rows and vice versa, or rotating around the descending diagonal), and a diagonal matrix $\mathbf{\Sigma}$, $$\mathbf{M} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^{T}$$ This does not immediately look useful, but it turns out that the three parts have very useful properties.
Of particular use is the pseudoinverse: $$\mathbf{M}^{+} = \mathbf{V} \mathbf{\Sigma}^{+} \mathbf{U}^T$$ where $\mathbf{\Sigma}^{+}$ is $\mathbf{\Sigma}$ but all nonzero entries replaced with their reciprocals, i.e. $$\Sigma_{k k}^{+} = \begin{cases} \frac{1}{\Sigma_{k k}}, & \Sigma_{k k} \ne 0 \\ 0, & \Sigma_{k k} = 0 \\ \end{cases}$$
The pseudoinverse is extremely useful, because if you know some matrix $\mathbf{M}$ and vector $\vec{y}$, and want to find vector $\vec{x}$, such that $$\mathbf{M} \vec{x} = \vec{y}$$ and you have the pseudoinverse $\mathbf{M}^{+}$, then $$\vec{x} = \mathbf{M}^{+} \vec{y}$$
In other words, when the elements of $\mathbf{M}$ are real numbers, you can solve the $\mathbf{M} \vec{x} = \vec{y}$ problem for $\vec{x}$ by first doing singular value decomposition of matrix $\mathbf{M}$, $$\mathbf{M} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^T$$ so that the solution is $$\vec{x} = \mathbf{M}^{+} \vec{y} = \mathbf{V} \mathbf{\Sigma}^{+} \mathbf{U}^T$$ where $\mathbf{\Sigma}^{+}$ is calculated as previously mentioned.
As a practical problem, there are some crucial key points that should affect ones approach:
Are the point sets ordered or not? That is, does point $\vec{p}_i$ correspond to $\vec{q}_i$, or to some $\vec{q}_k$, with the mapping between $i$ and $k$ indexes unknown? If the points are extracted from e.g. images, the method of extraction will dictate whether the index is the same in both original and rotated point clouds.
Is there translation in addition to rotation? This complicates the picture somewhat, adding three new variables (in addition to the three/four describing the pure rotation matrix) to be solved.
Is there additional per-point movement? That is, if the point cloud is not rigid, but may deform or change between the two states, the problem becomes much harder to solve. Iterative methods recommeded.
The most generic form of the point cloud problem is $$\vec{a}_i + \vec{\epsilon}_i = \vec{t} + \mathbf{R} \vec{b}_i$$ where $\vec{b}_i$ and $\vec{a}_i$ are the two known locations of point $i$, with $\vec{t}$ a translation, and $\vec{\epsilon}_i$ some point-specific error or movement between the two orientations, is hard. There are solutions (in e.g. molecular dynamics simulations, see "rotation detection" and "rotation elimination"), but they are approximate, and will benefit from iterative refinement. That is, do not try to get a perfect result from the get go, but rather ensure you always refine the result to slightly better.
If we consider only pure rotation, i.e. $$\vec{a}_i = \mathbf{R} \vec{b}_i$$ where $\vec{b}_i$ is the position before rotation, and $\vec{a}_i$ position after rotation, of point $i$, then we can look at Wahba's problem, minimizing $J(\mathbf{R})$, $$J(\mathbf{R}) = \frac{1}{2} \sum_{i=1}^{N} w_i \left\lVert \vec{a}_i - \mathbf{R} \vec{b}_i \right\rVert^2$$ except that $w_i = 2 ~ \forall ~ i$.
In other words, we're trying to find the solution $\mathbf{R}$ where the sum of squared errors in distances to $\vec{a}_i$ after rotation, is minimized. Or, writing the problem as $$\vec{a}_i + \vec{\epsilon}_i = \mathbf{R} \vec{b}_i, \quad i = 1 .. N$$ we are minimizing $\sum \lVert\vec{\epsilon}_i\rVert^2$. In a perfect world, we'd minimize it to zero.
The Wikipedia page explains that the solution is to first construct construct matrix $\mathbf{B}$, $$\mathbf{B} = \sum_{i=1}^{N} \vec{a}_i \vec{b}_i^T$$ where $\vec{b}_i^T$ means the position of point $i$ before rotation as a horizontal vector, $\vec{a}_i$ being a vertical vector, using matrix-matrix multiplication. In other words, $$\mathbf{B} = \sum_{i=1}^{N} \left [ \begin{matrix} X_i x_i & X_i y_i & X_i z_i \\ Y_i x_i & Y_i y_i & Y_i z_i \\ Z_i x_i & Z_i y_i & Z_i z_i \\ \end{matrix} \right ], \vec{a}_i = \left [ \begin{matrix} X_i \\ Y_i \\ Z_i \end{matrix} \right ], \vec{b}_i = \left [ \begin{matrix} x_i \\ y_i \\ z_i \end{matrix} \right ]$$ Then, obtain the singular value decomposition of $\mathbf{B}$: $$\mathbf{B} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^T$$ Calculate the (product of the) determinants of $\mathbf{U}$ and $\mathbf{V}$, and form a new $3 \times 3$ matrix $\mathbf{T}$: $$\mathbf{T} = \left [ \begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & \det(\mathbf{U})\det(\mathbf{V}) \end{matrix} \right ]$$ Then, the sought after rotation matrix $\mathbf{R}$ is $$\mathbf{R} = \mathbf{U} \mathbf{T} \mathbf{V}^T$$