Python `(N,M,M)` 矩阵的 `expm`

Question

让 A 成为一个 (N,M,M) 矩阵（N 非常大），我想为每个 n in range(N) 计算 scipy.linalg.expm(A[n,:,:])。我当然可以只使用 for 循环，但我想知道是否有一些技巧可以更好地做到这一点（比如 np.einsum）。

我对其他操作也有同样的问题，比如求逆矩阵（求逆在评论中解决）。

Answer 1

根据矩阵的大小和结构，您可以做得比循环更好。

假设您的矩阵可以对角化为 A = V D V^(-1)（其中 D 在其对角线上具有特征值，并且 V 包含相应的特征向量作为列），您可以计算矩阵指数作为

exp(A) = V exp(D) V^(-1)

其中 exp(D) 只包含 exp(lambda) 其对角线上的每个特征值 lambda。如果我们使用指数函数的幂级数定义，这真的很容易证明。如果矩阵 A 进一步正规，则矩阵 V 是酉矩阵，因此可以通过简单地取它的伴随来计算它的逆矩阵。

好消息是 numpy.linalg.eig 和 numpy.linalg.inv 都可以很好地处理堆叠矩阵：

import numpy as np
import scipy.linalg

A = np.random.rand(1000,10,10)

def loopy_expm(A):
    expmA = np.zeros_like(A)
    for n in range(A.shape[0]):
        expmA[n,...] = scipy.linalg.expm(A[n,...])
    return expmA

def eigy_expm(A):
    vals,vects = np.linalg.eig(A)
    return np.einsum('...ik, ...k, ...kj -> ...ij',
                     vects,np.exp(vals),np.linalg.inv(vects))

请注意，在对 einsum 的调用中指定操作顺序可能有一些优化空间，但我没有对此进行调查。

正在为随机数组测试以上内容：

In [59]: np.allclose(loopy_expm(A),eigy_expm(A))
Out[59]: True

In [60]: %timeit loopy_expm(A)
824 ms ± 55.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [61]: %timeit eigy_expm(A)
138 ms ± 992 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

已经很不错了。如果你足够幸运，你的矩阵都是正常的（比如，因为它们是实对称的）：

A = np.random.rand(1000,10,10)
A = (A + A.transpose(0,2,1))/2

def eigy_expm_normal(A):
    vals,vects = np.linalg.eig(A)
    return np.einsum('...ik, ...k, ...jk -> ...ij',
                     vects,np.exp(vals),vects.conj())

注意输入矩阵的对称定义和 einsum 模式内的转置。结果：

In [80]: np.allclose(loopy_expm(A),eigy_expm_normal(A))
Out[80]: True

In [79]: %timeit loopy_expm(A)
878 ms ± 89.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [80]: %timeit eigy_expm_normal(A)
55.8 ms ± 868 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

上述示例形状的速度提高了 15 倍。

需要注意的是 scipy.linalg.eigm 根据文档使用 Padé 近似。这可能意味着如果您的矩阵是 ill-conditioned，则特征值分解可能会产生与 scipy.linalg.eigm 不同的结果。我不熟悉这个功能是如何工作的，但我希望它对病态输入更安全。

Python `(N,M,M)` 矩阵的 `expm`

Python `expm` of an `(N,M,M)` matrix

python

numpy

linear-algebra

scipy

numpy-einsum