Быстрый метод нанесения синуса на матрицу

Question

Мне нужно оптимизировать этот код.

# K             sparse node coupling matrix
# P             vector of net power production at each node
# alpha         vector of damping at each node
# n             number of dimensions
# theta         vector of phases
# omega         vector of angular velocities


    def d_omega(K,P,alpha):
        def layer1(n):
            def layer2(theta,omega):
                T = theta[:, None] - theta
                T = np.sin(T)
                A = K*T
                A = A.dot(np.ones(n))
                return - alpha*omega + P - A
            return layer2
        return layer1

Я знаю, что большая часть времени вычислений используется в строке

T = np.sin(T)

. Есть ли более быстрый способ сделать это?

Answer 1

Используйте Формула Эйлера . Таким образом, мы можем уменьшить количество дорогостоящих вычислений функций с n^2 до 2n (если считать сложные функции в два раза). То, что мы делаем, выражаем синус внешнего различия как внешний продукт сложных экспонент. Поскольку умножение намного дешевле синуса, это приводит к хорошему net ускорению:

import numpy as np
from timeit import timeit

def pp():
    T = np.exp(1j*theta)
    T = np.outer(T,T.conj()).imag
    return T

def orig():
    T = theta[:, None] - theta
    T = np.sin(T)
    return T

rng = np.random.default_rng()

theta = rng.uniform(-np.pi,np.pi,1000)

print(np.allclose(orig(),pp()))
print(timeit(orig,number=10))
print(timeit(pp,number=10))

Пример выполнения:

True                     #  results are the same
0.2818465740419924       #  original method
0.04591922601684928      #  Euler

Ускорение (зависит от количества элементов тета, здесь 1000) ~ 6x

В зависимости от формы K (скаляр, 1D или 2D) мы можем оптимизировать далее:

def pp():
    T = np.exp(1j*theta)
    if K.ndim == 0:
        A = T*(T.sum().conj()*K)
    elif K.ndim == 1:
        A = T * (K@T.conj())
    else:
        A = T * (K@T.conj())
    return A.imag

def orig():
    T = theta[:, None] - theta
    T = np.sin(T)
    A = K*T
    A = A.dot(np.ones(1000))
    return A

rng = np.random.default_rng()

theta = rng.uniform(-np.pi,np.pi,1000)

for K in (rng.uniform(-10,10,shp) for shp in ((),(1000,),(1000,1000))):
    print(np.allclose(orig(),pp()))
    print(timeit(orig,number=10))
    print(timeit(pp,number=10))

Пример выполнения:

True
0.493153132032603
0.0012746050488203764
True
0.49636399815790355
0.0012419759295880795
True
0.47554834792390466
0.05685037490911782

Таким образом, с K скаляром или вектором мы получаем ускорение ~ 350x, в то время как K является матрицей, это всего ~ 8x.

Answer 2

Из задания

T0 = θ[:, None] - θ

мы имеем это T0[i,j] = θ[i]-θ[j], а из

T1 = sin(T0)

имеем T1[i,j] = sin(θ[i]-θ[j]), но из простой тригонометрии

T1[i,j] = sin(θ[i])*cos(θ[j]) - cos(θ[i])*sin(θ[j])

Признавая, что sin(θ[i])*cos(θ[j]) является элементом матрицы, которая является результатом внешнего произведения sin(θ) и cos(θ), мы можем написать небольшую программу, которую я интерактивно проверил

In [42]: from numpy import allclose, arange, cos, sin 
    ...: x = arange(5) 
    ...: s0 = sin(x[:,None]-x) 
    ...: s1 = sin(x[:,None]) * cos(x) 
    ...: s1 -= s1.T 
    ...: allclose(s0,s1)                                                                  
Out[42]: True