How to get orthogonal distances of vectors from plane in Numpy/Scipy?

Question 1

A more explicit way to achieve @Jaime's answer, you could be explicit about the construction of the projection operator:

def normalized(v):
    return v/np.sqrt( np.dot(v,v))
def ortho_proj_vec(v, uhat ):
    '''Returns the projection of the vector v  on the subspace
    orthogonal to uhat (which must be a unit vector) by subtracting off
    the appropriate multiple of uhat.
    i.e. dot( retval, uhat )==0
    '''
    return v-np.dot(v,uhat)*uhat

def ortho_proj_array( Varray, uhat ):
     ''' Compute the orhogonal projection for an entire array of vectors.
     @arg Varray:  is an array of vectors, each row is one vector
          (i.e. Varray.shape[1]==len(uhat)).
     @arg uhat: a unit vector
     @retval : an array (same shape as Varray), where each vector
               has had the component parallel to uhat removed.
               postcondition: np.dot( retval[i,:], uhat) ==0.0
               for all i. 
    ''' 
    return Varray-np.outer( np.dot( Varray, uhat), uhat )




# We need to ensure that the vectors defining the subspace are unit norm
v1hat=normalized( v1 )

# now to deal with v2, we need to project it into the subspace
# orhogonal to v1, and normalize it
v2hat=normalized( ortho_proj(v2, v1hat ) )
# TODO: if np.dot( normalized(v2), v1hat) ) is close to 1.0, we probably
# have an ill-conditioned system (the vectors are close to parallel)



# Act on each of your data vectors with the projection matrix,
# take the norm of the resulting vector.
result=np.zeros( M.shape[0], dtype=M.dtype )
for idx in xrange( M.shape[0] ):
    tmp=ortho_proj_vec( ortho_proj_vec(M[idx,:], v1hat), v2hat )             
    result[idx]=np.sqrt(np.dot(tmp,tmp))

 # The preceeding loop could be avoided via
 #tmp=orhto_proj_array( ortho_proj_array( M, v1hat), v2hat )
 #result=np.sum( tmp**2, axis=-1 )
 # but this results in many copies of matrices that are the same 
 # size as M, so, initially, I prefer the loop approach just on
 # a memory usage basis.

This is really just the generalization of Gram-Schmidt orthogonalization procedure. note that at the end of this process we have np.dot(v1hat, v1hat.T)==1, np.dot(v2hat,v2hat.T)==1, np.dot(v1hat, v2hat)==0 (to within numerical precision)

Question 2

You need to construct the unit-normal to the plane:

In three dimensions this can be easily done:

nhat=np.cross( v1, v2 )
nhat=nhat/np.sqrt( np.dot( nhat,nhat) )

and then dot this with each of your vectors; which I assume is an Nx3 matrix M

result=np.zeros( M.shape[0], dtype=M.dtype )
for idx in xrange( M.shape[0] ):
    result[idx]=np.abs(np.dot( nhat, M[idx,:] ))

so that result[idx] is the distance of the idx'th vector from the plane.

Question 3

EDIT The original code I wrote does not work properly, so I have removed it. But following the same idea, explained below, if you spend some time thinking at it, there is no need for Cramer's rule, and the code can be streamlined quite a bit as follows :

def distance(v1, v2, u) :
    u = np.array(u, ndmin=2)
    v = np.vstack((v1, v2))
    vv = np.dot(v, v.T) # shape (2, 2)
    uv = np.dot(u, v.T) # shape (n ,2)
    ab = np.dot(np.linalg.inv(vv), uv.T) # shape(2, n)
    w = u - np.dot(ab.T, v)
    return np.sqrt(np.sum(w**2, axis=1)) # shape (n,)

To make sure it works properly, I have packed Dave's code into a function as distance_3d and tried the following:

>>> d, n = 3, 1000
>>> v1, v2, u = np.random.rand(d), np.random.rand(d), np.random.rand(n, d)
>>> np.testing.assert_almost_equal(distance_3d(v1, v2, u), distance(v1, v2, u))

But of course it now works for any d:

>>> d, n = 1000, 3
>>> v1, v2, u = np.random.rand(d), np.random.rand(d), np.random.rand(n, d)
>>> distance(v1, v2, u)
array([ 10.57891286,  10.89765779,  10.75935644])

You have to decompose your vector, lets call it u, in the sum of two vectors, u = v + w, v is in the plane, and so can be decomposed as v = a * v1 + b * v2, while w is perpendicular to the plane, and thus np.dot(w, v1) = np.dot(w, v2) = 0.

If you write u = a * v1 + b * v2 + w and take the dot product of this expression with v1 and v2, you get two equations with two unknowns:

np.dot(u, v1) = a * np.dot(v1, v1) + b * np.dot(v2, v1)
np.dot(u, v2) = a * np.dot(v1, v2) + b * np.dot(v2, v2)

Since it is only a 2x2 system, we can solve it using Cramer's rule as:

uv1 = np.dot(u, v1)
uv2 = np.dot(u, v2)
v11 = np.dot(v1, v2)
v22 = np.dot(v2, v2)
v12 = np.dot(v1, v2)
v21 = np.dot(v2, v1)
det = v11 * v22 - v21 * v12
a = (uv1 * v22 - v21 * uv2) / det
b = (v11 * uv2 - uv1 * v12) / det

From here, you can get:

w = u - v = u - a * v1 - b * v2

and the distance to the plane is the modulus of w.