Разреженные проблемы собственных значений с ARPACK

Введение

ARPACK [1] является пакетом на Fortran, который предоставляет подпрограммы для быстрого поиска нескольких собственных значений/собственных векторов больших разреженных матриц. Для нахождения этих решений требуется только умножение слева на рассматриваемую матрицу. Эта операция выполняется через обратная связь интерфейс. Результат этой структуры в том, что ARPACK может находить собственные значения и векторы любой линейной функции, отображающей вектор в вектор.

Весь функционал, предоставляемый в ARPACK, содержится в двух высокоуровневых интерфейсах scipy.sparse.linalg.eigs и scipy.sparse.linalg.eigsh. eigs предоставляет интерфейсы для нахождения собственных значений/векторов вещественных или комплексных несимметричных квадратных матриц, в то время как eigsh предоставляет интерфейсы для вещественно-симметричных или комплексно-эрмитовых матриц.

Базовая функциональность

ARPACK может решать либо стандартные задачи на собственные значения вида

\[A \mathbf{x} = \lambda \mathbf{x}\]

или общие проблемы собственных значений вида

\[A \mathbf{x} = \lambda M \mathbf{x}.\]

Сила ARPACK заключается в том, что он может вычислять только указанное подмножество пар собственных значений/векторов. Это достигается с помощью ключевого слова which. Следующие значения which доступны:

• which = 'LM' : Собственные значения с наибольшей величиной (eigs, eigsh), то есть наибольшие собственные значения в евклидовой норме комплексных чисел.

• which = 'SM' : Собственные значения с наименьшей величиной (eigs, eigsh), то есть наименьшие собственные значения в евклидовой норме комплексных чисел.

• which = 'LR' : Собственные значения с наибольшей вещественной частью (eigs).

• which = 'SR' : Собственные значения с наименьшей вещественной частью (eigs).

• which = 'LI' : Собственные значения с наибольшей мнимой частью (eigs).

• which = 'SI' : Собственные значения с наименьшей мнимой частью (eigs).

• which = 'LA' : Собственные значения с наибольшим алгебраическим значением (eigsh), то есть наибольшие собственные значения, включая любой отрицательный знак.

• which = 'SA' : Собственные значения с наименьшим алгебраическим значением (eigsh), то есть наименьшие собственные значения, включая любой отрицательный знак.

• which = 'BE' : Собственные значения с обоих концов спектра (eigsh).

Заметим, что ARPACK обычно лучше находит экстремальные собственные значения, то есть собственные значения с большими величинами. В частности, используя which = 'SM' может привести к медленному времени выполнения и/или аномальным результатам. Лучший подход — использовать режим shift-invert.

Режим сдвига-инверсии

Режим сдвига-инверсии основан на следующем наблюдении. Для обобщённой задачи на собственные значения

\[A \mathbf{x} = \lambda M \mathbf{x},\]

можно показать, что

\[(A - \sigma M)^{-1} M \mathbf{x} = \nu \mathbf{x},\]

где

\[\nu = \frac{1}{\lambda - \sigma}.\]

Примеры

Представьте, что вы хотите найти наименьшие и наибольшие собственные значения и соответствующие собственные векторы для большой матрицы. ARPACK может обрабатывать множество форм ввода: плотные матрицы, такие как numpy.ndarray экземпляры, разреженные матрицы, такие как scipy.sparse.csr_matrix, или общий линейный оператор, полученный из scipy.sparse.linalg.LinearOperator. Для этого примера, для простоты, мы построим симметричную, положительно-определенную матрицу.

>>> import numpy as np
>>> from scipy.linalg import eig, eigh
>>> from scipy.sparse.linalg import eigs, eigsh
>>> np.set_printoptions(suppress=True)
>>> rng = np.random.default_rng()
>>>
>>> X = rng.random((100, 100)) - 0.5
>>> X = np.dot(X, X.T)  # create a symmetric matrix

Теперь у нас есть симметричная матрица X, с помощью которого тестируются процедуры. Сначала вычислите стандартное разложение по собственным значениям, используя eigh:

>>> evals_all, evecs_all = eigh(X)

По мере увеличения размерности X растёт, эта процедура становится очень медленной. Особенно, если требуется только несколько собственных векторов и собственных значений, ARPACK может быть лучшим вариантом. Сначала вычислим наибольшие собственные значения (which = 'LM') из X и сравнить их с известными результатами:

>>> evals_large, evecs_large = eigsh(X, 3, which='LM')
>>> print(evals_all[-3:])
[29.22435321 30.05590784 30.58591252]
>>> print(evals_large)
[29.22435321 30.05590784 30.58591252]
>>> print(np.dot(evecs_large.T, evecs_all[:,-3:]))
array([[-1.  0.  0.],       # may vary (signs)
       [ 0.  1.  0.],
       [-0.  0. -1.]])

