Hol*_*onk 4 python numpy linear-algebra scipy numpy-einsum
让A是一个(N,M,M)矩阵（N非常大），我想计算scipy.linalg.expm(A[n,:,:])每个n in range(N)。我当然可以只使用for循环，但我想知道是否有一些技巧可以以更好的方式做到这一点（例如np.einsum）。
我对矩阵求逆等其他操作也有同样的问题（在评论中解决了求逆问题）。
根据矩阵的大小和结构，您可以做得比循环更好。
\n\n假设您的矩阵可以对角化为A = V D V^(-1)（其中D对角线中有特征值并V包含相应的特征向量作为列），您可以将矩阵指数计算为
exp(A) = V exp(D) V^(-1)\n其中exp(D)仅包含其对角线中的exp(lambda)每个特征值。lambda如果我们使用指数函数的幂级数定义，这确实很容易证明。如果矩阵A进一步正规，则矩阵V是酉矩阵，因此可以通过简单地取其伴随矩阵来计算其逆矩阵。
好消息是numpy.linalg.eig，numpy.linalg.inv两者都可以很好地处理堆叠矩阵：
import numpy as np\nimport scipy.linalg\n\nA = np.random.rand(1000,10,10)\n\ndef loopy_expm(A):\n expmA = np.zeros_like(A)\n for n in range(A.shape[0]):\n expmA[n,...] = scipy.linalg.expm(A[n,...])\n return expmA\n\ndef eigy_expm(A):\n vals,vects = np.linalg.eig(A)\n return np.einsum('...ik, ...k, ...kj -> ...ij',\n vects,np.exp(vals),np.linalg.inv(vects))\n请注意，在调用 时指定操作顺序可能还有一些优化空间einsum，但我没有对此进行调查。
针对随机数组测试上述内容：
\n\nIn [59]: np.allclose(loopy_expm(A),eigy_expm(A))\nOut[59]: True\n\nIn [60]: %timeit loopy_expm(A)\n824 ms \xc2\xb1 55.7 ms per loop (mean \xc2\xb1 std. dev. of 7 runs, 1 loop each)\n\nIn [61]: %timeit eigy_expm(A)\n138 ms \xc2\xb1 992 \xc2\xb5s per loop (mean \xc2\xb1 std. dev. of 7 runs, 10 loops each)\n那已经很好了。如果你足够幸运，你的矩阵都是正常的（比如说，因为它们是真正对称的）：
\n\nA = np.random.rand(1000,10,10)\nA = (A + A.transpose(0,2,1))/2\n\ndef eigy_expm_normal(A):\n vals,vects = np.linalg.eig(A)\n return np.einsum('...ik, ...k, ...jk -> ...ij',\n vects,np.exp(vals),vects.conj())\n请注意输入矩阵的对称定义和 模式内的转置einsum。结果：
In [80]: np.allclose(loopy_expm(A),eigy_expm_normal(A))\nOut[80]: True\n\nIn [79]: %timeit loopy_expm(A)\n878 ms \xc2\xb1 89.7 ms per loop (mean \xc2\xb1 std. dev. of 7 runs, 1 loop each)\n\nIn [80]: %timeit eigy_expm_normal(A)\n55.8 ms \xc2\xb1 868 \xc2\xb5s per loop (mean \xc2\xb1 std. dev. of 7 runs, 10 loops each)\n对于上述示例形状来说，速度提高了 15 倍。
\n\n应该注意的是，scipy.linalg.eigm根据文档，它使用 Pad\xc3\xa9 近似值。这可能意味着，如果您的矩阵病态，特征值分解可能会产生与 不同的结果scipy.linalg.eigm。我不熟悉这个函数是如何工作的，但我希望它对于病理输入来说更安全。