许多矩阵指数的矢量化计算

Question

我有很多2×2矩阵，A，A.shape == (2, 2, 7324)，我必须计算所有这些矩阵的指数。不幸的是，scipy.linalg.expm一次只接受一个矩阵，所以我必须循环计算，

import numpy
import scipy.linalg
numpy.random.seed(0)
A = numpy.random.rand(2, 2, 7324)
out = numpy.array([scipy.linalg.expm(A[:, :, k]) for k in range(A.shape[-1])])
out = numpy.moveaxis(out, 0, -1)

有没有办法避免这种循环？

编辑：相应的 scipy 错误：#12838

Answer 1

从 SciPy 版本 1.9.0(发布于 7 月 28 日、22 日(开始，您可以传递scipy.linalg.expm最后两个维度为正方形的数组(即形状为(..., n, n)的数组(，它将计算矩阵指数 以矢量化的方式。

文档在这里。

Answer 2

通过浏览代码，https://github.com/scipy/scipy/blob/v1.4.1/scipy/sparse/linalg/matfuncs.py#L550-L595 似乎它实际上只假设一个矩阵。因此，如果你想避免python循环，你需要提取相关部分，然后自己对它们进行cython化或numpy-vector化。只有迎合密集矩阵才能使其更容易。

Answer 3

这是PADE近似的代码。这仅适用于 1 个数组，您必须对其进行矢量化/胞化。或者只是使用来自 numba 的@njit

def matrixexponen(n,a):
q = 6
a2 = a.copy ( )
a_norm = np.linalg.norm ( a2, np.inf )
ee = ( int ) ( np.log2 ( a_norm ) ) + 1
s = max ( 0, ee + 1 )
a2 = a2 / ( 2.0 ** s )
x = a2.copy ( )
c = 0.5
e = np.eye ( n, dtype = np.complex64 ) + c * a2
d = np.eye ( n, dtype = np.complex64 ) - c * a2
p = True
for k in range ( 2, q + 1 ):
c = c * float ( q - k + 1 ) / float ( k * ( 2 * q - k + 1 ) )
x = np.dot ( a2, x )
e = e + c * x

if ( p ):
d = d + c * x
else:
d = d - c * x
p = not p
#  E -> inverse(D) * E
e = np.linalg.solve ( d, e )
#  E -> E^(2*S)
for k in range ( 0, s ):
e = np.dot ( e, e )
return e