如何使用 numpy 或 scipy 优化非常大的数组构造

Question

我在 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))

我需要代码运行得更快，我知道它可能写得很糟糕，但我想学习如何编写更优化的代码来提高我的编码技能。

Answer 1

该代码难以理解，而且速度慢是可以理解的。稀疏矩阵上的迭代甚至比密集数组上的迭代还要慢。我几乎希望您从一个使用密集数组的小型工作示例开始，然后再担心使其适用于大型案例。我不打算尝试全面修复或加快速度，只是在这里和那里吃点东西。

第一个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 循环中的使用来看，它看起来像一个一维密集数组。看起来您正在设置 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)