ЛОБПКГ

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

$Ax=\lambda Bx,$

для данной пары $(A,B)$ комплексных эрмитовых или действительных симметричных матриц, гдематрица $B$ также предполагается положительно-определенным .

Фон

Канторович в 1948 году предложил вычислять наименьшее собственное значение. $\lambda _{1}$ симметричной матрицы $A$ по наискорейшему спуску по направлению $r=Ax-\lambda (x)x$ масштабированного градиента коэффициента Рэлея $\lambda (x)=(x,Ax)/(x,x)$ в скалярном произведении $(x,y)=x'y$, с размером шага, вычисляемым путем минимизации коэффициента Рэлея в линейной области векторов $x$ и $w$, то есть локально оптимальным образом. Самокиш предложил применить предобуславливатель $T$ вектору остатков $r$ для создания заранее обусловленного направления $w=Tr$ и выведена асимптотика, как $x$ приближается к собственному вектору , границам скорости сходимости. Дьяконов предложил спектрально эквивалентная предобусловливание и полученные неасимптотические границы скорости сходимости. Блочный локально оптимальный многошаговый наискорейший спуск для задач на собственные значения описан в . В 2007 году появилась локальная минимизация фактора Рэлея на подпространстве, натянутом текущим приближением, текущей невязкой и предыдущим приближением, а также его блочный вариант. Предварительно обусловленная версия была проанализирована в и.

Основные особенности

• Безматрицный , т.е. не требует явного хранения матрицы коэффициентов, но может получить доступ к матрице путем оценки произведений матрицы-вектора.
• Не требует факторизации , т.е. не требует какого-либо матричного разложения даже для обобщенной задачи на собственные значения .
• Затраты на итерацию и использование памяти конкурентоспособны по сравнению с методом Ланцоша , вычисляющим одну крайнюю собственную пару симметричной матрицы.
• Линейная сходимость теоретически гарантирована и практически наблюдается.
• Ускоренная сходимость за счет прямой предобуславливания , в отличие от метода Ланцоша , включая переменную и несимметричную, а также фиксированную и положительно определенную предобусловливание .
• Позволяет тривиально использовать эффективную декомпозицию домена и многосеточные методы посредством предварительной подготовки.
• Теплый старт и вычисление аппроксимации собственного вектора на каждой итерации.
• Более численно стабилен по сравнению с методом Ланцоша и может работать в компьютерной арифметике низкой точности.
• Легко реализовать, уже появилось много версий.
• Блокировка позволяет использовать высокоэффективные матрично-матричные операции, например BLAS 3.
• Размер блока можно настроить, чтобы сбалансировать скорость сходимости и компьютерные затраты на ортогонализацию и метод Рэлея-Ритца на каждой итерации.

Алгоритм

Одновекторная версия

Предварительные сведения: Градиентный спуск для задач на собственные значения

Метод выполняет итеративную максимизацию (или минимизацию) обобщенного коэффициента Рэлея.

$\rho (x):=\rho (A,B;x):={\frac {x^{T}Ax}{x^{T}Bx}},$

что приводит к нахождению наибольших (или наименьших) собственных пар $Ax=\lambda Bx.$

Направление наибольшего подъема, которое представляет собой градиент обобщенного фактора Рэлея , положительно пропорционально вектору

$r:=Ax-\rho (x)Bx,$

называется остатком собственного вектора . Если предобуславливатель $T$ доступен, он применяется к остатку и дает вектор

$w:=Tr,$

называется предобусловленным остатком. Без предобуславливания положим $T:=I$ и так $w:=r$. Итерационный метод

$x^{i+1}:=x^{i}+\alpha ^{i}T(Ax^{i}-\rho (x^{i})Bx^{i}),$

или, короче говоря,

$x^{i+1}:=x^{i}+\alpha ^{i}w^{i},\,$

$w^{i}:=Tr^{i},\,$

$r^{i}:=Ax^{i}-\rho (x^{i})Bx^{i},$

известен как предварительно обусловленный крутой подъем (или спуск), где скаляр $\alpha ^{i}$ называется размером шага. Оптимальный размер шага можно определить путем максимизации коэффициента Рэлея, т.е.

