memoryview ベクトルが与えられた関数があり、そのベクトルのノルムを計算したいと考えています。これまでは、memoryview を Numpy 配列に変換し、 を介してノルムを計算することでそれを達成しましたnp.sqrt(V.dot(V))。速度の理由からそのステップを取り除きたいのですが、次の実装のある時点でプログラムが失敗します。
cdef do_something(np.double_t[::1] M_mem):
cdef:
int i
np.double_t norm_mv = 0
np.double_t norm_np = 0
np.ndarray[np.double_t, ndim=1] V = np.copy(np.asarray(M_mem))
# Original implementation -- working
norm_np = np.sqrt(V.dot(V))
# My failed try with memoryview -- not working
for i in range(M_mem.shape[0]):
norm_mv += M_mem[i]**2
norm_mv = np.sqrt(norm_mv)
# norm_mv != norm_np
この理由は、十分に大きなベクトルの邪魔になる浮動小数点演算にあると思われます。Cython メモリビューのノルムを計算する数値的に安定した方法はありますか?
アップデート
チェックすると、丸め誤差はおそらく無意味であることがわかります。代わりに、本当に奇妙なことが起こっています。私の実際の機能は次のようになります。
@cython.boundscheck(False)
@cython.cdivision(True)
@cython.wraparound(False)
cdef np.double_t[:,::1] GS_coefficients(np.double_t[:,::1] M_mem):
cdef:
int n, i, k
int N_E = M_mem.shape[1]
np.ndarray[np.double_t, ndim=2] W = np.asarray(M_mem)
np.ndarray[np.double_t, ndim=2] V = np.copy(W)
np.double_t[:,::1] G = np.eye(N_E, dtype=np.float64)
np.longdouble_t norm = 0 # np.sqrt(V[:,0].dot(V[:,0]))
for i in range(M_mem.shape[0]):
norm += M_mem[i,0]**2
norm = sqrt(norm)
print("npx: ", np.sqrt(V[:,0].dot(V[:,0]))) # line 1
print("cp: ", norm) # line 2
V[:,0] /= norm
G[0,0] /= norm
for n in range(1, N_E):
for i in range(0, n):
G[n,i] = - (V[:,i].dot(W[:,n]))
V[:,n] += G[n,i] * V[:,i]
norm = np.sqrt(V[:,n].dot(V[:,n]))
V[:,n] /= norm
for i in range(n+1):
G[n,i] /= norm
return G
print結果が「どれだけ等しいか」をチェックするステートメントを挿入しましnormた。奇妙なことは、上記のコードが立っているように、すべてが正常に機能するようになったことです。しかし、最初の print ステートメント (1 行目) をコメントアウトすると、コードは関数を正常に実行しますが、プログラムの直後に失敗します。そこで何が起こっているのですか？これは、他の操作に影響を与えるべきではない単なるprintステートメントではありませんか?
更新 2
最小限の、完全で検証可能な例での私の試みは次のとおりです。
DEF N_E_cpt = 4
cimport cython
cimport numpy as np
import numpy as np
from libc.math cimport sqrt
@cython.boundscheck(False)
@cython.cdivision(True)
@cython.wraparound(False)
cdef np.double_t[:,::1] GS_coefficients(np.double_t[:,::1] M_mem):
"""Writes the coefficients, that the Gram-Schmidt procedure
provides in a Matrix and retruns it."""
cdef:
int n, i, k
int N_E = M_mem.shape[1]
np.ndarray[np.double_t, ndim=2] W = np.asarray(M_mem)
np.ndarray[np.double_t, ndim=2] V = np.copy(W)
np.double_t[:,::1] G = np.eye(N_E, dtype=np.float64)
np.longdouble_t norm = 0 # np.sqrt(V[:,0].dot(V[:,0]))
for i in range(M_mem.shape[0]):
norm += M_mem[i,0]**2
norm = sqrt(norm)
print("npx: ", np.sqrt(V[:,0].dot(V[:,0]))) # line 1
print("cp: ", norm) # line 2
V[:,0] /= norm
G[0,0] /= norm
for n in range(1, N_E):
for i in range(0, n):
G[n,i] = - (V[:,i].dot(W[:,n]))
V[:,n] += G[n,i] * V[:,i]
norm = np.sqrt(V[:,n].dot(V[:,n]))
V[:,n] /= norm
for i in range(n+1):
G[n,i] /= norm
return G
@cython.boundscheck(False)
@cython.cdivision(True)
@cython.wraparound(False)
cdef np.double_t[:,::1] G_mat(np.double_t[:,::1] M_mem):
"""Calls GS_coefficients and uses the coefficients to calculate
the entries of the transformation matrix G_ij"""
cdef:
np.double_t[:,::1] G_mem = GS_coefficients(M_mem)
int N_E = G_mem.shape[1]
np.double_t carr[N_E_cpt][N_E_cpt]
np.double_t[:,::1] G = carr
int n, i, j
# delete lower triangle in G
G[...] = G_mem
for i in range(N_E_cpt):
for j in range(0, i):
G[i,j] = 0.
for n in range(1, N_E):
for i in range(0, n):
for j in range(0, i+1):
G[n,j] += G_mem[n,i] * G[i,j]
return G
def run_test():
cdef:
np.double_t[:,::1] A_mem
np.double_t[:,::1] G
np.ndarray[np.double_t, ndim=2] A = np.random.rand(400**2, N)
int N = 4
A_mem = A
G = G_mat(A_mem)
X = np.zeros((400**2, N))
for i in range(0, N):
for j in range(0,i+1):
X[:,i] += G[i,j] * A[:,j]
print(X)
print("\n", X.T.dot(X))
run_test()
そのコードが何をするのかを理解する必要はないと思います。私にとっての謎は、そのprint声明がなぜ違いを生むのかということです。
したがって、このコードが想定しているのは、行列 A の列ベクトルとして記述された非正規直交ベクトル セットを取得し、次のようにベクトル セットを正規直交正規化する正規直交行列を返すことです。
したがって、A_{orthonormal} はコードの X 行列に相当します。直交正規行列の転置を直交正規行列自体で乗算すると、単一行列が得られます。これは、print# line1 ステートメントがそこにある限り得られるものです。それを削除するとすぐに、行列が直交していないことを意味する対角外のエントリも取得されます。なんで？