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

• Look up orthogonal distance regression. – J. M. is a poor mathematician May 11 '12 at 0:58
• @J.M. this discussion mathforum.org/library/drmath/view/63765.html seems to say final output is a plane, not a line. – Joe Stavitsky May 11 '12 at 1:46
• I haven't read through the link (sorry), but if you have a plane in ${\mathbb R}^3,$ then you can extract the unit vector in the direction of the normal to the plane: $\vec n.$ Indeed, if the plane is vertical or horizontal, then you can see that immediately from $\vec n.$ – user2468 May 11 '12 at 3:21
• I am afraid the discussion I posted is over my head mathwise; could sombody post some pseudocode maybe? Thanks again – Joe Stavitsky May 11 '12 at 3:40
• I think Principle component analysis is what you want. Essentially what you want is that the line passes through the centroid in the direction of maximal variance; the line you want corresponds to the first PCA vector, shifted to pass through the centroid. – guy May 11 '12 at 3:50

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.

• This might work. A bit later I will test this and mark accept if I have no additional problems. Thanks so much for the help. – Joe Stavitsky May 11 '12 at 3:34
• I am having trouble understanding the last ewuation - how is that a line? Thanks again – Joe Stavitsky May 11 '12 at 3:38
• @JoeStavitsky: $k$ is the parameter that expresses distance along the line. The points on the line are then $(k,km_y+b_y,km_z+b_z)$ where $k$ is any real number. One way to define a line in 3-space is to give a point and a vector along the line. The whole line is then that point plus any constant times the direction. – Ross Millikan May 11 '12 at 3:40

The core of the orthogonal regression is expressible by a "principal component" analysis as implemented in many statistical programs. You do a "PCA"-rotation on your centered dataset and use the first two "factors" (=first two principal components) as x and y-axes and the z-axis is then the measure of the vertical distance of the datapoints to that plane. The eigenvalue associated to that third component is then equivalent to the sum-of-squares of the vertical deviances and is minimized by the PC-rotation.

Here is a listing of a calculation (using my matrix-calculator "MatMate" and its own language, sorry, don't have Math'ica or Maple...). Also the occuring matrices are documented. The bracketed numbers are only statement numbers in the listing.

 set cldezweite=5  // decimal digits for printing
 n=8  // use 8 data points
 set randomstart=41
 udata  = randomu(3,n) '   // get 3 columns of n rows random coordinates
------------------------------------------
udata :
x-axis      y-axis          z-axis
------------------------------------------
28.65775     4.77896    86.38809
8.65334    62.09805    20.01654
88.59632    78.99541    56.30901
33.55478    18.60703    90.15799
65.22586    81.45162    63.61901
87.68102    39.27773     4.50016
45.83511    48.76525    70.40108
67.04584    85.29792    40.20077
------------------------------------------

 center= meansp(udata)  // determine the center as means columnwise
------------------------------------------
center :
53.15625    52.40900    53.94908
------------------------------------------

 abwdata= udata - center  // generate centered data
------------------------------------------
abwdata :
-24.49850   -47.63004   32.43901
-44.50292    9.68905    -33.93254
35.44007    26.58641     2.35993
-19.60147   -33.80197   36.20891
12.06961    29.04263     9.66992
34.52477    -13.13127   -49.44892
-7.32114    -3.64375    16.45200
13.88958    32.88893    -13.74831
------------------------------------------

 chk = sqsumsp(abwdata)/n   // check the square-sums of each column
------------------------------------------
chk :
725.63213   790.34643   814.84356    // we have s² = 814.84 in the z-axis
------------------------------------------

Get rotationmatrix for principal components rotation (columnwise)

 t = gettrans(abwdata,"pca")
; rotate centered data and also the center
 optabwdata=abwdata*t
 optcenter = center*t
------------------------------------------
optabwdata :
x-axis      y-axis           z-axis
------------------------------------------
-61.37125   -0.48959    -12.42644
2.16570    50.81196    25.28358
33.93348    -28.56924   -0.85030
-52.31193   -9.23432    -4.00819
18.95788    -20.78838   17.06163
38.42755    20.42223    -43.77032
-15.60503   -7.60057     6.02151
35.80360    -4.55208    12.68852

optcenter :
29.71799    -84.50733   21.40434
------------------------------------------

Translate the rotated data to the (rotated) center

 optdata = optabwdata + optcenter
------------------------------------------
optdata :
-31.65326   -84.99691    8.97790
31.88369    -33.69537   46.68792
63.65147    -113.07656  20.55405
-22.59394   -93.74165   17.39615
48.67587    -105.29571  38.46597
68.14554    -64.08510   -22.36598
14.11296    -92.10790   27.42586
65.52159    -89.05941   34.09286
------------------------------------------

Check squaresums in the three new axes. The z-axis(3'rd column) "is minimal" and contains a version of the least-squares-error of the related orthogonal regression

 chk1=sqsumsp(abwsp(optdata))/n
------------------------------------------
chk1 :
1377.57898  551.41041   401.83274  // we have s² = 401.83 in the z-axis
// this should be the least possible
// using rotations and translations only
------------------------------------------

; list the rotation-matrix for the columnwise rotation
------------------------------------------
t :
0.52648    -0.59715    -0.60517
0.62564    -0.20985     0.75136
-0.57567    -0.77419     0.26311
------------------------------------------

It is then possible to use the inverse(which is simply the transpose) of t for a function for the plane (or more correctly for the line as requested in your original question) .

For the line we had then for coordinates (x,y,z)' in the framework of the given points by a free real parameter x

(x,y,z)' = (x,0,0)* t' + center
= (x,0,0)* ( 0.52648   0.62564    -0.57567) + center
-0.59715  -0.20985    -0.77419)
-0.60517   0.75136     0.26311)
= (x*0.52648 ,x*0.62564, x* -0.57567) + center

For the plane we had then for coordinates (x,y,z)' in the framework of the given points by two free real parameters (x,y)

(x,y,z)' = (x,y,0)* t' + center
= (x,y,0)* ( 0.52648   0.62564    -0.57567) + center
-0.59715  -0.20985    -0.77419)
-0.60517   0.75136     0.26311)
= (x*0.52648 -y*0.59715,x*0.62564-y*0.20985, -x*0.57567-y*0.77419)
+ center

Here the formulae might be scaled by an arbitrary factor<>0 because (x,y) are arbitrary reals and we might factor out any number to make the decimal expansion of the coefficients at x any y in the equation more smooth.

(you can download MatMate and reproduce the computation)

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$.