scipy.sparse.linalg.

lobpcg#

scipy.sparse.linalg.lobpcg(A, X, B=None, M=None, Y=None, tol=None, maxiter=None, наибольшее=True, уровень детализации вывода=0, retLambdaHistory=False, retResidualNormsHistory=False, restartControl=20)[источник]#

Метод локально оптимального блочного предобусловленного сопряженного градиента (LOBPCG).

LOBPCG — это предобусловленный решатель собственных значений для больших вещественных симметричных и комплексных эрмитовых определённых обобщённых задач на собственные значения.

Параметры:

A{разреженная матрица, ndarray, LinearOperator, вызываемый объект}: Эрмитов линейный оператор задачи, обычно задаваемый разреженной матрицей. Часто называется "матрицей жесткости".
Xndarray, float32 или float64: Начальное приближение для k Собственные векторы (не разреженные). Если A имеет shape=(n,n) затем X должен иметь shape=(n,k).
B{разреженная матрица, ndarray, LinearOperator, вызываемый объект}: Опционально. По умолчанию B = None, что эквивалентно единичной матрице. Правая часть оператора в обобщённой задаче на собственные значения, если присутствует. Часто называется «матрицей масс». Должна быть эрмитовой положительно определённой.
M{разреженная матрица, ndarray, LinearOperator, вызываемый объект}: Опционально. По умолчанию M = None, что эквивалентно тождеству. Предобуславливатель, направленный на ускорение сходимости.
Yndarray, float32 или float64, по умолчанию: None: An n-by-sizeY ndarray ограничений с sizeY < n. Итерации будут выполняться в B-ортогональное дополнение пространства столбцов Y. Y должен быть полного ранга, если присутствует.
tolскаляр, опционально: Значение по умолчанию tol=n*sqrt(eps). Допуск решателя для критерия остановки.
maxiterint, по умолчанию: 20: Максимальное количество итераций.
наибольшееbool, по умолчанию: True: Если True, решать для наибольших собственных значений, иначе для наименьших.
уровень детализации выводаint, необязательный: По умолчанию verbosityLevel=0 нет вывода. Управляет стандартным/экранным выводом решателя.
retLambdaHistorybool, по умолчанию: False: Возвращать ли историю итеративных собственных значений.
retResidualNormsHistorybool, по умолчанию: False: Возвращать ли историю итераций норм невязок.
restartControlint, опционально.: Итерации перезапускаются, если остатки скачут 2**restartControl раз по сравнению с наименьшим записанным в retResidualNormsHistory. По умолчанию restartControl=20, делая перезапуски редкими для обратной совместимости.

Возвращает:

lambdandarray формы (k, ).: Массив k приближённые собственные значения.
vndarray той же формы, что и X.shape.: Массив k приближённые собственные векторы.
lambdaHistoryndarray, опционально.: История собственных значений, если retLambdaHistory является True.
уже не использовались; теперь они помечены как устаревшие.ndarray, опционально.: История норм невязок, если retResidualNormsHistory является True.

Примечания

Итерационный цикл выполняется maxit=maxiter (20 если maxit=None) итераций максимум и завершается раньше, если достигнута допустимая погрешность. Нарушая обратную совместимость с предыдущей версией, LOBPCG теперь возвращает блок итерационных векторов с наилучшей точностью, а не последний проитерированный, как средство от возможной расходимости.

Если X.dtype == np.float32 и пользовательские операции/умножения на A, B, и M все сохраняют np.float32 тип данных, все вычисления и вывод осуществляются в np.float32.

Размер вывода истории итераций равен количеству лучших (ограничено maxit) итераций плюс 3: начальная, конечная и постобработка.

Если оба retLambdaHistory и retResidualNormsHistory являются True, возвращаемый кортеж имеет следующий формат (lambda, V, lambda history, residual norms history).

В следующем n обозначает размер матрицы и k количество требуемых собственных значений (наименьших или наибольших).

