计算科学 - 使用 Python 加速线性变换 - 吾爱随笔录

使用 Python 加速线性变换

计算科学 Python 并行计算软件麻木的多核

2021-12-04 00:56:47

在光波前传播问题中，我需要进行过多的傅里叶类型计算：

import numpy as np

N=10000
input=np.random.random(N)

x=np.linspace(-1,1,N)
y=np.linspace(-2,2,N)
X,Y=np.meshgrid(x,y,sparse=True)

output=np.dot(np.exp(1j*2*np.pi*X*Y),input)

然而，真正的内核（这里是复指数）更复杂。特别是，它不是周期性的，因此不可能进行 FFT 类型的缩减！因此，简而言之，我需要在一个非常大的网格上评估三角函数，然后进行矩阵向量乘法。

使用纯numpy执行此操作需要一个多星期才能在中产阶级台式机（具有八个内核）上执行我的计算。使用numexpr可将运行时间减少大约 1/3。相反，使用并行 Cython可将运行时间进一步减少约 8/10 倍。

尽管如此，这仍然太慢了。我根本不是科学计算方面的专家。所以，对我来说，很难判断如何进一步加快速度。因此，我的问题是：

是否存在更适合此问题的其他软件框架（numba、weave、OpenCL...）？
什么样的硬件升级才有意义？（显然，主要的计算成本是函数评估。）

编辑：连续问题是一个广义的一维菲涅耳传播器

u (x) = \int_{- \infty}^{\infty} \exp (i π \frac{(x - y)^{2}}{λ z (x, y)}) \frac{u_{0} (y)}{\sqrt{i λ z (x, y)}} d y,

$u(x)= \int_{-\infty}^\infty \exp\left(i\pi\frac{(x-y)^2}{\lambda z(x,y)}\right) \frac{u_0(y)}{\sqrt{i\lambda z(x,y)}} dy,$ 在哪里

λ

$\lambda$ 是光的波长，z 是一个平滑的正函数，但或多或少是任意的

2个回答

假设您的内核有些平滑，请使用 low-rank approximation。

这是一个天真的例子：

import numpy as np

N=2000
input=np.random.random(N)

x=np.linspace(-1,1,N)
y=np.linspace(-2,2,N)
X,Y=np.meshgrid(x,y,sparse=True)
A = np.exp(1j*2*np.pi*X*Y)

output = np.dot(A, input)

U,S,V = np.linalg.svd(A)

# find truncation rank for given tolerance
k = np.nonzero(S < 1e-12)[0][0]
print('Truncation rank:', k)
Xk = U[:,:k].dot(np.diag(S[:k]))
Yk = V[:k,:]

output2 = Xk.dot(Yk.dot(input))

# check quality of approximation
print(np.linalg.norm(output - output2, np.inf))

%timeit np.dot(A, input)
%timeit Xk.dot(Yk.dot(input))

输出：

Truncation rank: 32
6.3535593536e-12
100 loops, best of 3: 3.9 ms per loop
10000 loops, best of 3: 61.4 µs per loop

上面所做的是计算矩阵的截断奇异值分解，并使用得到的 A 的低秩近似来加速矩阵向量乘积。所以，而不是使用 $A$ , 你使用一个近似值

A \approx X_{k} Y_{k}^{T}, X_{k}, Y_{k} \in R^{N \times k}

$A \approx X_k Y_k^T, \qquad X_k,Y_k \in \mathbb R^{N\times k}$ 和

k ≪ N

$k \ll N$ .

在示例中，这会导致超过 60 倍的加速，并且只会随着 N 变大而变得更好。那是因为完整的矩阵向量乘积具有复杂性 $\mathcal O(N^2)$ ，而低秩矩阵向量积具有复杂性 $\mathcal O(Nk)$ . 结果并不准确，但您可以通过使用截断等级来影响准确性。

当然，计算 SVD 非常慢。相反，有更快的替代方法来计算矩阵的低秩近似，例如：快速随机 SVD，或者我推荐的自适应交叉逼近 (ACA)。后者是一种黑盒算法，用于计算低秩近似，甚至不需要计算整个矩阵 A，因此非常适合非常大的问题。它也很容易实现。

关于 ACA 的原始文献（特别是 M. Bebendorf 从 2000 年开始撰写的）有些难以获得。取而代之的是，这是一篇以基本形式描述 ACA 的随机论文：（链接）。

如果您的内核相对简单，甚至可以通过分析得出一个可分离的近似值，这直接导致低秩近似值，而不必通过算法计算它。查阅有关“多极展开”和“快速多极方法”的文献，了解许多这方面的例子。

最后，原始程序中的瓶颈实际上是计算 A，而不是矩阵向量积。ACA 将消除这一点，因为它不需要完整的矩阵。

更新：

我刚刚意识到 scipy 已经包含了一种在鲜为人知的scipy.linalg.interpolative模块中逼近 SVD 的快速方法。这是一个例子：

import scipy.linalg.interpolative as sli

U,S,V = sli.svd(A, 1e-12)
print('Rank:', U.shape[1])

Xk = U.dot(np.diag(S))
Yk = V.conj().T

output3 = Xk.dot(Yk.dot(input))

# check quality of approximation
print(np.linalg.norm(output - output3, np.inf))

和输出：

Rank: 37
2.62815545927e-10

精确 SVD 和近似 SVD 之间的时间差异（同样，由于涉及的渐近，这将在更大的 N 下变得更好）：

%timeit np.linalg.svd(A)
%timeit sli.svd(A, 1e-12)

1 loop, best of 3: 11.4 s per loop
1 loop, best of 3: 297 ms per loop

但是，这仍然需要您首先计算整个矩阵 A 并且不能摆脱该瓶颈。然而，如果您只需要计算 A 一次然后有许多矩阵向量乘积，它可能是一个合适的选择。

矩阵向量乘法是受 RAM 限制的。您对读取的每个浮点数进行两次算术运算。因此，主要成本是将整个矩阵存储在慢速、非缓存的核心内存中并将其读回。

您应该考虑即时评估您的“复杂内核”，而不是将其保存到内存中：

output = np.zeros(n)
for j in range(n):
  for i in range(n):
    output[i] += kernel(i, j) * input[j]

（警告：在你对它进行 cythonize 之前，这个循环会非常慢）。

如果您必须使用每个矩阵制作单个矩阵向量乘积，就像在您的示例中一样，这可能会是一个主要的加速。

如果您必须使用相同的矩阵制造许多产品，那么这取决于：尝试两种方法并比较它们。（另外，如果您同时拥有所有可用的向量，请将它们存储在矩阵中，并将此操作转换为矩阵-矩阵乘法）。在硬件部门，首先要检查的是确保你有足够的 RAM 来存储整个矩阵——如果你必须打磁盘，那么事情就会变得很糟糕。 .

其它你可能感兴趣的问题

上一篇稀疏矩阵应该有多“稀疏”才能看到好处？下一篇用多项式逼近阶跃函数