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.
[14] set cldezweite=5 // decimal digits for printing [15] n=8 // use 8 data points [16] set randomstart=41 [17] 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 ------------------------------------------ [18] center= meansp(udata) // determine the center as means columnwise ------------------------------------------ center : 53.15625 52.40900 53.94908 ------------------------------------------
[19] 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 ------------------------------------------ [20] 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)
[21] t = gettrans(abwdata,"pca") ; rotate centered data and also the center [22] optabwdata=abwdata*t [23] 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
[24] 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
[25] 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 ------------------------------------------
[added subsequent computations] 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)