Используйте NumPy с умом, чтобы ускорить этот цикл

Question

Это вполне теоретический вопрос. Речь идет о l oop в конце прикрепленного кода.

У меня есть три совместимых массива: Z, Wt_1 и Wt_2.

Я использовал векторизованные столбцы для ускорения l oop, ie Я кодирую столбцы, как если бы они были скалярами, и весь процесс работает. Но этот l oop по-прежнему занимает вечность (2 мин +). Здесь N = 5000, а длина столбца - 100000.

Есть ли синтаксис, который ускоряет это?

M=100000

N=5000
t0=0.132
T=0.22
dt=T / N

Z0=8
Y0=7

ssqZ=0.04
ssqY=0.01
sigma_Z=np.sqrt(ssqZ)
sigma_Y=np.sqrt(ssqY)
rho=- 0.4

print("Genreating the Brownian Motions...")
Wt_1=rd.randn(M//2,N)
# shape is (50 000,5 000)
Wt_1=np.vstack((Wt_1,-Wt_1))
#print(np.shape(Wt_1))
Wt_2=rd.randn(M//2,N)
Wt_2=np.vstack((Wt_2,-Wt_2))
print("BM generated!\n")

Z=np.zeros((M,N+1))
# shape is (100 000, 5 001)
Y=np.zeros((M,N+1))
#in SciPy indexing starts from 0 :
Z[:,0]=Z0
Y[:,0]=Y0

for i in range(1,N+1) :
    Z[:,i]=Z[:,i-1] - 0.5*(sigma_Z**2)*dt + sigma_Z*np.sqrt(dt)*Wt_1[:,i-1]

    Y[:,i]=Y[:,i-1] + sigma_Y*np.sqrt(dt)* (rho*Wt_1[:,i-1] + np.sqrt(1 - rho**2)*Wt_2[:,i-1])

Я должен добавить, что та же процедура в Matlab занимает 19 секунд, тогда как этот Python нужно 137 с, что мне противоречит

Answer 1

Эту задачу довольно сложно эффективно векторизовать Numpy. В вашей версии у вас есть векторизованная команда, которая выполняет итерацию по значениям, которые не хранятся в памяти непрерывно Numpy / C мажор строки, Fortran / Matlab / Julia - мажор столбца.

Если вы хотите использовать Numpy, вероятно, самый простой способ - это транспонировать все массивы и операции или использовать массивы, упорядоченные по Фортрану (должен давать коэффициент 4).

Другой способ - записать циклы и сделать использование Numba Numba или Cython .

Numba Пример

@nb.njit(fastmath=True,parallel=True,cache=True)
def Numba(Z0,Y0,sigma_Y,sigma_Z,Wt_1,Wt_2,rho,dt):
    N,M=Wt_1.shape

    Z=np.empty((M,N+1))
    Y=np.empty((M,N+1))
    Z[:,0]=Z0
    Y[:,0]=Y0
    #the possibillity of negative values can lead to a slow wrap-around check
    for i in nb.prange(M):
        #only beneficial for very large 
        for j in range(N):
            Z[i,j+1]=Z[i,j] - 0.5*(sigma_Z**2)*dt + sigma_Z*np.sqrt(dt)*Wt_1[i,j]
            Y[i,j+1]=Y[i,j] + sigma_Y*np.sqrt(dt)* (rho*Wt_1[i,j] + np.sqrt(1 - rho**2)*Wt_2[i,j])
    return Z,Y

Ваш Numpy пример

Также было бы неплохо добавить пример того, как вы получили свои тайминги, поскольку генерация входных данных занимает больше времени, чем фактическое вычисление (по крайней мере, для версии Numba).

def Numpy_not_aligned(Z0,Y0,sigma_Y,sigma_Z,Wt_1,Wt_2,rho,dt):
    Z=np.zeros((M,N+1))
    # shape is (100 000, 5 001)
    Y=np.zeros((M,N+1))
    #in SciPy indexing starts from 0 :
    Z[:,0]=Z0
    Y[:,0]=Y0

    for i in range(1,N+1) :
        Z[:,i]=Z[:,i-1] - 0.5*(sigma_Z**2)*dt + sigma_Z*np.sqrt(dt)*Wt_1[:,i-1]
        Y[:,i]=Y[:,i-1] + sigma_Y*np.sqrt(dt)* (rho*Wt_1[:,i-1] + np.sqrt(1 - rho**2)*Wt_2[:,i-1])
    return Y,Z

Время

# M=M=50000
%timeit Numpy_not_aligned(Z0,Y0,sigma_Y,sigma_Z,Wt_1,Wt_2,rho,dt)
#25.9 s ± 82.1 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)    
%timeit Numba(Z0,Y0,sigma_Y,sigma_Z,Wt_1,Wt_2,rho,dt)
#1.03 s ± 127 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Answer 2

Попробуем вывести несколько повторяющихся вычислений за пределы l oop

sub_1 = 0.5*(sigma_Z**2)*dt
sub_2 = sigma_Z*np.sqrt(dt)
sub_3 = np.sqrt(1 - rho**2)
sub_4 = sigma_Y*np.sqrt(dt)

for i in range(1,N+1) :
    Z[:,i]=np.add(np.subtract(Z[:,i-1] ,sub_1), np.multiply(Wt_1[:,i-1], sub_2))

    Y[:,i]=np.add(Y[:,i-1] , np.multiply(np.add(np.multiply(Wt_2[:,i-1],sub_3),np.multiply(Wt_1[:,i-1],rho)), sub_4))

Answer 3

Вы не получите никаких реальных улучшений скорости, если не устраните l oop. Естественно, если мы присмотримся внимательнее, мы сможем сформулировать два уравнения по-разному. Я кратко продемонстрирую это с помощью первого, Z:

Z(i) = Z(i-1) - C + D*Wt_1(i-1), where C and D are constants.
Z(1) = Z(0) - C + D*Wt_1(0)
Z(2) = Z(1) - C + D*Wt_1(1) = Z(0) - 2*C + D*(Wt_1(0)+Wt_1(1))
i.e. Z(i) = Z(0) -i*C + D*Sum(Wt_1(0..i-1))

Как вы можете видеть, член Z (0) на самом деле является константой и присутствует во всех столбцах; мы можем просто установить его во время инициализации. Единственное реальное повторение - в Wt_1. В уравнении Y, аналогично, повторение зависит от Wt_1 и Wt_2.

И здесь у нас есть лазейка: это повторение W фактически представляет собой растущую по столбцам или кумулятивную сумму, позволяющую нам отбросить for l oop и использовать вместо него np.cumsum. Мы также можем транслировать и добавлять константу C в уравнение Z, умножая ее на np.arange(N+1).

Z = np.ones((M,N+1))*Z0
Y = np.ones((M,N+1))*Y0
W1 = np.cumsum(Wt_1,axis=1)
W2 = np.cumsum(Wt_2,axis=1)

# Calculate Z
Z -= (0.5*(sigma_Z**2)*dt)*np.arange(N+1)
Z[:,1:] = Z[:,1:] + W1*sigma_Z*np.sqrt(dt)

# Calculate Y
Y[:,1:] = Y[:,1:] + sigma_Y*np.sqrt(dt)*(rho*W1 + np.sqrt(1 - rho**2)*W2)

Answer 4

Я попытался переместить вычисление из l oop. Кажется, что нижнее выражение представляет собой простое вычитание матрицы, поэтому я полностью исключил его.

sqrt_dt = np.sqrt(dt)
A = 0.5*(sigma_Z**2)*dt
B = sigma_Z* sqrt_dt
C = sigma_Y * sqrt_dt * (rho*Wt_1 + np.sqrt(1 - rho**2)*Wt_2)

Y[:,1:N+1] = Y[:,0:N] - C[:,0:N]

for i in range(1,N+1) :
    Z[:,i]=Z[:,i-1] - A + B  * Wt_1[:,i-1]

Надеюсь, это поможет!