First, some background. Let's review the underlying math and notation.
Transformation $\mathbf{T}$ has four free parameters: scale $\lambda$, rotation $\theta$, and translation by $(t_x, t_y)$. Let's define the transformation of point $\vec{p} = (p_x , p_y)$ to point $\vec{q} = (q_x , q_y)$ as
$$\vec{q} = \mathbf{T} \vec{p} \quad \iff \quad
\left[ \begin{matrix} q_x \\ q_y \\ 1 \end{matrix} \right] = \left[ \begin{matrix}
\lambda \cos \theta & -\lambda \sin \theta & t_x \\
\lambda \sin \theta & \lambda \cos \theta & t_y \\
0 & 0 & 1 \end{matrix} \right ] \left [ \begin{matrix} p_x \\ p_y \\ 1 \end{matrix} \right ] \tag{1a}\label{G1a}$$
Defining the transform as a single matrix $\mathbf{T}$ as
$$\mathbf{T} = \left[ \begin{matrix} \lambda \cos \theta & -\lambda \sin \theta & t_x \\
\lambda \sin \theta & \lambda \cos \theta & t_y \\
0 & 0 & 1 \end{matrix} \right ] \tag{1b}\label{G1b}$$
is useful, because we can calculate the inverse ($\vec{p} = \mathbf{T}^{-1} \vec{q}$),
$$\mathbf{T}^{-1} = \left[ \begin{matrix}
\frac{\cos\theta}{\lambda} & \frac{\sin\theta}{\lambda} & - \frac{t_x \cos\theta + t_y \sin\theta}{\lambda} \\
- \frac{\sin\theta}{\lambda} & \frac{\cos\theta}{\lambda} & \frac{t_x \sin\theta - t_y \cos\theta}{\lambda} \\
0 & 0 & 1 \end{matrix} \right] \tag{1c}\label{G1c}$$
and combine transformations via matrix multiplication:
$$\mathbf{T}_\text{combined} = \mathbf{T}_\text{after} \mathbf{T}_\text{before} \tag{1d}\label{G1d}$$
where the rightmost transformation is applied first, and leftmost transformation last. Note that the identity transform, "no transform" or "as-is", is
$$\mathbf{1} = \left[ \begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{matrix} \right] \tag{1e}\label{G1e}$$
This way we can also "ignore" the current transform ($\mathbf{T}_\text{before}$) for the purposes here, and just calculate the transform $\mathbf{T}_\text{after}$ caused by the gesture here; and get the combined transform using $\eqref{G1d}$.
(Note that such $3 \times 3$ matrices (with only $2 \times 3$ elements actually stored) are commonly used in 2D computer graphics, including SVG transform matrix() attribute, with the third component ($1$) of vectors not stored, and just baked in to the functions doing matrix-vector multiplication. Similarly for $4 \times 4$ matrices (with only $3 \times 4$ elements actually stored) for 3D computer graphics, including OpenGL.)
The problem is to find the transformation $\mathbf{T}$ that transforms $\vec{p}_1$ and $\vec{p}_2$ to $\vec{q}_1$ and $\vec{q}_2$,
$$\left\lbrace \, \begin{aligned}
\vec{q}_1 &= \mathbf{T} \vec{p}_1 \\
\vec{q}_2 &= \mathbf{T} \vec{p}_2 \\
\end{aligned} \right. \tag{2}\label{G2}$$
which is actually a set of four equations (two Cartesian components per vector), and four unknowns (that define $\mathbf{T}$), and therefore should have a single solution (except for degenerate cases, where $\vec{p}_1 = \vec{p}_2$ or something similar).
Note that because of finite precision of floating-point numbers, we do not want to apply the transformation incrementally. That is, we should remember the very initial touch points $\vec{p}_1$ and $\vec{p}_2$, keeping them unchanged for the duration of the gesture, also remembering the original transformation $\mathbf{T}_\text{before}$, compute a temporary gesture transformation $\mathbf{T}_\text{gesture}$ based on the initial touch points $\vec{p}_1$ and $\vec{p}_2$ and current touch points $\vec{q}_1$ and $\vec{q}_2$, and a temporary currently applied transformation $\mathbf{T}_\text{temp} = \mathbf{T}_\text{gesture}$. Only when the gesture is finished, will the final $\mathbf{T}_\text{temp}$ become the currently applied and stored transformation.
If you apply the transformation incrementally, i.e. $\vec{p}_1$, $\vec{p}_2$ corresponding to the previous touch locations, and $\vec{q}_1$, $\vec{q}_2$ corresponding to current touch locations, and continously apply $\mathbf{T}_\text{current} = \mathbf{T}_\text{gesture} \mathbf{T}_\text{previous}$, not only do your rounding errors compound (making the gesture control erratic), but it will also be jittery, because small movement of one finger when close to the other will cause large rotations; the faster your touch pad update rate, the more likely such small movements are. The incremental moment-to-moment transformation approach just does not work in practice.
Let $\vec{p}_1$ and $\vec{p}_2$ be the initial touch points, and $\vec{q}_1$ and $\vec{q}_2$ the corresponding final positions:
$$\vec{p}_1 = \left[ \begin{matrix} x_1 \\ y_1 \\ 1 \end{matrix} \right], \;
\vec{p}_2 = \left[ \begin{matrix} x_2 \\ y_2 \\ 1 \end{matrix} \right], \;
\vec{q}_1 = \left[ \begin{matrix} \chi_1 \\ \gamma_1 \\ 1 \end{matrix} \right], \;
\vec{q}_2 = \left[ \begin{matrix} \chi_2 \\ \gamma_2 \\ 1 \end{matrix} \right]$$
Let $C = \lambda \cos \theta$, $S = \lambda \sin \theta$, and translation is by $(X, Y)$, so that the transformation matrix is
$$\mathbf{T} = \left[ \begin{matrix}
C & -S & X \\
S & C & Y \\
0 & 0 & 1 \\
\end{matrix} \right] \tag{3a}\label{G3a}$$
Note that as long $C S \ne 0$, the upper left $2 \times 2$ matrix is a proper rotation and scaling matrix, for all $C \in \mathbb{R}$, $S \in \mathbb{R}$. The vector $(C, S)$ is the new $x$ axis basis vector, and $(-S, C)$ the new $y$ axis basis vector, and the two are always perpendicular (with positive $y$ axis 90° counterclockwise from positive $x$ acis) and have the same length.
Substituting these into $\eqref{G2}$ we have four equations and four unknowns ($C$, $S$, $X$, and $Y$), which has exactly one solution:
$$\left\lbrace \, \begin{aligned}
L^2 &= (x_2 - x_1)^2 + (y_2 - y_1)^2 \\
C &= \frac{ (\chi_2 - \chi_1)(x_2 - x_1) + (\gamma_2 - \gamma_1)(y_2 - y_1) }{L^2} \\
S &= \frac{ (\gamma_2 - \gamma_1)(x_2 - x_1) - (\chi_2 - \chi_1)(y_2 - y_1) }{L^2} \\
X &= \frac{ (x_2 - x_1)(\chi_1 x_2 - \chi_2 x_1) + (y_2 - y_1)(\chi_1 y_2 - \chi_2 y_1) + (\gamma_2 - \gamma_1)(x_2 y_1 - x_1 y_2) }{L^2} \\
Y &= \frac{ (\chi_2 - \chi_1)(x_1 y_2 - x_2 y_1) + (y_2 - y_1)(\gamma_1 y_2 - \gamma_2 y_1) + (x_2 - x_1)(\gamma_1 x_2 - \gamma_2 x_1) }{L^2} \\
\end{aligned} \right. \tag{3b}\label{G3b}$$
If we do need the angle $\theta$ and the scale factor $\lambda$, they are obviously
$$\left\lbrace ~ \begin{aligned}
\lambda &= \sqrt{C^2 + S^2} \\
\theta &= \operatorname{atan2}\left(S, C\right) \\
\cos\theta &= \displaystyle \frac{C}{\sqrt{C^2 + S^2}} \\
\sin\theta &= \displaystyle \frac{S}{\sqrt{C^2 + S^2}} \\
\end{aligned} \right. \tag{3c}\label{G3c}$$
where $\operatorname{atan2}$ is the two-parameter form of arcus tangent, provided by almost all programming languages. ($\operatorname{atan2}(y, x) = \arctan(y/x)$ for $x \gt 0$, but is valid for all $x, y \in \mathbb{R}$ except $x=0, y=0$, as it takes the signs of both $x$ and $y$ into account in other quadrants. Typically, it returns the angle in radians between $-\pi$ (-180°) and $+\pi$ (+180°), but it may vary between programming languages.)
We can verify the above using
$$\vec{o} = \left[\begin{matrix} o_x \\ o_y \\ 1 \end{matrix}\right], ~
\vec{t} = \left[\begin{matrix} t_x \\ t_y \\ 1 \end{matrix}\right], ~
\vec{b} = \left[\begin{matrix} b_x \\ b_y \\ 1 \end{matrix}\right], ~
\vec{a} = \left[\begin{matrix} a_x \\ a_y \\ 1 \end{matrix}\right]$$
and
$$\vec{p}_1 = \vec{o} + \vec{b}, ~ ~
\vec{p}_2 = \vec{o} - \vec{b}, ~ ~
\vec{q}_1 = \vec{o} + \vec{t} + \vec{a}, ~ ~
\vec{q}_2 = \vec{o} + \vec{t} - \vec{a}$$
so that $\vec{o}$ represents the middle point between the two initial touch points, $\vec{b}$ represents the vector from the middle point to the first initial touch point and $-\vec{b}$ the vector from the middle point to the second initial touch point, $\vec{t}$ represents the translation, and $\vec{a}$ represents the vector from the middle point between the final touch points ($\vec{o}+\vec{t}$) to the first final touch point, and $-\vec{a}$ the vector from the middle point between the final touch points to the second final touch point. The transformation matrix is then
$$\mathbf{T} = \left[ \begin{matrix}
\displaystyle \frac{a_x b_x + a_y b_y}{b_x^2 + b_y^2} & \displaystyle \frac{a_x b_y - a_y b_x}{b_x^2 + b_y^2} & \displaystyle t_x + o_x - \frac{o_x (a_x b_x + a_y b_y ) + o_y (a_x b_y - a_y b_x )}{b_x^2 + b_y^2} \\
\displaystyle \frac{a_y b_x - a_x b_y}{b_x^2 + b_y^2} & \displaystyle \frac{a_x b_x + a_y b_y}{b_x^2 + b_y^2} & \displaystyle t_y + o_y - \frac{ o_x ( a_y b_x - a_x b_y ) + o_y (a_x b_x + a_y b_y ) }{b_x^2 + b_y^2} \\
0 & 0 & 1 \end{matrix} \right]$$
and if we do the math, we indeed find that $\mathbf{T} \vec{p}_1 = \vec{q}_1$ and $\mathbf{T} \vec{p}_2 = \vec{q}_2$.
Note that we can extract rotation angle, scale factor, and translation from such a $3 \times 3$ matrix (assuming it is orthogonal, not skewed) with entries
$$\mathbf{M} = \left[ \begin{matrix} m_{11} & m_{12} & m_{13} \\ m_{21} & m_{22} & m_{23} \\ 0 & 0 & 1 \end{matrix} \right] \tag{4}\label{G4a}$$
via
$$\left\lbrace ~ \begin{aligned}
\theta &= \operatorname{atan2}\left( m_{21} - m_{22} ,\, m_{11} + m_{22} \right) \\
\lambda &= \frac{1}{2}\left(\sqrt{m_{11}^2 + m_{21}^2} + \sqrt{m_{12}^2 + m_{22}^2}\right) \\
t_x &= X = m_{13} \\
t_y &= Y = m_{23} \\
\end{aligned} \right. \tag{4b}\label{G4b}$$
In particular, if you have an existing transfrom (by $\theta_1$, $\lambda_1$, $(X_1, Y_1)$), and you wish to further transform it (by $\theta_2$, $\lambda_2$, $(X_2, Y_2)$), first construct the two matrices using
$$\mathbf{M} = \left[ \begin{matrix}
\lambda \cos\theta & -\lambda \sin\theta & X \\
\lambda \sin\theta & \lambda \cos\theta & Y \\
0 & 0 & 1 \end{matrix} \right ] \tag{4c}\label{G4c}$$
then calculate their product via matrix multiplication, original orientation rightmost, latest applied transform leftmost,
$$\mathbf{M}_\text{final} = \mathbf{M}_2 \mathbf{M}_1$$
and then extract the rotation angle $\theta$, scale factor $\lambda$, and translation $(X, Y)$ as specified in $\eqref{G4b}$.
Best Answer
In principle you need two rotations and a matrix-multiplication. Call the recentered matrices Q ad S and P as T.
Then find the rotation s, which rotates S to "triangular" form of its coordinates, so that the coordinates of the first point becomes $[x_{s,1},0,0]$, of the second becomes $[x_{s,2},y_{s,2},0]$. This can be done with your ansatz of using unknowns, if you assume three rotations (for each pair of planes 1 rotation, and gives a 3x3-matrix whose 3'rd column is zero and describes thus a triangle in a plane)
Then do the same with another rotation t which rotes T in the same manner.
Then $T = S \cdot s \cdot t^{-1} = S \cdot r$ where r is the matrix for the complete rotation.
A bit of pseudocode (code in my proprietary MatMate-tool):
Example:
We assume S and T being centered. Let
$ \small \text{ S =} \begin{bmatrix} 14.469944&22.964690&-7.581631\\ -15.275348&5.923432&23.720255\\ 0.805404&-28.888122&-16.138624 \end{bmatrix} $
and
$ \small \text{ T =} \begin{bmatrix} 22.808501&2.515200&16.361035\\ 8.393637&-5.071089&-27.109127\\ -31.202138&2.555889&10.748092 \end{bmatrix} $
Now you can solve for a rotation in y/z-plane, which makes the entry in $S_{1,3}=0$. The rotation-parameters are some cos/sin-values. Apply this and you get
$\small \text{ S}^{(1)} = \begin{bmatrix} 14.469944&24.183840&0.000000\\ -15.275348&-1.811476&24.381471\\ 0.805404&-22.372364&-24.381471 \end{bmatrix}$ Now you can solve for a rotation in x/y-plane, which makes the entry in $S_{1,2}=0$. The rotation-parameters are some other cos/sin-values. Apply this and you get
$\small \text{S}^{(2)} = \begin{bmatrix} 28.182218&0.000000&0.000000\\ -9.397482&12.178056&24.381471\\ -18.784736&-12.178056&-24.381471 \end{bmatrix}$
After that a third set of rotation-parameters cos/sin-values make S triangular. It looks then like this
$\small \text{ S}^{(3)} = \begin{bmatrix} 28.182218&0.000000&0.000000\\ -9.397482&27.253645&0.000000\\ -18.784736&-27.253645&0.000000 \end{bmatrix}$
Because the center of your original triangle was moved to the origin, the last column (the z-coordinates) are zero/not needed, since 3 points can always be placed in a plane.
From the three rotation with their cos/sin-values you can construct a 3x3-rotation-matrix, say s.
The same can be done using the matrix T leading to a rotation-matrix t. If S and T describe in fact the same triangle, only rotated, the results are equal: $T^{(3)}=S^{(3)}$. Then you can use the nice fact, that the inverse of a rotation-matrix is just its transpose, such that with
$\small \text{ tr =} \begin{bmatrix} 0.205215&0.860645&0.466022\\ 0.946329&-0.295966&0.129867\\ 0.249696&0.414359&-0.875190 \end{bmatrix}$ we get $ T = S \cdot tr $
In principle this can also be solved using the concept of pseudoinverses: we demand
$ tr = S^{-1} \cdot T $ . But because S has reduced rank the inverse means to divide by zero. Using SVD-decomposition (see wikipedia) the pseudoinverse can be computed if the inverse of the diagonal matrix of the SVD-factors is used (where zeros are simply left zero). This should lead to the same solution.