我有一组向量作为 numpy 数组。我需要从由 2 个向量 v1 和 v2 定义的平面中获取每个平面的正交距离。我可以使用 Gram-Schmidt 过程轻松获得单个向量。有没有一种方法可以非常快速地处理许多向量,而无需遍历每个向量或使用 np.vectorize?
谢谢!
我有一组向量作为 numpy 数组。我需要从由 2 个向量 v1 和 v2 定义的平面中获取每个平面的正交距离。我可以使用 Gram-Schmidt 过程轻松获得单个向量。有没有一种方法可以非常快速地处理许多向量,而无需遍历每个向量或使用 np.vectorize?
谢谢!
您需要构造平面的单位法线:
在三个维度上,这很容易做到:
nhat=np.cross( v1, v2 )
nhat=nhat/np.sqrt( np.dot( nhat,nhat) )
然后用你的每个向量点这个;我假设是一个Nx3
矩阵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,:] ))
所以这result[idx]
就是idx'th
向量到平面的距离。
实现@Jaime 答案的更明确的方法是,您可以明确说明投影运算符的构造:
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.
这实际上只是对 Gram-Schmidt 正交化过程的概括。请注意,在此过程结束时,我们有np.dot(v1hat, v1hat.T)==1
, np.dot(v2hat,v2hat.T)==1
, np.dot(v1hat, v2hat)==0
(在数值精度范围内)
编辑我写的原始代码不能正常工作,所以我把它删除了。但是遵循相同的想法,下面解释,如果你花一些时间思考它,就不需要克莱默规则,代码可以简化如下:
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,)
为了确保它正常工作,我将 Dave 的代码打包到一个函数中,distance_3d
并尝试了以下操作:
>>> 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))
但当然它现在适用于任何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])
你必须分解你的向量,我们称之为u
,在两个向量之和中,u = v + w
,v
在平面内,因此可以分解为v = a * v1 + b * v2
,而w
垂直于平面,因此np.dot(w, v1) = np.dot(w, v2) = 0
。
如果你用 和 写并u = a * v1 + b * v2 + w
取这个表达式的点积,你会得到两个有两个未知数的方程:v1
v2
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)
由于它只是一个 2x2 系统,我们可以使用Cramer 规则解决它:
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
从这里,您可以获得:
w = u - v = u - a * v1 - b * v2
到平面的距离是 的模数w
。