Код LOBPCG внутренне решает задачи на собственные значения размера 3k на каждой итерации, вызывая плотный решатель собственных значений eigh, поэтому если k не является достаточно малым по сравнению с n, нет смысла вызывать код LOBPCG. Более того, если вызвать алгоритм LOBPCG для 5k > n, это, вероятно, приведет к внутреннему сбою, поэтому код вызывает стандартную функцию eigh вместо. Не то, чтобы n должно быть большим для работы LOBPCG, а скорее соотношение n / k должно быть большим. Если вы вызываете LOBPCG с k=1 и n=10, он работает, хотя n мал. Метод предназначен для чрезвычайно больших n / k.

Скорость сходимости в основном зависит от трех факторов:

Качество начальных приближений X Элементы для слияния.
Относительное отделение желаемых собственных значений от остальных собственных значений. Можно варьировать k для улучшения разделения.
Правильное предобуславливание для уменьшения спектрального разброса. Например, задача теста вибрации стержня (в каталоге tests) плохо обусловлена для больших n, поэтому сходимость будет медленной, если не используется эффективное предобуславливание. Для этой конкретной задачи хорошей простой функцией предобуславливания будет линейное решение для A, что легко кодируется, поскольку A является трёхдиагональной.

Ссылки

[1]

A. V. Knyazev (2001), Toward the Optimal Preconditioned Eigensolver: Locally Optimal Block Preconditioned Conjugate Gradient Method. SIAM Journal on Scientific Computing 23, no. 2, pp. 517-541. DOI:10.1137/S1064827500366124

[2]

А. В. Князев, И. Лашук, М. Е. Аргентиати и Е. Овчинников (2007), Block Locally Optimal Preconditioned Eigenvalue Xolvers (BLOPEX) в hypre и PETSc. arXiv:0705.2626

[3]

Реализации A. V. Knyazev на C и MATLAB: lobpcg/blopex

Примеры

Наш первый пример минималистичен — найти наибольшее собственное значение диагональной матрицы, решая необобщённую задачу на собственные значения A x = lambda x без ограничений или предобусловливания.

>>> import numpy as np
>>> from scipy.sparse import diags_array
>>> from scipy.sparse.linalg import LinearOperator, aslinearoperator
>>> from scipy.sparse.linalg import lobpcg

Размер квадратной матрицы равен

>>> n = 100

и его диагональные элементы равны 1, …, 100, определенные как

>>> vals = np.arange(1, n + 1).astype(np.int16)

Первый обязательный входной параметр в этом тесте — разреженная диагональная матрица A задачи на собственные значения A x = lambda x для решения.

>>> A = diags_array(vals, offsets=0, shape=(n, n))
>>> A = A.astype(np.int16)
>>> A.toarray()
array([[  1,   0,   0, ...,   0,   0,   0],
       [  0,   2,   0, ...,   0,   0,   0],
       [  0,   0,   3, ...,   0,   0,   0],
       ...,
       [  0,   0,   0, ...,  98,   0,   0],
       [  0,   0,   0, ...,   0,  99,   0],
       [  0,   0,   0, ...,   0,   0, 100]], shape=(100, 100), dtype=int16)

Второй обязательный входной параметр X является двумерным массивом с размерностью строк, определяющей количество запрашиваемых собственных значений. X является начальным приближением для целевых собственных векторов. X должны иметь линейно независимые столбцы. Если начальные приближения недоступны, случайно ориентированные векторы обычно работают лучше, например, с компонентами, нормально распределёнными вокруг нуля или равномерно распределёнными на интервале [-1 1]. Установка начальных приближений в dtype np.float32 приводит все итеративные значения к dtype np.float32 ускорение выполнения при сохранении точности вычисления собственных значений.

>>> k = 1
>>> rng = np.random.default_rng()
>>> X = rng.normal(size=(n, k))
>>> X = X.astype(np.float32)

