Ускорение цикла

Question

Очень медленная скорость 3d цикла: 43 с. Можно ли заметно ускорить такой код?

import numpy as np
import scipy.constants
import time

start_time = time.time()

c = scipy.constants.c #speed of light

N = 1024
Nx = 128
Ny = 128
x_size = 50
y_size = 50
U = np.zeros(shape=(N, Nx, Ny), dtype=complex)

fx=np.linspace( - Nx / ( 2 * x_size ), Nx / ( 2 * x_size ), Nx )
fy=np.linspace( - Ny / ( 2 * y_size ), Ny / ( 2 * y_size ), Ny )
nu=np.linspace( - x_size, x_size, N )

for k in range( N ):
    for i in range( Nx ):
        for j in range( Ny ):
            U[k, i, j] = 1 - ((c * fx[i] / nu[k]) ** 2) - ((c * fy[j] / nu[k]) ** 2 )

print("--- %s seconds ---" % (time.time() - start_time))

Answer 1

Ускорение в 430 раз

Задача поставлена некорректно

Мне удалось добиться решения за 0.07 секунды с помощью оптимизации вычислений и алгоритмизации.

Сразу следует отметить, что данный процесс не сходится и для различных оптимизаций мы можем получить очень сильные вариации. Это связано с неустойчивостью численного решения, который приводит автор. В частности, для конечного решения это означает, что математические действия над данным выражениям могут носить фатальный характер и конечный результат будет очень сильно меняться. Пример детальный разбор представлен последним пунктом.

1

Во-первых, перенесите начало отсчёта на момент "перед циклом" (в вашем случае, это ничего не изменит, конечно, но, всё-таки)

import numpy as np
import scipy.constants
import time

c = scipy.constants.c #speed of light

N = 1024
Nx = 128
Ny = 128
x_size = 50
y_size = 50
U = np.zeros(shape=(N, Nx, Ny), dtype=complex)

fx=np.linspace( - Nx / ( 2 * x_size ), Nx / ( 2 * x_size ), Nx )
fy=np.linspace( - Ny / ( 2 * y_size ), Ny / ( 2 * y_size ), Ny )
nu=np.linspace( - x_size, x_size, N )


start_time = time.time()
for k in range( N ):
    for i in range( Nx ):
        for j in range( Ny ):
            U[k, i, j] = 1 - ((c * fx[i] / nu[k]) ** 2) - ((c * fy[j] / nu[k]) ** 2 )

print("--- %s seconds ---" % (time.time() - start_time))

2

Всё сильно зависит от вашего процессора. Например, у меня получается 30 секунд.

3

1024 * 128 * 128 -- это не то чтобы очень мало. 16 млн. Кажется, что 30 секунд вполне разумная цена. Более того, внутри цикла, у вас есть довольно тяжёлые операции

4

Если хотите ускорения, то давайте оптимизируем цикл:

c2 = c ** 2
for k in range( N ):
    nu_k2 = nu[k] ** 2
    for i in range( Nx ):
        fx_i = fx[i] ** 2
        for j in range( Ny ):
            U[k, i, j] = 1 - c2 / nu_k2 * (fx_i + ((fy[j]) ** 2 ))

Основная оптимизация происходит за счёт того, что мы не будем повторно вычислять значения. Таким образом, после того, как мы вынесли nu_k2 и fx_i, получаем 18 секунд.

Откажемся от деления на каждом шаге:

c2 = c ** 2
for k in range( N ):
    nu_k2 = nu[k] ** 2
    d = c2 / nu_k2
    for i in range( Nx ):
        fx_i = fx[i] ** 2
        for j in range( Ny ):
            U[k, i, j] = 1 - d * (fx_i + (fy[j]) ** 2 )

и получим прирост ещё в 1-2 секунды.

Так как numpy написан на C, то можем сделать следующее:

fx=np.linspace( - Nx / ( 2 * x_size ), Nx / ( 2 * x_size ), Nx ) ** 2
fy=np.linspace( - Ny / ( 2 * y_size ), Ny / ( 2 * y_size ), Ny ) ** 2
nu=np.linspace( - x_size, x_size, N )

c2 = c ** 2
for k in range( N ):
    nu_k2 = nu[k] ** 2
    d = c2 / nu_k2
    for i in range( Nx ):
        for j in range( Ny ):
            U[k, i, j] = 1 - d * (fx[i] + fy[j])

Сократим расходы до 13 секунд.

6 (0.07 секунды)

