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

Question

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

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

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

Спасибо за помощь.

Answer 1

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

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

Собирая все это вместе, я бы посоветовал вам попробовать пакет матрицы R для вашей матрицы. Попробуйте все доступные подписи и убедитесь, что у вас установлена ​​высокопроизводительная реализация 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.