python - 将正弦应用于矩阵的快速方法

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

使用Euler's formula 。这样我们就可以将昂贵的函数求值次数从 n^2 减少到 2n(如果我们将复杂函数加倍)。我们正在做的是将外差的正弦表示为复指数的外积。由于乘法比正弦运算便宜得多，因此可以实现良好的网络加速:

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

加速比(取决于 theta 的元素数量，此处为 1000)~6 倍

根据 K 的形状(标量、一维或二维)，我们可以进一步优化:

def pp():
    T = np.exp(1j*theta)
    if K.ndim == 0:
        A = T*(T.sum().conj()*K)
    elif K.ndim == 1:
        A = T * (<a href="https://stackoverflow.com/cdn-cgi/l/email-protection" class="__cf_email__" data-cfemail="8bc0cbdfa5e8e4e5e1" rel="noreferrer noopener nofollow">[email protected]</a>())
    else:
        A = T * (<a href="https://stackoverflow.com/cdn-cgi/l/email-protection" class="__cf_email__" data-cfemail="a7ece7f389c4c8c9cd" rel="noreferrer noopener nofollow">[email protected]</a>())
    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 标量或向量，我们可以获得约 350 倍的加速，而如果 K 是矩阵，则仅为约 8 倍。