total least squares derivation with matrices Taken from a computer vision book: "to minimize the sum of the perpendicular distances between points and lines, we need to minimize $$ \sum_i (ax_i + by_i +c)^2$$ subject to $a^2 +b^2 =1$. Now using a Lagrangian multiplier $\lambda$, we have a solution if $$ \left( \begin{array}{ccc}
\overline{x^2} & \overline{xy} & \overline{x} \\
\overline{xy} & \overline{y} & \overline{y} \\
\overline{x} & \overline{y} & 1 \end{array} \right)\left[ \begin{array}{c}
   a \\
   b \\
   c
  \end{array} \right]  = \lambda \left[ \begin{array}{cc}
   2a\\
   2b \\
   0
  \end{array} \right]$$
How is the book getting these matrices? 
Also, the notion is that $\overline{u} = \frac{\sum u_i}{k}$. (Yeah, I don't know what $k$ stands for. I can only assume this is an average.)
It goes onto say that $c = -a\overline{x} - b\overline{y}$, and that we can substitute this back to get the eigenvalue problem $$\left[ \begin{array}{cc}
   \overline{x^2} -\overline{x}~\overline{x} &  \overline{xy} -  \overline{x}\overline{y}\\
    \overline{xy} -  \overline{x}\overline{y} &  \overline{y^2} -  \overline{y} ~\overline{y} \\
  \end{array} \right] \left[\begin{array}{cc}
   a\\
   b
  \end{array} \right] = \mu \left[ \begin{array}{cc}
   a\\
   b
  \end{array} \right].$$
I don't see what they substituted into, and how the answer is derived.

We want to minimize
$$
S=\sum_i(ax_i+by_i+c)^2\tag{1}
$$
Varying $(1)$ yields
$$
\frac12\delta S=
\begin{bmatrix}\delta a&\delta b&\delta c\end{bmatrix}
\left(\sum_i\begin{bmatrix}x_i\\y_i\\1\end{bmatrix}\begin{bmatrix}x_i&y_i&1\end{bmatrix}\right)
\begin{bmatrix}a\\b\\c\end{bmatrix}\tag{2}
$$
Since, $a^2+b^2=1$, we restrict the variations to those so that
$$
0=\begin{bmatrix}\delta a&\delta b&\delta c\end{bmatrix}\begin{bmatrix}a\\b\\0\end{bmatrix}\tag{3}
$$
We want to find $a,b,c$ that cancel $(2)$ under the constraints in $(3)$. Orthogonality implies that we need
$$
\left(\sum_i\begin{bmatrix}x_i\\y_i\\1\end{bmatrix}\begin{bmatrix}x_i&y_i&1\end{bmatrix}\right)
\begin{bmatrix}a\\b\\c\end{bmatrix}
=\lambda\begin{bmatrix}a\\b\\0\end{bmatrix}\tag{4}
$$
Equation $(4)$ is the first equation in the question.
The bottom row of $(4)$ says
$$
\left(\sum_i\begin{bmatrix}x_i&y_i&1\end{bmatrix}\right)
\begin{bmatrix}a\\b\\c\end{bmatrix}=0\tag{5}
$$
The top two rows of $(4)$ say
$$
\left(\sum_i\begin{bmatrix}x_i\\y_i\end{bmatrix}\begin{bmatrix}x_i&y_i&1\end{bmatrix}\right)
\begin{bmatrix}a\\b\\c\end{bmatrix}
=\lambda\begin{bmatrix}a\\b\end{bmatrix}\tag{6}
$$
we can rewrite $(6)$ as
$$
\left(\sum_i\begin{bmatrix}x_i\\y_i\end{bmatrix}\begin{bmatrix}x_i&y_i\end{bmatrix}\right)
\begin{bmatrix}a\\b\end{bmatrix}
+\left(\sum_i\begin{bmatrix}x_i\\y_i\end{bmatrix}\right)c
=\lambda\begin{bmatrix}a\\b\end{bmatrix}\tag{7}
$$
and $(5)$ as
$$
\left(\sum_i\begin{bmatrix}x_i&y_i\end{bmatrix}\right)
\begin{bmatrix}a\\b\end{bmatrix}+\left(\sum_i1\right)c=0\tag{8}
$$
If we set $n=\sum\limits_i1$ and substitute $c$ computed from $(8)$ into $(7)$, we get
$$
\left(\sum_i\begin{bmatrix}x_i\\y_i\end{bmatrix}\begin{bmatrix}x_i&y_i\end{bmatrix}-\frac1n\sum_i\begin{bmatrix}x_i\\y_i\end{bmatrix}\sum_i\begin{bmatrix}x_i&y_i\end{bmatrix}\right)
\begin{bmatrix}a\\b\end{bmatrix}
=\lambda\begin{bmatrix}a\\b\end{bmatrix}\tag{9}
$$
Equation $(9)$ is the last equation in the question.
