减少输入矩阵系数所需的时间

Question

使用 Python，我想构建一个方阵，其系数是在某些给定点评估的函数。矩阵是下三角矩阵，但在大小为 1000 的情况下，输入所有系数需要 200 多秒。这让我很惊讶，因为我已经处理过大小大于 100 万的方阵。

我猜想时间的损失来自定义矩阵所有非零系数的函数：它是四个项的乘积，每个项都涉及一个复幂（请参见下面的最小示例）。 如果有人可以帮助我，我会非常高兴：目的是在几秒钟内达到 100 万的规模（希望这是可行的......）。也许使用 ** 对于复数幂不是最优的；或者也许有一种更快的方法来填充下三角矩阵？任何帮助将不胜感激！

对于下面的代码，Ns 是方阵的大小。

import math
from math import *
import numpy as np
from numpy import exp, log, sqrt, cos, sin, arange
import time


tmps1 = time.time()


Ns = 100
h = 1/float(Ns)
Ns = int(Ns)

mu = 0.01

rn = -6.34
rc = 0.86
r1 = 1.32
r2 = 4.16

P = np.zeros((Ns + 1, 1))
for j in range(0, Ns + 1):
    P[j] = r1 + j*(r2-r1)*h

kappan = 0.24
kappac = -0.24
kappa1 = 0.095
kappa2 = -0.095


z = complex(0.01,0.01)


def E(exponent,p):
    return (  ( ((p-rn)/(2-rn))**((exponent)/(2*kappan)) )*( ((p-rc)/(2-rc))**((exponent)/(2*kappac)) )*( ((p-r1)/(2-r1))**((exponent)/(2*kappa1)) )*( ((r2-p)/(r2-2))**((exponent)/(2*kappa2)) )  )

def D(p,r):
    return (      ( 1/(2*complex(0.0,1.0)*(z-mu/r1)) )*(   ( 1/(2*(kappa1-complex(0.0,1.0)*(z-mu/r1))) )*( E(2*(kappa1-complex(0.0,1.0)*(z-mu/r1)),r)*E(2*complex(0.0,1.0)*(z-mu/r1),p) )   -   ( 1/(2*kappa1) )*E(2*kappa1,r)   )      )


A = np.zeros((Ns-1, Ns-1), dtype=np.complex_)

for j in range(1, Ns-1):
    for k in range(1, j+1):
        A[j,k] = D(P[j+1], P[k+1]) - D(P[j+1], P[k])


tmps2 = time.time()-tmps1
print "\n\nExecution time = %f\n\n" %tmps2

Answer 1

另一种节省多次计算函数值的方法是memoization 装饰器。这个技巧将为您带来相同的结果，而无需手动推出您的代码。

def memoize(f):
    results = {}
    def helper(x,y):
        key = (float(x), float(y))
        if key not in results:
            results[key] = f(x,y)
        return results[key]
    return helper

@memoize
def D(p,r):
    return (      ( 1/(2*C01*(z-mu/r1)) )*(   ( 1/(2*(kappa1-C01*(z-mu/r1))) )*( E(2*(kappa1-C01*(z-mu/r1)),r)*E(2*C01*(z-mu/r1),p) )   -   ( 1/(2*kappa1) )*E(2*kappa1,r)   )      )

Answer 2

除了ckedar之前的回答，可以通过使用numpy向量化改进P向量的计算来加快一点速度。在以下代码中，我使用了 numpy 的 fromfunction 函数：

const = (r2-r1)*h
P = np.fromfunction(lambda j, i: r1 + j*const, (Ns + 1, 1), dtype=float)

在我的机器上使用此代码，Ns = 1000，生成 P 大约需要 15.2 µs，而使用您的代码生成 P 需要 421 µs（速度提高了 96%）。

Answer 3

优化，第一轮：

避免重新计算 D。

d = np.zeros((Ns+1, Ns+1), dtype=np.complex_)
for j in range(2, Ns):
    for k in range(1, j+1):
        d[j,k] = D(P[j], P[k])

A = np.zeros((Ns-1, Ns-1), dtype=np.complex_)
for j in range(1, Ns-1):
    for k in range(1, j+1):
        A[j,k] = d[j+1, k+1] - d[j+1, k]

时间减半

优化，第二轮：

向量化 D 的计算。
将 d 的计算替换为：

d = np.zeros((Ns, Ns), dtype=np.complex_)
for shift in range(0, Ns-1):
    x = D(P, np.roll(P,shift))
    for j in range(shift+1, Ns):
        d[j,j-shift] = x[j]

优化，第三轮：

向量化 A 的计算。
将 A 的计算替换为：

d = np.roll(d, -1, axis=(0,1))
A = d - np.roll(d, 1, axis=1)
A *= np.tri(*A.shape)

Ns=1000 大约需要 1.2 秒

编辑：

第四轮：

P = np.fromfunction(lambda j, i: r1 + j*(r2-r1)*h, (Ns+1, 1), dtype=float)
Q = np.fromfunction(lambda j, i: r1 + ((j-i)%(Ns+1))*(r2-r1)*h, (Ns+1, Ns), dtype=float)
# P = Q[:,0].reshape(11,1)
d = D(P, Q)
d = np.fromfunction(lambda r, c: d[r,(r-c)], (Ns,Ns), dtype=int)

d = np.roll(d, -1, axis=(0,1))
A = d - np.roll(d, 1, axis=1)
A = np.tril(A)

感谢 Fabio Lipreri 对 fromfunction 的建议，我已经充分利用了它！

耗时：400 毫秒，Ns=1000