Финальное решение. Вынесем общую подматрицу в предподсчёт и будем её переиспользовать. И, наконец, 0.07 секунды

a = np.zeros((fx.shape[0], fy.shape[0]))
for i in range( Nx ):
    for j in range( Ny ):
        a[i, j] = fx[i] + fy[j]

c2 = c ** 2
for k in range( N ):
    nu_k2 = nu[k] ** 2
    d = c2 / nu_k2
    U[k, :, :] = 1 - d * a

7

Запишем проверку, которая покажет, что решения идентичные

import numpy as np
import scipy.constants

c = scipy.constants.c #speed of light

N = 8
Nx = 4
Ny = 4
x_size = 50
y_size = 50
U = np.zeros(shape=(N, Nx, Ny), dtype=complex)

fx=np.linspace( - Nx / ( 2 * x_size ), Nx / ( 2 * x_size ), Nx )
fy=np.linspace( - Ny / ( 2 * y_size ), Ny / ( 2 * y_size ), Ny )
nu=np.linspace( - x_size, x_size, N )


for k in range( N ):
    for i in range( Nx ):
        for j in range( Ny ):
            U[k, i, j] = 1 - ((c * fx[i] / nu[k]) ** 2) - ((c * fy[j] / nu[k]) ** 2 )

U0 = np.zeros(shape=(N, Nx, Ny), dtype=complex)

c2 = c ** 2
fx=np.linspace( - Nx / ( 2 * x_size ), Nx / ( 2 * x_size ), Nx ) ** 2
fy=np.linspace( - Ny / ( 2 * y_size ), Ny / ( 2 * y_size ), Ny ) ** 2
nu=c2 / np.linspace( - x_size, x_size, N ) ** 2

a = np.zeros((fx.shape[0], fy.shape[0]))
for i in range( Nx ):
    for j in range( Ny ):
        a[i, j] = fx[i] + fy[j]


for k in range(N):
    for i in range( Nx ):
        for j in range( Ny ):
            U0[k, :, :] = 1 - nu[k] * a

print(abs(U0 - U).sum())

Норма равно 0.01

8

Далее, более формально можно привести к матричному виду и ещё ускорить. Но это уже совсем другая задача.

Анализ сходимости

Рассмотрим небольшие начальные данные. Сделаем простое преобразование. В результирующей формуле вынесем минус за скобки. Тогда получим выражение:

1 - ((c * fx[i] / nu[k]) ** 2) + ((c * fy[j] / nu[k]) ** 2 )

Итоговый код, который сравивает получающиеся решения:

import numpy as np
import scipy.constants
import time

c = scipy.constants.c #speed of light

N = 64
Nx = 32
Ny = 32
x_size = 50
y_size = 50
U = np.zeros(shape=(N, Nx, Ny), dtype=complex)

fx=np.linspace( - Nx / ( 2 * x_size ), Nx / ( 2 * x_size ), Nx )
fy=np.linspace( - Ny / ( 2 * y_size ), Ny / ( 2 * y_size ), Ny )
nu=np.linspace( - x_size, x_size, N )


for k in range( N ):
    for i in range( Nx ):
        for j in range( Ny ):
            U[k, i, j] = 1 - ((c * fx[i] / nu[k]) ** 2) - ((c * fy[j] / nu[k]) ** 2 )

U0 = np.zeros(shape=(N, Nx, Ny), dtype=complex)

fx=np.linspace( - Nx / ( 2 * x_size ), Nx / ( 2 * x_size ), Nx )
fy=np.linspace( - Ny / ( 2 * y_size ), Ny / ( 2 * y_size ), Ny )
nu=np.linspace( - x_size, x_size, N )


for k in range( N ):
    for i in range( Nx ):
        for j in range( Ny ):
            x = ((c * fx[i] / nu[k]) ** 2) + ((c * fy[j] / nu[k]) ** 2 )
            U0[k, i, j] = 1 - x

print(abs(U0 - U).sum())

Если мы посмотрим на результат, то заметим, что исходное решение и предложенное отличается на 1058. Таким образом, делаем вывод, что данный процесс не сходится. И здесь вопрос к автору. Что это за задача

Answer 2

