scipy - 科学中的逆矩阵积的痕迹

Question

我希望有效地计算逆矩阵乘积的轨迹，即Tr(A ^ -1 B)。我在许多其他地方都使用A的逆，因此我可以访问Cholesky分解或显式逆。

幼稚的事情将是以下两项之一:

np.trace(np.dot(invA, B))
np.trace(linalg.cho_solve(choA, B))

但是，这似乎是浪费的，因为它仅在使用对角线计算轨迹之前就在O(N ^ 3)处计算完整矩阵A ^ -1B。

给定显式逆，则O(N ^ 2)解决方案将是:

np.sum(invA.T * B)

尽管这需要显式逆，但这是不希望的。

我认为这样做的理想方法是仅根据给定的Cholesky分解计算A ^ -1 B的对角元素，然后简单地求和。可以使用scipy吗？还是考虑到Cholesky分解，还有另一种以数值稳定的方式计算Tr(invA * B)的方法吗？

Answer 1

@DSM直观建议的einsum似乎就是fastest way，因为它只会计算所需的条件。在python中手动执行此操作会比较慢，并且除非numpy / scipy中有专门的例程(一种数学技巧)，或者没有用较低级别的语言编写某些东西，否则我找不到更好的方法。为了检查时间，我们设置了一个虚拟数组，

import numpy as np
from scipy import linalg

invA = np.random.rand(100,100)
B = np.random.rand(100,100)
choA = linalg.cholesky(a)

尝试幼稚的建议，

%timeit np.sum(invA.T * B)
10000 loops, best of 3: 38.5 us per loop

使用逆的方法

%timeit np.sum(invA.T * B)
10000 loops, best of 3: 39.1 us per loop

似乎差不多。最后使用 einsum，

%timeit np.einsum("ij,ji", invA, B)
100000 loops, best of 3: 17.4 us per loop

大约快一倍。