inv(A'*A)*A'*b looks like a least squares approximation. You can approximate
$$ A\cdot x \approx b $$
by solving
$$ \left(A^T\cdot A\right)\cdot x = A^T\cdot b $$
which can be solved for $x$ as
$$ x = \left(A^T\cdot A\right)^{-1}\cdot A^T\cdot b $$
You can also fiond this forumla in the section Model verification by linear least squares in the Wikipedia article on SIFT.
Let's also translate these matrices into math notation:
$$
M_1= \begin{bmatrix}
1&0&0&0&0&0\\
0&1&0&0&0&0\\
0&0&0&0&0&0\\
0&0&0&0&0&0\\
0&0&0&0&1&0\\
0&0&0&0&0&0
\end{bmatrix}
\qquad
M_2= \begin{bmatrix}
0&0&0&0&0&0\\
0&0&0&0&0&0\\
0&0&1&0&0&0\\
0&0&0&1&0&0\\
0&0&0&0&0&0\\
0&0&0&0&0&1
\end{bmatrix}
$$
The matrices pt1 and pt2 seem to look something like this:
$$
\textit{pt}_1 = \begin{bmatrix}
x_1 & y_1 & 1 & 1 & 1 & 1 \\
x_2 & y_2 & 1 & 1 & 1 & 1 \\
\vdots & & & & & \vdots \\
x_n & y_n & 1 & 1 & 1 & 1
\end{bmatrix}
\qquad
\textit{pt}_2 = \begin{bmatrix}
1 & 1 & x_1 & y_1 & 1 & 1 \\
1 & 1 & x_2 & y_2 & 1 & 1 \\
\vdots & & & & & \vdots \\
1 & 1 & x_n & y_n & 1 & 1
\end{bmatrix}
$$
The above formulation disregards the difference between pointand pts. Multiplying those $M$ matrices from the right sets certain columns to zero, so you end up with
$$
t_1 = \begin{bmatrix}
x_1 & y_1 & 0 & 0 & 1 & 0 \\
x_2 & y_2 & 0 & 0 & 1 & 0 \\
\vdots & & & & & \vdots \\
x_n & y_n & 0 & 0 & 1 & 0
\end{bmatrix}
\qquad
t_2 = \begin{bmatrix}
0 & 0 & x_1 & y_1 & 0 & 1 \\
0 & 0 & x_2 & y_2 & 0 & 1 \\
\vdots & & & & & \vdots \\
0 & 0 & x_n & y_n & 0 & 1
\end{bmatrix}
$$
The reshape then stacks these matrices.
$$
A = \begin{bmatrix}
x_1 & y_1 & 0 & 0 & 1 & 0 \\
x_2 & y_2 & 0 & 0 & 1 & 0 \\
\vdots & & & & & \vdots \\
x_n & y_n & 0 & 0 & 1 & 0 \\
0 & 0 & x_1 & y_1 & 0 & 1 \\
0 & 0 & x_2 & y_2 & 0 & 1 \\
\vdots & & & & & \vdots \\
0 & 0 & x_n & y_n & 0 & 1
\end{bmatrix}
$$
These are the same lines you also find in the Wikipedia article, although Wikipedia has them interleaved, so they are in a different order.
Now the right side is built from npoint. Despite the comment, I assume that this is not point indices, but instead the new positions of the points, after the transformation was applied. Wikipedia calls this $u$ and $v$. So you have
$$
b = \begin{bmatrix}
u_1 \\
u_2 \\
\vdots \\
u_n \\
v_1 \\
v_2 \\
\vdots \\
v_n
\end{bmatrix}
$$
So what you are doing is you are solving the following least squares approximation problem:
$$
\begin{bmatrix}
x_1 & y_1 & 0 & 0 & 1 & 0 \\
x_2 & y_2 & 0 & 0 & 1 & 0 \\
\vdots & & & & & \vdots \\
x_n & y_n & 0 & 0 & 1 & 0 \\
0 & 0 & x_1 & y_1 & 0 & 1 \\
0 & 0 & x_2 & y_2 & 0 & 1 \\
\vdots & & & & & \vdots \\
0 & 0 & x_n & y_n & 0 & 1
\end{bmatrix}
\cdot
\begin{bmatrix}
m_1 \\ m_2 \\ m_3 \\ m_4 \\ t_x \\ t_y
\end{bmatrix}
\approx
\begin{bmatrix}
u_1 \\
u_2 \\
\vdots \\
u_n \\
v_1 \\
v_2 \\
\vdots \\
v_n
\end{bmatrix}
$$
which is just another way of writing
$$
\begin{bmatrix}m_1 & m_2 \\ m_3 & m_4\end{bmatrix}
\cdot\begin{bmatrix}x_i \\ y_i\end{bmatrix}
+\begin{bmatrix}t_x \\ t_y\end{bmatrix}
\approx\begin{bmatrix}u_i \\ v_i\end{bmatrix}
$$
So in the end, your transformation will be a vector of six elements, namely $(m_1, m_2, m_3, m_4, t_x, t_y)$, which define the affine transformation
$$
\begin{bmatrix}x \\ y\end{bmatrix} \mapsto
\begin{bmatrix}m_1 & m_2 \\ m_3 & m_4\end{bmatrix}
\cdot\begin{bmatrix}x \\ y\end{bmatrix}
+\begin{bmatrix}t_x \\ t_y\end{bmatrix}
$$
Note that I wouldn't program things that way. I'd rather treat this as two systems of equations, since rows related to $u$ don't influence those related to $v$. So I'd write something like
$$
\begin{bmatrix}
x_1 & y_1 & 1 \\
x_2 & y_2 & 1 \\
\vdots & & \vdots \\
x_n & y_n & 1
\end{bmatrix}
\cdot
\begin{bmatrix}
m_1 & m_3 \\
m_2 & m_4 \\
t_x & t_y
\end{bmatrix}
\approx
\begin{bmatrix}
u_1 & v_1 \\
u_2 & v_2 \\
\vdots & \vdots \\
u_n & v_n
\end{bmatrix}
$$
This should translate into matlab like this:
pts = [x;y]';
A = [pts, ones(size(pts,1),1)];
transformation = (inv(A'*A)*A'*npoint)'