$x^{i+1}:=\arg \max _{y\in span\{x^{i},w^{i}\}}\rho (y)$

(или $\arg \min$ в случае минимизации), в этом случае метод называется локально оптимальным.

Трехкратный рецидив

Чтобы значительно ускорить сходимость локально оптимального заранее обусловленного наискорейшего подъема (или спуска), к двухчленному рекуррентному соотношению можно добавить один дополнительный вектор , чтобы сделать его трехчленным:

$x^{i+1}:=\arg \max _{y\in span\{x^{i},w^{i},x^{i-1}\}}\rho (y)$

(использовать $\arg \min$ в случае минимизации). Максимизация/минимизация фактора Рэлея в трехмерном подпространстве может быть выполнена численно с помощью метода Рэлея–Ритца . Добавление дополнительных векторов, см., например, экстраполяцию Ричардсона , не приводит к значительному ускорению. но увеличивает затраты на вычисления, поэтому обычно не рекомендуется.

Улучшения численной стабильности

По мере сходимости итераций векторы $x^{i}$ и $x^{i-1}$ становятся почти линейно зависимыми , что приводит к потере точности и делает метод Рэлея-Ритца численно нестабильным при наличии ошибок округления. Потери точности можно избежать, заменив вектор $x^{i-1}$ с вектором $p^{i}$, который может быть дальше от $x^{i}$, в базисе трехмерного подпространства $span\{x^{i},w^{i},x^{i-1}\}$, сохраняя при этом подпространство неизменным и избегая ортогонализации или любых других дополнительных операций. Кроме того, ортогонализация базиса трехмерного подпространства может потребоваться для плохо обусловленных задач на собственные значения для повышения устойчивости и достижимой точности.

Krylov subspace analogs

Это одновекторная версия метода LOBPCG — одно из возможных обобщений линейных решателей с предварительно обусловленными сопряженными градиентами на случай симметричных задач на собственные значения . Даже в тривиальном случае $T=I$ и $B=I$ полученное приближение с $i>3$ будет отличаться от полученного алгоритмом Ланцоша , хотя обе аппроксимации будут принадлежать одному и тому же подпространству Крылова .

Сценарии практического использования

Чрезвычайная простота и высокая эффективность одновекторной версии LOBPCG делают ее привлекательной для приложений, связанных с собственными значениями, в условиях жестких аппаратных ограничений, начиная от спектральной кластеризации в реальном времени на основе обнаружения аномалий через разделение графа на встроенных ASIC или FPGA до моделирования физических явлений записи. сложность вычислений на экзафлопсных суперкомпьютерах TOP500 .

Блок версия

Краткое содержание

Последующие собственные пары могут быть вычислены одна за другой с помощью одновекторного LOBPCG, дополненного ортогональным дефляцией, или одновременно как блок. В первом подходе неточности в уже вычисленных приближенных собственных векторах аддитивно влияют на точность вычисленных впоследствии собственных векторов, тем самым увеличивая ошибку с каждым новым вычислением. Итерация нескольких приближенных собственных векторов вместе в блоке локально оптимальным способом в блочной версии LOBPCG. позволяет быстро, точно и надежно вычислить собственные векторы, в том числе те, которые соответствуют почти множеству собственных значений, где одновекторный LOBPCG страдает от медленной сходимости. Размер блока можно настроить, чтобы сбалансировать численную стабильность, скорость сходимости и компьютерные затраты на ортогонализацию и метод Рэлея-Ритца на каждой итерации.

Основная конструкция

Блочный подход в LOBPCG заменяет одновекторный подход. $x^{i},\,w^{i},$ и $p^{i}$ с блок-векторами, т.е. матрицами $X^{i},\,W^{i},$ и $P^{i}$, где, например, каждый столбец $X^{i}$ аппроксимирует один из собственных векторов. Все столбцы повторяются одновременно, и следующая матрица приближенных собственных векторов $X^{i+1}$ определяется методом Рэлея–Ритца на подпространстве, натянутом на все столбцы матриц $X^{i},\,W^{i},$ и $P^{i}$. Каждый столбец $W^{i}$ вычисляется просто как предварительно обусловленный остаток для каждого столбца $X^{i}.$ Матрица $P^{i}$ определяется так, что подпространства, охватываемые столбцами $[X^{i},\,P^{i}]$ и из $[X^{i},\,X^{i-1}]$ одинаковы.