Результаты соответствуют ожиданиям. ARPACK восстанавливает нужные собственные значения, и они совпадают с ранее известными результатами. Кроме того, собственные векторы ортогональны, как и ожидалось. Теперь попробуем найти собственные значения с наименьшей величиной:

>>> evals_small, evecs_small = eigsh(X, 3, which='SM')
Traceback (most recent call last):       # may vary (convergence)
...
scipy.sparse.linalg._eigen.arpack.arpack.ArpackNoConvergence:
ARPACK error -1: No convergence (1001 iterations, 0/3 eigenvectors converged)

Упс. Мы видим, что, как упоминалось выше, ARPACK не так хорошо справляется с поиском малых собственных значений. Есть несколько способов решения этой проблемы. Мы могли бы увеличить допуск (tol) для ускорения сходимости:

>>> evals_small, evecs_small = eigsh(X, 3, which='SM', tol=1E-2)
>>> evals_all[:3]
array([0.00053181, 0.00298319, 0.01387821])
>>> evals_small
array([0.00053181, 0.00298319, 0.01387821])
>>> np.dot(evecs_small.T, evecs_all[:,:3])
array([[ 0.99999999  0.00000024 -0.00000049],    # may vary (signs)
       [-0.00000023  0.99999999  0.00000056],
       [ 0.00000031 -0.00000037  0.99999852]])

Это работает, но мы теряем точность в результатах. Другой вариант — увеличить максимальное количество итераций (maxiter) с 1000 до 5000:

>>> evals_small, evecs_small = eigsh(X, 3, which='SM', maxiter=5000)
>>> evals_all[:3]
array([0.00053181, 0.00298319, 0.01387821])
>>> evals_small
array([0.00053181, 0.00298319, 0.01387821])
>>> np.dot(evecs_small.T, evecs_all[:,:3])
array([[ 1.  0.  0.],           # may vary (signs)
       [-0.  1.  0.],
       [ 0.  0. -1.]])

Мы получаем ожидаемые результаты, но время вычислений значительно больше. К счастью, ARPACK содержит режим, который позволяет быстро определить невнешние собственные значения: режим shift-invert. Как упомянуто выше, этот режим включает преобразование задачи на собственные значения в эквивалентную задачу с другими собственными значениями. В этом случае мы надеемся найти собственные значения около нуля, поэтому выберем sigma = 0. Преобразованные собственные значения будут удовлетворять \(\nu = 1/(\lambda - \sigma) = 1/\lambda\), поэтому наши малые собственные значения \(\lambda\) становятся большими собственными значениями \(\nu\).

>>> evals_small, evecs_small = eigsh(X, 3, sigma=0, which='LM')
>>> evals_all[:3]
array([0.00053181, 0.00298319, 0.01387821])
>>> evals_small
array([0.00053181, 0.00298319, 0.01387821])
>>> np.dot(evecs_small.T, evecs_all[:,:3])
array([[ 1.  0.  0.],    # may vary (signs)
       [ 0. -1. -0.],
       [-0. -0.  1.]])

Мы получаем ожидаемые результаты с гораздо меньшими вычислительными затратами. Обратите внимание, что преобразование из \(\nu \to \lambda\) происходит полностью в фоновом режиме. Пользователю не нужно беспокоиться о деталях.

Режим shift-invert предоставляет больше, чем просто быстрый способ получения нескольких малых собственных значений. Допустим, вы хотите найти внутренние собственные значения и собственные векторы, например, ближайшие к \(\lambda = 1\). Просто установите sigma = 1 и ARPACK позаботится об остальном:

