使用 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
最佳答案
优化,第 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
关于python - 减少输入矩阵系数所需的时间,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/57810112/