Численная стабильность против эффективности

Результат метода Рэлея – Ритца определяется подпространством, охватываемым всеми столбцами матриц. $X^{i},\,W^{i},$ и $P^{i}$, где базис подпространства теоретически может быть произвольным. Однако в неточной компьютерной арифметике метод Рэлея-Ритца становится численно неустойчивым, если некоторые из базисных векторов приблизительно линейно зависимы. Численная нестабильность обычно возникает, например, если некоторые из собственных векторов в итерационном блоке уже достигают достижимой точности для заданной компьютерной точности и особенно заметны при низкой точности, например, одинарной точности .

Искусство множественной реализации LOBPCG заключается в обеспечении численной устойчивости метода Рэлея-Ритца при минимальных вычислительных затратах за счет выбора хорошего базиса подпространства. Вероятно, наиболее стабильный подход к ортогонализации базисных векторов, например, с помощью процесса Грама – Шмидта , также является самым дорогостоящим в вычислительном отношении. Например, реализации LOBPCG, использовать нестабильное, но эффективное разложение Холецкого нормальной матрицы , которое выполняется только на отдельных матрицах $W^{i}$ и $P^{i}$, а не на всем подпространстве. Постоянно увеличивающийся объем компьютерной памяти позволяет в настоящее время использовать типичные размеры блоков. $10^{3}-10^{4}$ диапазон, в котором процент вычислительного времени, затрачиваемый на ортогонализацию, и метод Рэлея-Ритца начинает доминировать.

Блокировка ранее сошедшихся собственных векторов

Блочные методы для задач собственных значений, которые перебирают подпространства, обычно имеют некоторые из итеративных собственных векторов, сходящихся быстрее, чем другие, что мотивирует блокировку уже сошедшихся собственных векторов, т. е. удаление их из итеративного цикла, чтобы исключить ненужные вычисления и улучшить численную стабильность. Простое удаление собственного вектора, вероятно, может привести к образованию его дубликата в повторяющихся векторах. Тот факт, что собственные векторы симметричных задач на собственные значения попарно ортогональны, предполагает сохранение всех итерационных векторов ортогональными заблокированным векторам.

Блокировку можно реализовать по-разному, сохраняя числовую точность и стабильность при минимизации вычислительных затрат. Например, реализации LOBPCG, следовать, отделение жесткой блокировки, т. е. дефляции путем ограничения, когда заблокированные собственные векторы служат входными данными кода и не изменяются, от мягкой блокировки, когда заблокированные векторы не участвуют в обычно самом дорогостоящем итеративном этапе вычисления остатков, однако полностью участвуют в методе Рэлея-Ритца и, таким образом, могут быть изменены методом Рэлея-Ритца.

Модификации, LOBPCG II

LOBPCG включает все столбцы матриц. $X^{i},\,W^{i},$ и $P^{i}$ в метод Рэлея-Ритца, что приводит к $3k$-к- $3k$ проблему собственных значений, которую необходимо решить, и до $9k^{2}$ скалярные произведения для вычисления на каждой итерации, где $k$ обозначает размер блока — количество столбцов. Для блоков большого размера $k$ это начинает преобладать над затратами на вычисления и ввод-вывод и ограничивать распараллеливание, когда несколько вычислительных устройств работают одновременно.

Оригинальная бумага LOBPCG. описывает модификацию, названную LOBPCG II, для решения такой проблемы, использующую одновекторную версию метода LOBPCG для каждой желаемой собственной пары с решением процедуры Рэлея-Ритца. $k$ прогнозируемых проблем собственных значений 3х3. Глобальная процедура Рэлея-Ритца для всех $k$ собственные пары есть на каждой итерации, но только в столбцах матрицы $X^{i}$, тем самым уменьшая количество необходимых скалярных произведений до $k^{2}+3k$ от $9k^{2}$ и размер глобальной прогнозируемой проблемы собственных значений $k$-к- $k$ от $3k$-к- $3k$ на каждой итерации. Ссылка идет дальше, применяя алгоритм LOBPCG к каждому аппроксимированному собственному вектору отдельно, т. е. запуская разблокированную версию метода LOBPCG для каждой желаемой собственной пары в течение фиксированного числа итераций. Процедуры Рэлея-Ритца в этих запусках должны решить только набор задач на проецируемые собственные значения 3 × 3. Глобальная процедура Рэлея-Ритца для всех желаемых собственных пар применяется только периодически в конце фиксированного числа разблокированных итераций LOBPCG.