>>> eigenvalues, _ = lobpcg(A, X, maxiter=60)
>>> eigenvalues
array([100.], dtype=float32)

lobpcg требует только доступа к матричному произведению с A а не саму матрицу. Поскольку матрица A диагонален в этом примере, можно написать функцию матричного произведения A @ X используя диагональные значения vals только, например, поэлементным умножением с трансляцией в лямбда-функции

>>> A_lambda = lambda X: vals[:, np.newaxis] * X

или обычная функция

>>> def A_matmat(X):
...     return vals[:, np.newaxis] * X

и используйте указатель на один из этих вызываемых объектов в качестве входных данных

>>> eigenvalues, _ = lobpcg(A_lambda, X, maxiter=60)
>>> eigenvalues
array([100.], dtype=float32)
>>> eigenvalues, _ = lobpcg(A_matmat, X, maxiter=60)
>>> eigenvalues
array([100.], dtype=float32)

Традиционный вызываемый объект LinearOperator больше не требуется, но всё ещё поддерживается как вход для lobpcg. Указание matmat=A_matmat явно улучшает производительность.

>>> A_lo = LinearOperator((n, n), matvec=A_matmat, matmat=A_matmat, dtype=np.int16)
>>> eigenvalues, _ = lobpcg(A_lo, X, maxiter=80)
>>> eigenvalues
array([100.], dtype=float32)

Наименее эффективный вариант в виде вызываемого объекта — aslinearoperator:

>>> eigenvalues, _ = lobpcg(aslinearoperator(A), X, maxiter=80)
>>> eigenvalues
array([100.], dtype=float32)

Теперь переключаемся на вычисление трёх наименьших собственных значений, указывая

>>> k = 3
>>> X = np.random.default_rng().normal(size=(n, k))

и largest=False параметр

>>> eigenvalues, _ = lobpcg(A, X, largest=False, maxiter=90)
>>> print(eigenvalues)  
[1. 2. 3.]

Следующий пример иллюстрирует вычисление 3 наименьших собственных значений той же матрицы A заданный функциональным указателем A_matmat но с ограничениями и предобусловливанием.

Ограничения - необязательный входной параметр, представляющий собой 2D массив, состоящий из векторов-столбцов, которым собственные векторы должны быть ортогональны

>>> Y = np.eye(n, 3)

Предобуславливатель действует как обратный к A в этом примере, но в пониженной точности np.float32 хотя начальное X и, следовательно, все итерации и выходные данные находятся в полном np.float64.

>>> inv_vals = 1./vals
>>> inv_vals = inv_vals.astype(np.float32)
>>> M = lambda X: inv_vals[:, np.newaxis] * X

Давайте теперь решим задачу на собственные значения для матрицы A сначала без предобусловливания, запрашивая 80 итераций

>>> eigenvalues, _ = lobpcg(A_matmat, X, Y=Y, largest=False, maxiter=80)
>>> eigenvalues
array([4., 5., 6.])
>>> eigenvalues.dtype
dtype('float64')

С предобуславливанием нам нужно всего 20 итераций из тех же X

>>> eigenvalues, _ = lobpcg(A_matmat, X, Y=Y, M=M, largest=False, maxiter=20)
>>> eigenvalues
array([4., 5., 6.])

Обратите внимание, что векторы, передаваемые в Y являются собственными векторами 3 наименьших собственных значений. Возвращённые выше результаты ортогональны им.

Основная матрица A может быть неопределённой, например, после сдвига vals на 50 от 1, …, 100 до -49, …, 50, мы все еще можем вычислить 3 наименьших или наибольших собственных значения.

>>> vals = vals - 50
>>> X = rng.normal(size=(n, k))
>>> eigenvalues, _ = lobpcg(A_matmat, X, largest=False, maxiter=99)
>>> eigenvalues
array([-49., -48., -47.])
>>> eigenvalues, _ = lobpcg(A_matmat, X, largest=True, maxiter=99)
>>> eigenvalues
array([50., 49., 48.])