>>> evals_mid, evecs_mid = eigsh(X, 3, sigma=1, which='LM')
>>> i_sort = np.argsort(abs(1. / (1 - evals_all)))[-3:]
>>> evals_all[i_sort]
array([0.94164107, 1.05464515, 0.99090277])
>>> evals_mid
array([0.94164107, 0.99090277, 1.05464515])
>>> print(np.dot(evecs_mid.T, evecs_all[:,i_sort]))
array([[-0.  1.  0.],     # may vary (signs)
       [-0. -0.  1.],
       [ 1.  0.  0.]]

Собственные значения выводятся в другом порядке, но все присутствуют. Обратите внимание, что режим shift-invert требует внутреннего решения матричного обращения. Это автоматически обрабатывается eigsh и eigs, но операция также может быть указана пользователем. См. строку документации scipy.sparse.linalg.eigsh и scipy.sparse.linalg.eigs подробности.

Использование LinearOperator

Теперь рассмотрим случай, когда вы хотите избежать создания плотной матрицы и использовать scipy.sparse.linalg.LinearOperator вместо этого. Наш первый линейный оператор применяет поэлементное умножение между входным вектором и вектором \(\mathbf{d}\) предоставленные пользователем оператору. Этот оператор имитирует диагональную матрицу с элементами \(\mathbf{d}\) вдоль главной диагонали, и его основное преимущество в том, что прямые и сопряжённые операции являются простыми поэлементными умножениями, отличными от умножения матрицы на вектор. Для диагональной матрицы мы ожидаем, что собственные значения будут равны элементам вдоль главной диагонали, в данном случае \(\mathbf{d}\). Собственные значения и собственные векторы, полученные с помощью eigsh сравниваются с полученными при использовании eigh при применении к плотной матрице:

>>> from scipy.sparse.linalg import LinearOperator
>>> class Diagonal(LinearOperator):
...     def __init__(self, diag, dtype='float32'):
...         self.diag = diag
...         self.shape = (len(self.diag), len(self.diag))
...         self.dtype = np.dtype(dtype)
...     def _matvec(self, x):
...         return self.diag*x
...     def _rmatvec(self, x):
...         return self.diag*x

>>> N = 100
>>> rng = np.random.default_rng()
>>> d = rng.normal(0, 1, N).astype(np.float64)
>>> D = np.diag(d)
>>> Dop = Diagonal(d, dtype=np.float64)

>>> evals_all, evecs_all = eigh(D)
>>> evals_large, evecs_large = eigsh(Dop, 3, which='LA', maxiter=1e3)
>>> evals_all[-3:]
array([1.53092498, 1.77243671, 2.00582508])
>>> evals_large
array([1.53092498, 1.77243671, 2.00582508])
>>> print(np.dot(evecs_large.T, evecs_all[:,-3:]))
array([[-1.  0.  0.],     # may vary (signs)
       [-0. -1.  0.],
       [ 0.  0. -1.]]

В этом случае мы создали быстрый и простой Diagonal оператор. Внешняя библиотека PyLops предоставляет схожие возможности в Диагональ оператор, а также несколько других операторов.

Наконец, мы рассматриваем линейный оператор, имитирующий применение шаблона первой производной. В этом случае оператор эквивалентен вещественной несимметричной матрице. Снова мы сравниваем оценённые собственные значения и собственные векторы с теми, что получены из плотной матрицы, применяющей ту же первую производную к входному сигналу:

>>> class FirstDerivative(LinearOperator):
...     def __init__(self, N, dtype='float32'):
...         self.N = N
...         self.shape = (self.N, self.N)
...         self.dtype = np.dtype(dtype)
...     def _matvec(self, x):
...         y = np.zeros(self.N, self.dtype)
...         y[1:-1] = (0.5*x[2:]-0.5*x[0:-2])
...         return y
...     def _rmatvec(self, x):
...         y = np.zeros(self.N, self.dtype)
...         y[0:-2] = y[0:-2] - (0.5*x[1:-1])
...         y[2:] = y[2:] + (0.5*x[1:-1])
...         return y

>>> N = 21
>>> D = np.diag(0.5*np.ones(N-1), k=1) - np.diag(0.5*np.ones(N-1), k=-1)
>>> D[0] = D[-1] = 0 # take away edge effects
>>> Dop = FirstDerivative(N, dtype=np.float64)

>>> evals_all, evecs_all = eig(D)
>>> evals_large, evecs_large = eigs(Dop, 4, which='LI')
>>> evals_all_imag = evals_all.imag
>>> isort_imag = np.argsort(np.abs(evals_all_imag))
>>> evals_all_imag = evals_all_imag[isort_imag]
>>> evals_large_imag = evals_large.imag
>>> isort_imag = np.argsort(np.abs(evals_large_imag))
>>> evals_large_imag = evals_large_imag[isort_imag]
>>> evals_all_imag[-4:]
array([-0.95105652, 0.95105652, -0.98768834, 0.98768834])
>>> evals_large_imag
array([0.95105652, -0.95105652, 0.98768834, -0.98768834]) # may vary

Обратите внимание, что собственные значения этого оператора все мнимые. Более того, ключевое слово which='LI' of scipy.sparse.linalg.eigs выдаёт собственные значения с наибольшей абсолютной мнимой частью (как положительной, так и отрицательной). Опять же, более продвинутая реализация оператора первой производной доступна в PyLops библиотека под названием FirstDerivative оператор.

Ссылки

[1]

opencollab/arpack-ng