Efficiently solving a 2D affine transformation For an affine transformation in two dimensions defined as follows:
$$ 
p_i'=\mathbf{A}p_i \Leftrightarrow \\
\left[
\begin{matrix}
x_i' \\ y_i'
\end{matrix}
\right]
=
\left[
\begin{matrix}
a & b & e \\
c & d & f
\end{matrix}
\right]
\left[
\begin{matrix}
x_i \\ y_i \\ 1
\end{matrix}
\right]
$$
Where $(x_i,y_i), (x_i',y_i')$ are corresponding points, how can I find the parameters $\mathbf A$ efficiently?  
Rewriting this as a system of linear equations, given three points (six knowns, six unknowns):
$$
\textbf{P}\alpha=\textbf{P}' \Leftrightarrow \\
\left[
\begin{matrix}
x_0 & y_0 & 0 & 0 & 1 & 0 \\
0 & 0 & x_0 & y_0 & 0 & 1 \\
x_1 & y_1 & 0 & 0 & 1 & 0 \\
0 & 0 & x_1 & y_1 & 0 & 1 \\
x_2 & y_2 & 0 & 0 & 1 & 0 \\
0 & 0 & x_2 & y_2 & 0 & 1 \\
\end{matrix}
\right]
\left[
\begin{matrix}
a \\ b \\ c \\
d \\ e \\ f
\end{matrix}
\right]
=
\left[
\begin{matrix}
x_0' \\ y_0' \\x_1' \\ y_1' \\x_2' \\ y_2'
\end{matrix}
\right]
$$
Allows the use of an LU decomposition, which can be computed in $O(M(n))$ time, where $M(n)$ is the time to multiply two n×n matrices (according to 1).  
Can the specific structure of the $\mathbf P$ matrix be exploited to utilize the Gaussian elimination to reach the reduced row echelon form (thus solving the system) more efficiently?
Is there a way to symbolically derive the required operations? By hand seems rather cumbersome
Thanks
 A: From three point pairs in general position, we can derive an explicit expression for $A$, namely, $$A = \begin{bmatrix}x_0'&x_1'&x_2'\\y_0'&y_1'&y_2'\\1&1&1\end{bmatrix} \begin{bmatrix}x_0&x_1&x_2\\y_0&y_1&y_2\\1&1&1\end{bmatrix}^{-1}.$$ There might be some efficiencies to be gained by examining ways to compute that.
A: The suggested solutions are certainly much more efficient than the naive matrix inversion, but the result of the Gauss-Jordan elimination seems more efficient than one inversion and one multiplication of a $n=3$ matrix (please correct me if I'm wrong):
$$
\begin{bmatrix}
 1 & 0 & 0 & 0 & 0 & 0 & \frac{x_1' y_0-x_2' y_0 -x_0' y_1 +x_2' y_1 +x_0' y_2 -x_1' y_2}{x_1 y_0-x_2 y_0-x_0 y_1+x_2 y_1+x_0 y_2-x_1 y_2} \\
 0 & 1 & 0 & 0 & 0 & 0 & \frac{x_1' x_0-x_2' x_0 -x_0' x_1 +x_2' x_1 +x_0' x_2 -x_1' x_2}{-x_1 y_0+x_2 y_0+x_0 y_1-x_2 y_1-x_0 y_2+x_1 y_2} \\
 0 & 0 & 1 & 0 & 0 & 0 & \frac{y_1' y_0-y_2' y_0 -y_0' y_1 +y_2' y_1 +y_0' y_2 -y_1' y_2}{x_1 y_0-x_2 y_0-x_0 y_1+x_2 y_1+x_0 y_2-x_1 y_2} \\
 0 & 0 & 0 & 1 & 0 & 0 & \frac{y_1' x_0-y_2' x_0 -y_0' x_1 +y_2' x_1 +y_0' x_2 -y_1' x_2}{-x_1 y_0+x_2 y_0+x_0 y_1-x_2 y_1-x_0 y_2+x_1 y_2} \\
 0 & 0 & 0 & 0 & 1 & 0 & \frac{x_2' x_1 y_0-x_1' x_2 y_0 -x_2' x_0  y_1 +x_0' x_2 y_1 +x_1' x_0 y_2-x_0' x_1 y_2}{x_1 y_0-x_2 y_0-x_0 y_1+x_2 y_1+x_0 y_2-x_1 y_2} \\ 0 & 0 & 0 & 0 & 0 & 1 & \frac{y_2' x_1 y_0-y_1' x_2 y_0 -y_2' x_0  y_1 +y_0' x_2 y_1 +y_1' x_0 y_2-y_0' x_1 y_2}{x_1 y_0-x_2 y_0-x_0 y_1+x_2 y_1+x_0 y_2-x_1 y_2}
\end{bmatrix}
$$
A: One thing that can be done is to block factor by multiplying both sides with permutation matrix so that we get:
$$\begin{bmatrix}x_0&y_0&1&&&\\x_1&y_1&1&&&\\x_2&y_2&1&&&\\&&&x_0&y_0&1\\&&&x_1&y_1&1\\&&&x_2&y_2&1\end{bmatrix}$$
And corresponding right hand side
$$\begin{bmatrix}x_0'\\x_1'\\x_2'\\y_0'\\y_1'\\y_2'\end{bmatrix}$$
And then to utilize $$\begin{bmatrix}M_1&0\\0&M_2\end{bmatrix}^{-1}= \begin{bmatrix}M_1^{-1}&0\\0&M_2^{-1}\end{bmatrix}$$
in some suitable way. Then we have reduced to solve two $3\times 3$ systems which in general is much nicer than one $6 \times 6$.
