python numpy 求逆_numpy-用于快速矩阵求逆的伍德伯里身份-比...

Woodbury matrix identity指出,可以通过对原始矩阵的逆进行秩k校正来计算某些矩阵的逆k校正.

如果A是通过UCV进行秩校正的p×p满秩矩阵,其中U为p×k,C为k×k,V为k×p,则伍德伯里恒等式为：

(A + UCV)^{-1} = A^{-1} - A^{-1} U (C^{-1} + V A^{-1} U)^{-1} V A^{-1}

关键在于,您不需反转p×p矩阵,而需反转k×k矩阵.在许多应用中,我们可以假设k ＜k. p.在某些情况下,例如,如果A是对角矩阵,则A的反相可能会很快.

我在这里实现了这一点,假设A是对角线,而C是标识：

def woodbury(A, U, V, k):

A_inv = np.diag(1./np.diag(A)) # Fast matrix inversion of a diagonal.

B_inv = np.linalg.inv(np.eye(k) + V @ A_inv @ U)

return A_inv - (A_inv @ U @ B_inv @ V @ A_inv)

理智通过验证

n = 100000

p = 1000

k = 100

A = np.diag(np.random.randn(p))

U = np.random.randn(p, k)

V = U.T

M = U @ V + A

M_inv = woodbury(A, U, V, k)

assert np.allclose(M @ M_inv, np.eye(p))

但是,当我将其与numpy.linalg.inv进行实际比较时,我的Woodbury函数并没有我期望的那么快.我希望时间反转以三次方p增长.但是我的结果是：

我的问题是：为什么伍德伯里方法这么慢？仅仅是因为我正在将Python代码与LAPACK进行比较,还是发生了其他事情？

编辑：我对einsum()与广播的实验

我实现了三个版本：(1)使用einsum和einsum_path,(2)使用每个接受的答案的广播,以及(3)使用两者.这是我使用einsum的实现,使用einsum_path进行了优化：

def woodbury_einsum(A, U, V, k):

A_inv = np.diag(1./np.diag(A))

tmp = np.einsum('ab,bc,cd->ad',

V, A_inv, U,

optimize=['einsum_path', (1, 2), (0, 1)])

B_inv = np.linalg.inv(np.eye(k) + tmp)

tmp = np.einsum('ab,bc,cd,de,ef->af',

A_inv, U, B_inv, V, A_inv,

optimize=['einsum_path', (0, 1), (0, 1), (0, 2), (0, 1)])

return A_inv - tmp

结果在这里：

因此,避免使用对角矩阵的矩阵乘法的计算成本比使用einsum()优化矩阵乘法的顺序和内存占用要快.