Well, symbolically I could rewrite this as $$\left(\begin{array}{ccc} a & b & (\mathbf{Q_1} - \mathbf{Q_2})\\ c & d & (\mathbf{Q_2} - \mathbf{Q_3}) \end{array}\right) \left(\begin{array}{c} \mathbf{P_1}\\ \mathbf{P_2}\\ \lambda \end{array}\right) = \left(\begin{array}{c} e\mathbf{Q_1}\\ f\mathbf{Q_2} \end{array}\right),$$ though I'm not sure how useful this is as the matrix now contains both scalars and vectors (and re-writing the vectors as $3$ scalars each does not seem to make much sense, unless I'm missing something).

Comment was too long. Wanted to add the following — I'm curious about your approach (though additional hints are also welcome).

@Ailurus scalar multiplication on matrices/vectors involves a "hidden"/implicit identity matrix. That is often how multiplication by scalar is defined.

Ah, so I replace $a$ and other scalars by $aI$ in the matrix, turning it into a $6 \times 7$ matrix? That should work :)

However, using the approach based on normal equations now also requires an initial `guess' for $\lambda$, which influences the eventual solution (it finds the solution closest — in the sense of the 2-norm — to the initial guess). Any ideas on how to avoid this?

It seems linear in all aspects. I don't see why an initial guess would be needed. It is not like the right hand side is $0$ or anything.

Right, but the question is how to get $\mathbf{P_i}$ closest to the reference points $\mathbf{R_i}$ by optimizing $\lambda$ (or in some cases, multiple scalar variables). Which approach would you then recommend for this?

Ah wait I see now. Then you can add a term $\|\sum_j s_{ji} R_j-P_i\|_2, \|\sum_j s_{ji} -1\|_2$ to the minimization and the $s_{ji}$ to the vectorization. Probably you will need to reweight to force $s_{ji}$ to be 0 or 1

Hi, MathSE proposed to move to a chat instead. Although I (more or less) understand the building blocks you're suggesting to use, I don't have an overview yet. Are the $s_{ji}$ comparable to Lagrange multipliers, somehow resulting in a penalty term? If you could add a short example, or some names of the (numerical) optimization method(s) you have in mind, that'd be great :)

maybe i can spend some time on it later after work

Ok, cheers!