3D Linear Regression

Question

I want to write a program that, given a list of points in 3D-space, represented as an array of x,y,z coordinates in floating point, outputs a best-fit line in this space. The line can/should be in the form of a unit vector and a point on the line.

The problem is that I don't know how this is to be done. The closest thing I found was this link, though quite honestly I did not understand how he went from equation to equation and by the time we got to matrices I was pretty lost.

Is there a generalization of simple 2D linear regression that I can use/can someone explain (mathematically) if/how the above linked-to method works (and what one would have to do to compute the best-fit line using it)?

Answer 1

Great answer by @WaTeim

here is my contribution in python for those who need it. works with the numerical example provided

def regr(X):
    y= np.average(X, axis=0)
    Xm = X-y
    u, s, v = np.linalg.svd((1./X.shape[0])*np.matmul(Xm.T,Xm))

    # Extra Credit: Going back
    z= np.matmul(u[:,0].T, Xm.T)
    c = np.array([z*n for n in u[:,0]])
    d = np.array(y.tolist()*c.shape[1]).reshape(c.shape[1],-1).T
    e = (c+d).T
    return u,s,v

regr(np.array([[6, 4, 11],[8,5,15],[12,9,25],[2,1,3]]))

btw. Can anyone tell me why numpy's np.cov() gives different result than 1./X.shape[0])*np.matmul(Xm.T,Xm) ???