我有类似的示例代码
import numpy as np
nbins = 1301
init = np.ones(nbins*2+1)*(-1)
init[0] = 0
init[1::2] = 0
z0 = np.array(init, dtype=np.complex128, ndmin=1)
initn = z0.view(np.float64)
deltasf = np.linspace(-10, 10, nbins)
gsf = np.linspace(1, 2, nbins)
def jacobian_mbes_test(Y, t):
leny = len(Y)
J = np.zeros((leny, leny))
J[0,0] = -10
J[0,1] = 0
J[0, 2::4] = gsf
J[1,0] = 0
J[1,1] = -10
J[1,3::4] = gsf
J[2::4, 0] = gsf*Y[4::4]
J[3::4, 0] = gsf*Y[5::4]
J[4::4, 0] = -4*gsf*Y[2::4]
J[2::4, 1] = -gsf*Y[5::4]
J[3::4, 1] = gsf*Y[4::4]
J[4::4, 1] = -4*gsf*Y[3::4]
J[(range(2, leny, 4), range(2, leny, 4))] = -0.8
J[(range(3, leny, 4), range(3, leny, 4))] = -0.8
J[(range(4, leny, 4), range(4, leny, 4))] = -0.001
J[(range(5, leny, 4), range(5, leny, 4))] = -0.001
J[(range(2, leny, 4), range(3, leny, 4))] = deltas
J[(range(3, leny, 4), range(2, leny, 4))] = -deltas
J[(range(2, leny, 4), range(4, leny, 4))] = gsf*Y[0]
J[(range(2, leny, 4), range(5, leny, 4))] = -gsf*Y[1]
J[(range(3, leny, 4), range(4, leny, 4))] = gsf*Y[1]
J[(range(3, leny, 4), range(5, leny, 4))] = gsf*Y[0]
J[(range(4, leny, 4), range(2, leny, 4))] = -4*gsf*Y[0]
J[(range(4, leny, 4), range(3, leny, 4))] = -4*gsf*Y[1]
return J
它使用scipy.integrate.deint 为一组微分方程计算雅可比
不幸的是，函数jacobian_mbes_test的性能不是很好。%timeit jacobian_mbes_test(initn, 0)给出了12.2毫秒，尽管我试图通过切片(？(高效地完成所有工作。
将函数更改为(因为我想更复杂的索引分类需要花费大部分时间(
def jacobian_mbes_test2(Y, t):
leny = len(Y)
J = np.zeros((leny, leny))
J[0,0] = -10
J[0,1] = 0
J[0, 2::4] = gsf
J[1,0] = 0
J[1,1] = -10
J[1,3::4] = gsf
J[2::4, 0] = gsf*Y[4::4]
J[3::4, 0] = gsf*Y[5::4]
J[4::4, 0] = -4*gsf*Y[2::4]
J[2::4, 1] = -gsf*Y[5::4]
J[3::4, 1] = gsf*Y[4::4]
J[4::4, 1] = -4*gsf*Y[3::4]
return J
将这个时间减少到4.6ms，这仍然不是很大。有趣的是，如果我在这样的函数之外计时这些任务
J = np.zeros((len(initn), len(initn))
%timeit J[4::4, 1] = -4*gsf*initn[3::4]
我最多得到10µs，所以我希望第二个函数最多(！(100µs。在我不知道的这个函数中，内存中是否有一些复制？
有没有更有效的方法来为某些索引赋值？
问题来自矩阵大小(206 MiB(。事实上，矩阵相当大，每次调用函数时(实际上(都会填充零。假设所有值将在12.2毫秒内物理写入内存，则吞吐量将为5206**2*8/12.2e-3/1024**3 = 16.5 GiB/s，这对于顺序工作负载来说非常好(尤其是在普通个人计算机上，其顺序通常几乎不超过25GiB/s(。在实践中，内存不会设置为零，并且大部分计算时间来自虚拟内存的管理。
np.zeros((leny, leny))在内部请求您的操作系统(OS(保留一些空间，并将其初始化为零。您的操作系统实际上会分配许多页面(通常为4 KiB的块(，并将它们标记为稍后设置为零(延迟初始化(。因此，这种工艺非常便宜。然而，当稍后执行J[0,0] = -10时，它将数据写入需要初始化的给定第一个触摸页中。这被称为页面错误。这个过程是昂贵的(尤其是在顺序过程中(，因为它会停止进程的执行，在物理内存中找到一些位置，填充整个页面，然后恢复进程，对于每个页面！在您的情况下，许多要初始化的数组单元格会导致页面被填充，这比实际填充单元格的成本高得多。
CPU缓存在这段代码中也很重要，因为当在RAM中进行计算时，页面填充会慢得多。当数据量太大而无法放入最后一级缓存(通常为1-64MIB(时，就会发生这种情况。
内存访问模式强烈影响相当大的数据结构的性能(即不适合CPU缓存(。事实上，随机访问模式或大步跨步访问确实对性能不利，因为处理器无法轻松预测要提前加载的RAM块的位置(更不用说RAM延迟相当大(。连续访问要快得多。
当你在更简单的情况下衡量时间时，所有这些都应该考虑在内。
查看虚拟内存影响的一个简单方法是使用np.full((leny, leny), 0.0)而不是np.zeros(leny, leny)。这在我的机器上慢了3倍。在这种情况下，所有的页面都需要填写。您还可以创建测量时间来多次填充J矩阵。第一次应该慢很多(在我的机器上慢2.2倍(。
一种解决方案是一次分配和初始化数组，只填充/回收一些值，但这有点棘手。在最坏的情况下，大多数矩阵需要归零，这在您的情况下太慢了，因为它会受到内存吞吐量的限制。
一个简单的解决方案是使用空间矩阵。可以使用模块scipy.sparse来创建空间矩阵。它们的计算速度通常较慢，但如果有很多零(这是你的情况(，它们会更快、更紧凑。这里有一个例子：
from scipy.sparse.lil import lil_matrix
def jacobian_mbes_test(Y, t):
leny = len(Y)
J = lil_matrix((leny, leny))
J[0,0] = -10
J[0,1] = 0
J[0, 2::4] = gsf
J[1,0] = 0
J[1,1] = -10
J[1,3::4] = gsf
J[2::4, 0] = gsf*Y[4::4]
J[3::4, 0] = gsf*Y[5::4]
J[4::4, 0] = -4*gsf*Y[2::4]
J[2::4, 1] = -gsf*Y[5::4]
J[3::4, 1] = gsf*Y[4::4]
J[4::4, 1] = -4*gsf*Y[3::4]
J[(range(2, leny, 4), range(2, leny, 4))] = -0.8
J[(range(3, leny, 4), range(3, leny, 4))] = -0.8
J[(range(4, leny, 4), range(4, leny, 4))] = -0.001
J[(range(5, leny, 4), range(5, leny, 4))] = -0.001
J[(range(2, leny, 4), range(3, leny, 4))] = deltas
J[(range(3, leny, 4), range(2, leny, 4))] = -deltas
J[(range(2, leny, 4), range(4, leny, 4))] = gsf*Y[0]
J[(range(2, leny, 4), range(5, leny, 4))] = -gsf*Y[1]
J[(range(3, leny, 4), range(4, leny, 4))] = gsf*Y[1]
J[(range(3, leny, 4), range(5, leny, 4))] = gsf*Y[0]
J[(range(4, leny, 4), range(2, leny, 4))] = -4*gsf*Y[0]
J[(range(4, leny, 4), range(3, leny, 4))] = -4*gsf*Y[1]
return J
这在我的机器上快了2倍。有很多不同类型的稀疏矩阵，每种矩阵都适用于一些特定的用例。可能存在比LIL更好的稀疏矩阵类型。
LIL空间矩阵对于具有动态结构的空间矩阵是非常好的。CSR对于静态CSR非常好。因此，例如，您可以在计划写入一些数据的位置提前构建一个CSR矩阵，其中填充NaN值，然后在需要时复制它并用新值填充新的空间矩阵。这比我的机器快10-15%。使用更好的访问模式可能会进一步减少时间。