У меня есть функция, которой присваивается вектор обзора памяти, и я хочу вычислить норму этого вектора.До сих пор я добивался этого путем преобразования обзора памяти в массив 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.Странность в том, что теперь все работает нормально, как показано выше.Но когда я закомментирую первый оператор печати (строка 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.Как только вы удалите его, вы также получите недиагональные записи, что означает, что матрица даже не ортогональна.Почему?