Идентичность Вудбери для быстрой инверсии матрицы - медленнее, чем ожидалось

Question

Тождество матрицы Вудбери гласит, что обратное значение поправки на ранг-k некоторой матрицы можно вычислить, выполнив поправку на ранг-k к инверсии исходной матрицы.

Если A является p × p матрицей полного ранга, ранг которой скорректирован на UCV, где 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 < 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, моя функция Вудбери не так быстро, как я ожидал,Я ожидаю, что время инвертирования будет расти кубически с размерностью p.Но мои результаты:

Мой вопрос: почему метод Вудбери такой медленный?Это просто потому, что я сравниваю код Python с LAPACK или что-то еще происходит?

РЕДАКТИРОВАТЬ: Мои эксперименты с einsum () против трансляции

Iреализованы три версии: (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().

Answer 1

Как вы упоминаете, инвертирование A + UCV может быть выполнено быстрее с использованием техники Вудбери в случае, когда A имеет диагональ .То есть в формуле Вудбери ваши умножения на A^{-1} должны происходить за O(p x m) время вместо O(p x m x p), поскольку все, что вы делаете, это масштабирование строк / столбцов правого / левого члена.

Однако, это не то, что вы делаете в следующем коде!

def woodbury(A, U, V, k):
    A_inv = np.diag(1./np.diag(A))
    B_inv = np.linalg.inv(np.eye(k) + V @ A_inv @ U)
    return A_inv - (A_inv @ U @ B_inv @ V @ A_inv)

Ваш A_inv является полной p x p матрицей!Да, диагональ - это единственная часть, которая содержит ненулевые элементы, но арифметика со всеми нулевыми элементами все равно будет выполняться в этой плотной матрице!Вместо этого вы должны использовать возможности вещания Numpys, чтобы избежать этой ненужной работы.(Или, в качестве разреженной диагональной матрицы, используя модуль Scipy sparse .)

Например,

def efficient_woodbury(A, U, V, k):
    A_inv_diag = 1./np.diag(A)  # note! A_inv_diag is a vector!
    B_inv = np.linalg.inv(np.eye(k) + (V * A_inv_diag) @ U)
    return np.diag(A_inv_diag) - (A_inv_diag.reshape(-1,1) * U @ B_inv @ V * A_inv_diag)

Продукт V * A_inv_diag эквивалентен вашему V @ A_inv но работает в O(p x k) времени, а не O(p x k x p).Аналогично для других замененных продуктов.

Я перезапустил время на моей (немного более быстрой) машине и создал следующий график:

Намного более четкое различие в производительности!