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

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

优化，第 1 轮:

避免重新计算 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]

将时间减少一半

优化，第 2 轮:

向量化 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]

优化，第 3 轮:

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

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

绕圈 1.2 秒 Ns=1000

编辑:

第 4 轮:

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 ，我已经把它利用到了极致!

耗时:400ms , 对于 Ns=1000