我正在处理python 64位上的非常大的数据集，需要一些帮助来优化我的插值代码。
我习惯于 使用 numpy 来避免循环，但这里有 2 个循环我找不到避免的方法。
主要问题还在于，当我使用 numpy 时，我需要计算的数组大小会出现内存错误，所以我切换到 scipy 稀疏数组，它可以工作，但计算 2 个左循环需要太多时间......
我尝试使用 numpy.fromfunction 迭代构建我的矩阵，但它无法运行，因为数组的大小太大。
我已经阅读了很多关于构建大型数组的文章，但是与我必须构建的数组相比，被问到的数组太简单了，因此解决方案在这里不起作用。
我无法减小数据集的大小，因为它是我已经在 10x10 图块中平铺的点云。
这是我的插值代码：
z_int = ss.dok_matrix((x_int.shape))
n,p = x_obs.shape
m = y_obs.shape[0]
a1 = ss.coo_matrix( (n, 3), dtype=np.int64 )
a2 = ss.coo_matrix( (3, 3), dtype=np.int64 )
a3 = ss.dok_matrix( (n, m))
a4 = ss.coo_matrix( (3, n), dtype=np.int64)
b = ss.vstack((z_obs, ss.coo_matrix( (3, 1), dtype=np.int64 ))).tocoo()
a1 = ss.hstack((ss.coo_matrix(np.ones((n,p))), ss.coo_matrix(x_obs), ss.coo_matrix(y_obs)))
shape_a3 = a3.shape[0]
for l in np.arange(0, shape_a3):
for c in np.arange(0, shape_a3) :
if l == c:
a3[l, c] = rho
else:
a3[l, c] = phi(x_obs[l] - x_obs[c], y_obs[l] - y_obs[c])
a4 = a1.transpose()
a12 = ss.vstack((a1, a2))
a34 = ss.vstack((a3, a4))
a = ss.hstack((a12, a34)).tocoo()
x = spsolve(a, b)
for i in np.arange(0, z_int.shape[0]):
for j in np.arange(0, z_int.shape[0]):
z_int[i, j] = x[0] + x[1] * x_int[i, j] + x[2] * y_int[i, j] + np.sum(x[3:] * phi(x_int[i, j] - x_obs, y_int[i, j] - y_obs).T)
return z_int.todense()
其中 dist() 是一个计算距离的函数，phi 是以下：
return dist(dx, dy) ** 2 * np.log(dist(dx, dy))
我需要代码运行得更快，我知道它可能写得很糟糕，但我想学习如何编写更优化的东西来提高我的编码技能。
该代码很难遵循，而且速度很慢，这是可以理解的。 稀疏矩阵的迭代甚至比密集数组的迭代慢。 我几乎希望您从使用密集数组的小型工作示例开始，然后再担心使其适用于大情况。 我不会尝试全面的修复或加速，只是在这里和那里啃咬。
这第一个a1创造对你没有任何作用(除了浪费时间)。 Python 不是一种编译语言，您可以在开始时定义变量的类型。 第二个赋值之后a1是一个稀疏矩阵，因为这是hstack创建的，而不是因为前一个coo赋值。
a1 = ss.coo_matrix( (n, 3), dtype=np.int64 )
...
a1 = ss.hstack((ss.coo_matrix(np.ones((n,p))), ss.coo_matrix(x_obs), ss.coo_matrix(y_obs)))
初始化dok矩阵、zint和a3是正确的，因为您需要迭代以填充值。 但我喜欢看到这种初始化更接近循环，而不是回到顶部。 我会使用lil而不是dok，但我不确定这是否更快。
for l in np.arange(0, shape_a3):
for c in np.arange(0, shape_a3) :
if l == c:
a3[l, c] = rho
else:
a3[l, c] = phi(x_obs[l] - x_obs[c], y_obs[l] - y_obs[c])
l==c测试确定主对角线。 有一些方法可以制作对角矩阵。 但看起来您正在设置a3的所有元素。 如果是这样，为什么要使用较慢的稀疏方法？
什么是phi它需要标量输入吗？x_obs[:,None]-x_obs应该直接给出一个 (n，n) 数组。
spsolve产生什么？x，稀疏或密集。 从您在z_int循环中的使用来看，它看起来像一个 1D 密集数组。 看起来您正在设置z_int的所有值。
如果phi需要一个(n，n)数组，我认为
x[0] + x[1] * x_int[i, j] + x[2] * y_int[i, j] + np.sum(x[3:] * phi(x_int[i, j] - x_obs, y_int[i, j] - y_obs).T)
x[0] + x[1] * x_int + x[2] * y_int + np.sum(x[3:]) * phi(x_int-x_obs, y_int-y_obs).T)