Такие модификации могут быть менее надежными по сравнению с исходным LOBPCG. Индивидуально выполняемые ветви одновекторной LOBPCG не могут следовать непрерывным итеративным путям, вместо этого переворачиваясь и создавая дублированные аппроксимации одного и того же собственного вектора. Одновекторная LOBPCG может быть непригодна для кластеризованных собственных значений, но отдельные прогоны LOBPCG с малыми блоками требуют автоматического определения размеров их блоков в процессе итераций, поскольку количество кластеров собственных значений и их размеры могут быть априори неизвестны.

Теория и практика конвергенции

LOBPCG по конструкции гарантируется минимизировать коэффициент Рэлея не медленнее, чем блок наикрутейшего градиентного спуска , который имеет всеобъемлющую теорию сходимости. Каждый собственный вектор является стационарной точкой фактора Рэлея , где градиент исчезает. Таким образом, градиентный спуск может замедлиться в окрестности любого собственного вектора , однако он гарантированно либо сходится к собственному вектору с линейной скоростью сходимости, либо, если этот собственный вектор является седловой точкой , итеративный коэффициент Рэлея с большей вероятностью упадет. ниже соответствующего собственного значения и начните линейно сходиться к следующему собственному значению ниже. Определено наихудшее значение скорости линейной сходимости. и зависит от относительного разрыва между собственным значением и остальной частью матричного спектра и качества предобуславливателя , если он присутствует.

Для общей матрицы, очевидно, не существует способа предсказать собственные векторы и, таким образом, сгенерировать начальные приближения, которые всегда работают хорошо. Итерационное решение с помощью LOBPCG может быть чувствительным к начальным аппроксимациям собственных векторов, например, для сходимости требуется больше времени, а при прохождении промежуточных собственных пар оно замедляется. Более того, теоретически нельзя гарантировать сходимость обязательно к наименьшей собственной паре, хотя вероятность промаха равна нулю. хорошего качества Случайная функция Гаусса с нулевым средним значением обычно используется в LOBPCG по умолчанию для генерации начальных приближений. Чтобы зафиксировать начальные приближения, можно выбрать фиксированное начальное число для генератора случайных чисел .

В отличие от метода Ланцоша редко демонстрирует асимптотическую суперлинейную сходимость , LOBPCG на практике .

Частичный анализ главных компонент (PCA) и разложение по сингулярным значениям (SVD)

LOBPCG можно тривиально адаптировать для вычисления нескольких наибольших сингулярных значений и соответствующих сингулярных векторов (частичных SVD), например, для итеративного вычисления PCA для матрицы данных D с нулевым средним значением, без явного вычисления ковариационной матрицы D. ^Т D , то есть безматрицным способом . Основным расчетом является оценка функции произведения D ^Т (DX) ковариационной матрицы D ^Т D и блочный вектор X , который итеративно аппроксимирует нужные сингулярные векторы. PCA требует наибольших собственных значений ковариационной матрицы, тогда как LOBPCG обычно реализуется для расчета наименьших. Простой обходной путь — отменить функцию, заменив -D ^Т (DX) для D ^Т (DX) и, таким образом, меняя порядок собственных значений на противоположный, поскольку LOBPCG не заботится о том, является ли матрица задачи собственных значений положительно определенной или нет.

LOBPCG для PCA и SVD реализован в SciPy начиная с версии 1.4.0.

Общие реализации программного обеспечения