Зачем тип complex?? Не хочешь перейти на float? Тгда вычисления будут побыстрее )
Оптимизируй цикл. Попробуй так:

    import numpy as np
    import scipy.constants
    import time

    start_time = time.time()
    
    c = scipy.constants.c #speed of light
    
    N = 1024
    Nx = 128
    Ny = 128
    x_size = 50
    y_size = 50
    U = np.zeros(shape=(N, Nx, Ny), dtype=complex)
    
    c_fx_i_div_nu_k = np.zeros(shape=(N, Nx), dtype=complex)
    c_fy_j_div_nu_k = np.zeros(shape=(N, Ny), dtype=complex)
    
    fx=np.linspace( - Nx / ( 2 * x_size ), Nx / ( 2 * x_size ), Nx )
    fy=np.linspace( - Ny / ( 2 * y_size ), Ny / ( 2 * y_size ), Ny )
    nu=np.linspace( - x_size, x_size, N )
    
    c_div_nu = c / nu
    
    for k in range( N ):
        for i in range( Nx ):
            c_fx_i_div_nu_k[k,i] = (c_div_nu[k] * fx[i])**2
            
    for k in range( N ):
        for j in range( Ny ):
            c_fy_j_div_nu_k[k,j] = (c_div_nu[k] * fy[j])**2        
           
    for k in range( N ):
        for i in range( Nx ):
            for j in range( Ny ):
                U[k, i, j] = 1 - c_fx_i_div_nu_k[k,i] - c_fy_j_div_nu_k[k,j]
    
    print("--- %s seconds ---" % (time.time() - start_time))

так более чем в 2 раза быстрее.
а насколько быстрее хочется?
в 10 раз?

тогда надо подумать о вычислениях на видеокарте или при помощи, многопоточности numba
Вот пример с float + numba

import numpy as np
import scipy.constants
import time
import numba 

start_time = time.time()

@numba.njit
def calc():
    c = scipy.constants.c #speed of light

    N = 1024
    Nx = 128
    Ny = 128
    x_size = 50
    y_size = 50
    UU = np.zeros(shape=(N, Nx, Ny), dtype=np.float64)

    c_fx_i_div_nu_k = np.zeros(shape=(N, Nx), dtype=np.float64)
    c_fy_j_div_nu_k = np.zeros(shape=(N, Ny), dtype=np.float64)

    fx=np.linspace( - Nx / ( 2 * x_size ), Nx / ( 2 * x_size ), Nx )
    fy=np.linspace( - Ny / ( 2 * y_size ), Ny / ( 2 * y_size ), Ny )
    nu=np.linspace( - x_size, x_size, N )

    c_div_nu = c / nu

    for k in range( N ):
        for i in range( Nx ):
            c_fx_i_div_nu_k[k,i] = (c_div_nu[k] * fx[i])**2

    for k in range( N ):
        for j in range( Ny ):
            c_fy_j_div_nu_k[k,j] = (c_div_nu[k] * fy[j])**2        

    for k in range( N ):
        for i in range( Nx ):
            for j in range( Ny ):
                UU[k, i, j] = 1 - c_fx_i_div_nu_k[k,i] - c_fy_j_div_nu_k[k,j]
    return UU

U = calc()
print("--- %s seconds ---" % (time.time() - start_time))
print(U[3, 3, 3], U[-1, -1, -1])

но - проверь сам. Правильно ли посчиталось... ))
numba - она такая...

На моем стареньком компьюторе получилось:

--- 1.2685229778289795 seconds ---

Что в 85 раз быстрее исходного варианта

Answer 3

Я просто добавил использование Numba. Она очень любит циклы в циклах в циклах. Результат не проверял.

import numpy as np
import numba as nb
import scipy.constants
import time

@nb.njit
def calc_loops():
    c = scipy.constants.c #speed of light

    N = 1024
    Nx = 128
    Ny = 128
    x_size = 50
    y_size = 50
    U = np.zeros(shape=(N, Nx, Ny), dtype=nb.types.complex128)

    fx=np.linspace( - Nx / ( 2 * x_size ), Nx / ( 2 * x_size ), Nx )
    fy=np.linspace( - Ny / ( 2 * y_size ), Ny / ( 2 * y_size ), Ny )
    nu=np.linspace( - x_size, x_size, N )

    for k in range( N ):
        for i in range( Nx ):
            for j in range( Ny ):
                U[k, i, j] = 1 - ((c * fx[i] / nu[k]) ** 2) - ((c * fy[j] / nu[k]) ** 2 )
    return U

start_time = time.time()
U = calc_loops()
print("--- %s seconds ---" % (time.time() - start_time))

Вывод:

--- 0.4318578243255615 seconds ---

БЛОГ НА HUSL

Ускорение цикла

Ответы (3 шт):

Ускорение в 430 раз

1

2

3

4

6 (0.07 секунды)

7

8

Анализ сходимости