Есть ли возможность инвертировать симметричную полосу (7 диагоналей) в линейное время?

Question

Мне нужно инвертировать p x p симметричную ленточную матрицу H, которая имеет 7 диагоналей. p может быть очень высоким (= 1000 или 10000).

H^{- 1} можно рассматривать как полосу, и, следовательно, мне не нужно вычислять полную обратную матрицу, а скорее ее приближение. (Например, можно предположить, что у вас 11 или 13 диагоналей.) Я ищу метод, который не предполагает распараллеливания.

Есть ли возможность построить такой алгоритм с R в линейном времени?

Благодарим за помощь.

Answer 1

Для этого, насколько мне известно, нет алгоритма линейного времени. Но вы не совсем без надежды:

1. Ваша матрица не так уж велика, поэтому использование относительно оптимизированной реализации может быть достаточно быстро для p < 10K. Например, плотный Разложение LU требует не более O(p^3), с p = 1000, что, вероятно, займет меньше секунды. На практике реализация для разреженных матриц должна обеспечивать гораздо лучшую производительность, используя преимущества разреженности;
2. Вы действительно, действительно, действительно нужно вычислить инверсию? Очень часто явное вычисление обратного может быть заменено решением эквивалентной линейной системы; с некоторыми методами, такими как итеративные решатели (например, conjugate gradient), решение линейной системы значительно более эффективно, поскольку сохраняется модель разреженности исходной матрицы, что приводит к уменьшению работы; при вычислении инверсного значения, даже если вы знаете, что это нормально для аппроксимации с помощью ленточной матрицы, все равно будет значительная сумма заполнения (добавлены ненулевые значения)

Собираем все это вместе Предлагаю вам попробовать из R matrix package для вашей матрицы. Попробуйте все доступные подписи и убедитесь, что установлена ​​высокопроизводительная реализация BLAS. Кроме того, попробуйте переписать вызов вычислить обратное:

# e.g. rewrite... 
A_inverse = solve(A) 
x = y * A_inverse 
# ... as 
x = solve(A, y) 

Это может быть более тонким для ваших целей, но есть очень большой шанс, что вы должны быть в состоянии сделать это, как это было предложено в пакете документации:

solve(a, b, ...) ## *the* two-argument version, almost always preferred to 
solve(a)   ## the *rarely* needed one-argument version 

Если все остальное терпит неудачу, то вы можете попробовать более эффективные реализации имеющихся в: Matlab, Suite Sparse, PetSC, Eigen или Intel MKL.