ベクトルのセットを numpy 配列として持っています。2 つのベクトル v1 と v2 で定義される平面からそれぞれの直交距離を取得する必要があります。Gram-Schmidt プロセスを使用して、単一のベクトルに対してこれを簡単に取得できます。それらのそれぞれをループしたり、np.vectorize を使用したりせずに、多くのベクトルに対して非常に高速に実行する方法はありますか?
ありがとう！
質問する
6485 次
3 に答える
4
平面に垂直な単位を構築する必要があります。
3 次元では、これは簡単に実行できます。
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
于 2013-01-31T18:45:09.870 に答える
4
@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.
これは、グラム・シュミット直交化手順の一般化にすぎません。np.dot(v1hat, v1hat.T)==1このプロセスの最後に、 、np.dot(v2hat,v2hat.T)==1、np.dot(v1hat, v2hat)==0(数値精度内)があることに注意してください。
于 2013-02-01T21:51:51.073 に答える
2
編集私が書いた元のコードは正しく動作しないため、削除しました。しかし、以下で説明する同じ考え方に従って、時間をかけて考えれば、Cramer の規則は必要なく、コードは次のようにかなり合理化できます。
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 のコードを as 関数にパック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ベクトルを分解する必要があります。2 つのベクトルの合計でu = v + wと呼びましょう。は平面にあるため、vとして分解できます。v = a * v1 + b * v2wnp.dot(w, v1) = np.dot(w, v2) = 0
u = a * v1 + b * v2 + wと を使用してこの式のドット積を作成して取得するv1とv2、2 つの未知数を含む 2 つの方程式が得られます。
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。
于 2013-02-01T06:21:44.493 に答える