least squares regression in 3space robjohn is giving me a hand with this, but in case anybody else knows...
I need to do a least-squares regression for linearity on a set of coordinates in 3space. If the dataset is linear, I need to see if it is close to vertical or horizontal. How could I do this?
Many thanks in advance
Joe Stavitsky
 A: Typically, vertical would say all the $(x,y)$ coordinates are the same and horizontal would say all the $z$ coordinates are the same.  So you could just look at the standard deviations of the coordinates of all the points to assess vertical or horizontal.  That doesn't check if the points lie on an arbitrary line.
The discussion you cite was indeed to give a plane-that was the hypothesis.  If you find a relation like $ax+by+cz=k$ you get a plane, as one equation reduces the dimension of the space by one.  If you believe your points lie on a line in 3-space, you have two options.  A standard linear regression (where one coordinate is fixed and you minimize the sum squared error in the other direction) may well work for you, and you can just do two regressions, one $x$ vs $y$ and another $x$ vs $z$.  If the equations you get are $y=m_yx+b_y, z=m_zxb_z$, the line is then $(0,b_y,b_z)+k(1,m_y,m_z)$.  If you want orthogonal distance regression, I believe (but didn't follow the derivation to confirm) the link you have is applicable, but you want the line to be in the direction of the maximum singular value, still through the centroid.
A: I posted this this morning in chat to Joe, but since the question is here, too, I will post it here.
Uniquely, up to the sign of $u$, specify a line as $p+ut$ where $|u|=1$, $p\cdot u=0$, and $t\in\mathbb{R}$. The square of the distance of $x$ from this line is
$$
|x-p|^2-((x-p)\cdot u)^2 = |x-p|^2-(x\cdot u)^2\tag{1}
$$
Thus, for a given $\{x_k\}$, we wish to minimize
$$
\Delta=\frac1n\sum_{k=1}^n\left[|x_k-p|^2-(x_k\cdot u)^2\right]\tag{2}
$$
over all $p$ and $u$ so that
$$
|u|=1\text{ and }p\cdot u=0\tag{3}
$$
Thus, we need that for any variation of $u$ and $p$ that maintain $(3)$, the variation of $(2)$ must be $0$.  That is, whenever
$$
u\cdot\delta u=0\quad\text{and}\quad p\cdot\delta u+u\cdot\delta p=0\tag{4}
$$
we want the variation of $\Delta$ to vanish
$$
0=\sum_{k=1}^n\left[(x_k-p)\cdot\delta p+u\cdot x_kx_k\cdot\delta u\right]\tag{5}
$$
Define
$$
\bar{x}=\frac1n\sum_{k=1}^nx_k\quad\text{and}\quad X=\frac1n\sum_{k=1}^nx_kx_k^T\tag{6}
$$
and $(5)$ becomes
$$
0=(\bar{x}-p)\cdot\delta p+(Xu)\cdot\delta u\tag{7}
$$
Considering $(4)$, if we hold $u$ constant and vary $p$, we get that for any $\delta p$ so that $u\cdot\delta p=0$
$$
0=(\bar{x}-p)\cdot\delta p\tag{8}
$$
Because $(8)$ holds for any $\delta p$ which is perpendicular to $u$, we get that $\bar{x}-p\,||\,u$, which leads to
$$
\bar{x}-p=(\bar{x}-p)\cdot uu=\bar{x}\cdot uu\tag{9}
$$
Combining $(4)$ and $(9)$ yields
$$
\bar{x}\cdot\delta u+u\cdot\delta p=0\tag{10}
$$
Plugging $(9)$ into $(7)$ and applying $(10)$ yields
$$
\begin{align}
0
&=\bar{x}\cdot uu\cdot\delta p+(Xu)\cdot\delta u\\
&=-\bar{x}\cdot u\bar{x}\cdot \delta u+(Xu)\cdot\delta u\\
&=\left(\left(X-\bar{x}\bar{x}^T\right)u\right)\cdot\delta u\tag{11}
\end{align}
$$
Since $(11)$ holds for any $\delta u$ perpendicular to $u$, we must have that
$$
\left(X-\bar{x}\bar{x}^T\right)u\,||\,u\tag{12}
$$
That is, $u$ is an eigenvector of $X-\bar{x}\bar{x}^T$.
Plugging $(9)$ into $(2)$ yields $\Delta$ as a function of $u$,
$$
\begin{align}
\Delta
&=\frac1n\sum_{k=1}^n\left[|x_k-p|^2-(x_k\cdot u)^2\right]\\
&=\frac1n\sum_{k=1}^n|(x_k-\bar{x})+\bar{x}\cdot uu|^2-u^TXu\\
&=\frac1n\sum_{k=1}^n\left(|x_k-\bar{x}|^2+2(x_k-\bar{x})\cdot u\bar{x}\cdot u+(\bar{x}\cdot u)^2\right)-u^TXu\\
&=\operatorname{Var}[\{x_k\}]-u^T\left(X-\bar{x}\bar{x}^T\right)u\tag{13}
\end{align}
$$
and since
$$
X-\bar{x}\bar{x}^T=\frac1n\sum_{k=1}^n(x_k-\bar{x})(x_k-\bar{x})^T\tag{14}
$$
we see that $X-\bar{x}\bar{x}^T$ is positive semi-definite. Therefore, the maximum eigenvalue corresponds to the minimum of $\Delta$.
Once we have found $u$, a unit eigenvector of $X-\bar{x}\bar{x}^T$ with the greatest eigenvalue, we can use $(9)$ to get $p=\bar{x}-\bar{x}\cdot uu$.