Изобретатель LOBPCG, Андрей Князев , опубликовал эталонную реализацию под названием Block Locally Optimal Preconditional Eigenvalue Xolvers (BLOPEX). с интерфейсами к PETSc , hypre и параллельному иерархическому адаптивному многоуровневому методу (PHAML). Другие реализации доступны, например, в GNU Octave . MATLAB (в том числе для распределенных или мозаичных массивов), Ява , Анасази ( Трилинос ), СЛЕПк , SciPy , Юлия , МАГМА, Пайторч , Ржавчина , OpenMP и OpenACC , CuPy ( NumPy -совместимая библиотека массивов, ускоряемая с помощью CUDA ), Гугл Джакс , и NVIDIA AMGX. LOBPCG реализован, но не включен в TensorFlow .

Приложения

Интеллектуальный анализ данных

Пакеты программного обеспечения scikit-learn и Megaman используйте LOBPCG для масштабирования спектральной кластеризации и разнообразное обучение через собственные карты Лапласа к большим наборам данных. NVIDIA реализовала LOBPCG в своей библиотеке nvGRAPH, представленной в CUDA 8. Sphynx, гибридный параллельный разделитель графов с распределенной и общей памятью — первый инструмент разделения графов, который работает на графических процессорах с настройками распределенной памяти — использует спектральную кластеризацию для разделения графа , вычисляя собственные векторы на матрице Лапласа графа с использованием LOBPCG из Анасази упаковка.

Материальные науки

LOBPCG реализован в ABINIT. (включая CUDA версию ) и Octopus . использовали его для создания матриц размером в несколько миллиардов Финалисты премии Гордона Белла на Earth Simulator суперкомпьютере в Японии. Модель Хаббарда для сильно коррелированных электронных систем для понимания механизма сверхпроводимости использует LOBPCG для расчета основного состояния гамильтониана . на K-компьютере и системы с несколькими графическими процессорами.

Есть МАТЛАБ и Джулия версии LOBPCG для уравнений Кона-Шэма и теории функционала плотности (DFT) с использованием базиса плоских волн.Последние реализации включают TTPY, Утконос-QM, МФД, ACE-Молекула, ЛАКОНИЧНЫЙ.

Механика и жидкости

LOBPCG от BLOPEX используется для настройки прекондиционера в библиотеке решателей многоуровневого балансирующего доменного разложения по ограничениям (BDDC) BDDCML, которая включена в OpenFTL (открытая библиотека шаблонов конечных элементов ) и симулятор Flow123d потока подземных вод, переноса растворенных веществ и тепла в трещиноватых пористых средах. . LOBPCG реализован в LS-DYNA и косвенно в ANSYS .

Уравнения Максвелла

LOBPCG — один из основных решателей собственных значений в PYFEMax и высокопроизводительном мультифизическом программном обеспечении конечных элементов Netgen/NGSolve. LOBPCG от hypre включен в с открытым исходным кодом облегченную масштабируемую C++ библиотеку конечных элементов для методов MFEM , которая используется во многих проектах, включая BLAST , XBraid, VisIt , xSDK, институт FASTMath в SciDAC и Центр совместного проектирования эффективной экзафлопсной дискретизации. (CEED) в проекте экзафлопсных вычислений .

шумоподавление

на основе LOBPCG итеративный приблизительный фильтр нижних частот можно использовать Для шумоподавления ; видеть, например, для ускорения общего шумоподавления .

Сегментация изображений

Сегментация изображения посредством спектральной кластеризации выполняет низкоразмерное встраивание с использованием матрицы сродства между пикселями с последующей кластеризацией компонентов собственных векторов в низкомерном пространстве, например, с использованием графа Лапласа для двустороннего фильтра . Сегментация изображений спектральных посредством разделения графов с помощью LOBPCG с многосеточной предварительной обработкой была впервые предложена в и фактически протестировано в и. Последний подход был позже реализован в Python scikit-learn. который использует LOBPCG из SciPy с алгебраической многосеточной предварительной обуславливанием для решения проблемы собственных значений для лапласиана графа.

Внешние ссылки

• LOBPCG в MATLAB
• ЛОБПКГ в Октаве
• LOBPCG в SciPy
• LOBPCG на Java в Google Code
• LOBPCG в блоке локально оптимальных предварительно подготовленных Xolvers собственных значений (BLOPEX) на GitHub и заархивирован в Google Code.