Doing a "linear regression" with 2d points under a linear transformation Let's say I have a list of 2D points $\vec v_i = (x_i,y_i)$ and the list of points after some linear transformation $\vec v_i'= A\vec v + b + \epsilon_i$ where $\epsilon_i$ is some deviation assumed to be small compared to the position. Is there a way to find the "best fit" matrix $A$ and displacement $b$, like there is in normal linear regression?
Example in R:
points <- matrix(runif(2000), nrow = 2)
plot.points <- function(p) {
  transpose <- t(p)
  qplot(x = transpose[,1], y = transpose[,2]) + xlim(c(-5,5)) + ylim(c(-5,5))
}

set.seed(34)
A = matrix(runif(4, -1, 1), ncol = 2)

           [,1]       [,2]
[1,] -0.1104629  0.7697880
[2,]  0.9970808 -0.5231481

## b is just (2,2)
b = matrix(rep(2, 2000), nrow =2)

library(ggplot2)
plot.points(points)
plot.points(A%*%points + b)

## is there a way to reverse engineer A and b 
## assuming we have points and A%*%points + b?

 A: A simple way that avoids having to solve equations or invert matrices is by defining the linearly independent functions $f_1(x,y) = 1$, $f_2(x,y) = x$, and $f_3(x,y) = y$ and then using Gram-Schmidt to create an orthogonal basis for the 3-dimensional vector space spanned by the 3 functions. If we take the first basis vector $e_1(x,y)$ to be $f_1(x,y) = 1$, then the next one $e_2(x,y)$ becomes:
$$e_2(x,y) = f_2(x,y) - \frac{\left\langle f_2,e_1\right\rangle}{\left\langle e_1,e_1\right\rangle}e_1(x,y) = x - \frac{1}{N}\sum_{i=1}^{N}x_i = x - \left\langle x\right\rangle$$
The third basis vector is:
$$\begin{split}
e_3(x,y) &= f_3(x,y) - \frac{\left\langle f_3,e_1\right\rangle}{\left\langle e_1,e_1\right\rangle}e_1(x,y) - \frac{\left\langle f_3,e_2\right\rangle}{\left\langle e_2,e_2\right\rangle}e_2(x,y)\\ &= y - \left\langle y\right\rangle - \frac{\left\langle xy\right\rangle - \left\langle x\right\rangle \left\langle y\right\rangle }{\left\langle x^2\right\rangle - \left\langle x\right\rangle^2}\left(x - \left\langle x\right\rangle\right)
\end{split}
$$
The best fit to the data is given by the projection on the 3-dimensional linear space of each of the two components of $\vec{v}' = (q_1,q_2)$. If the errors were exactly zero, we would have:
$$q_1(x,y) = \frac{\left\langle q_1,e_1\right\rangle}{\left\langle e_1,e_1\right\rangle}e_1(x,y) + \frac{\left\langle q_1,e_2\right\rangle}{\left\langle e_2,e_2\right\rangle}e_2(x,y) +\frac{\left\langle q_1,e_3\right\rangle}{\left\langle e_3,e_3\right\rangle}e_3(x,y)  $$
and
$$q_2(x,y) = \frac{\left\langle q_2,e_1\right\rangle}{\left\langle e_1,e_1\right\rangle}e_1(x,y) + \frac{\left\langle q_2,e_2\right\rangle}{\left\langle e_2,e_2\right\rangle}e_2(x,y) +\frac{\left\langle q_2,e_3\right\rangle}{\left\langle e_3,e_3\right\rangle}e_3(x,y)$$
in case of random errors this yields the least squarest fit. The inner products are readily expressed in terms of averages of the measured quantities, so, all the coefficients are easily obtained.
