6

I have a set of vectors as a numpy array. I need to get the orthogonal distances of each of those from a plane defined by 2 vectors v1 and v2. I can get this easily for a single vector using the Gram-Schmidt process. Is there a way to do it very fast in for many vectors without having to loop through each of them, or using np.vectorize?

Thanks!

costaz
  • 63
  • 1
  • 4

3 Answers3

4

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.

Jaime
  • 65,696
  • 17
  • 124
  • 159
Dave
  • 7,555
  • 8
  • 46
  • 88
  • 6
    You don't need the loop. `result = np.abs(M.dot(nhat))` will work. – Warren Weckesser Jan 31 '13 at 20:09
  • 1
    Thanks for the response folks! However, my vectors are more than 3 dimensions (50 dimensions). I run into the error - ValueError: incompatible dimensions for cross product (dimension must be 2 or 3). What do i do? – costaz Feb 01 '13 at 00:43
  • I infered 3-D from you question since you did not specify the dimensionality of the problem explicitly. – Dave Feb 01 '13 at 14:32
4

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)

Dave
  • 7,555
  • 8
  • 46
  • 88
  • Nice! For a large number of vectors you would want to make sure that your function can operate on arrays of row (or column) vectors, so that you can skip the python loop. – Jaime Feb 02 '13 at 15:29
  • @Jaime I don't usually think about how to vectorize things, so if you see a way to do it, feel free to edit this. – Dave Feb 02 '13 at 18:12
2

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.

Jaime
  • 65,696
  • 17
  • 124
  • 159