LOBPCG - LOBPCG

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

${ displaystyle Ax = lambda Bx,}$

для данной пары ${ Displaystyle (А, В)}$ сложных Эрмитский или реальный симметричный матрицы, где матрица ${ displaystyle B}$ также предполагается положительно определенный.

Фон

Канторович в 1948 г. предложил рассчитывать наименьшие собственное значение ${ displaystyle lambda _ {1}}$ симметричной матрицы ${ displaystyle A}$ к крутой спуск используя направление ${ displaystyle r = Ax- lambda (x) x}$ масштабного градиент из Фактор Рэлея ${ Displaystyle лямбда (х) = (х, Ах) / (х, х)}$ в скалярное произведение ${ Displaystyle (х, у) = х'у}$, с размером шага, вычисленным путем минимизации частного Рэлея в линейный пролет векторов ${ displaystyle x}$ и ${ displaystyle w}$, т.е. локально оптимальным образом. Самокиш^[1] предложил применить предварительный кондиционер ${ displaystyle T}$ к остаточному вектору ${ displaystyle r}$ для генерации предварительно обусловленного направления ${ displaystyle w = Tr}$ и получили асимптотику, как ${ displaystyle x}$ приближается к собственный вектор, границы скорости сходимости. Дьяконов предложенный^[2] спектрально эквивалентный предварительная подготовка и получены неасимптотические границы скорости сходимости. Блочный локально оптимальный многошаговый наискорейший спуск для задач на собственные значения описан в.^[3] Локальная минимизация отношения Рэлея на подпространстве, охватываемом текущим приближением, текущим остатком и предыдущим приближением, а также его блочной версией, появилась в.^[4] Предварительно подготовленная версия была проанализирована в ^[5] и.^[6]

Основные особенности^[7]

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

Алгоритм

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

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

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

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

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

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

${ displaystyle r: = Ax- rho (x) Bx,}$

называется собственным вектором остаточный. Если предварительный кондиционер ${ displaystyle T}$ доступен, он применяется к невязке и дает вектор

${ displaystyle w: = Tr,}$

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

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

или, короче,

${ Displaystyle х ^ {я + 1}: = х ^ {я} + альфа ^ {я} ш ^ {я}, ,}$

${ Displaystyle ш ^ {я}: = Тр ^ {я}, ,}$

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

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

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

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

Трехместное повторение

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

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

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

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

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

Аналоги подпространств Крылова

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

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

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

Версия блока

Резюме

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

Основной дизайн

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

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

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

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

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

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

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

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

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

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

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

Частичное Анализ главных компонентов (PCA) и Разложение по сингулярным значениям (СВД)

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

LOBPCG для PCA и SVD реализован в SciPy, начиная с версии 1.4.0.^[12]

Общие программные реализации

Изобретатель LOBPCG, Андрей Князев, опубликовал эталонную реализацию под названием Block Locally Optimal Preconditioned Eigenvalue Xolvers (BLOPEX)^[13][14] с интерфейсами к PETSc, гипре, и параллельный иерархический адаптивный многоуровневый метод (PHAML).^[15] Другие реализации доступны, например, в GNU Octave,^[16] MATLAB (в том числе для распределенных или мозаичных массивов),^[9] Ява,^[17] Анасази (Трилинос ),^[18] SLEPc,^[19][20] SciPy,^[10] Юля,^[21] МАГМА,^[22] Pytorch,^[23] Ржавчина,^[24] OpenMP и OpenACC,^[25] RAPIDS cuGraph^[26] и NVIDIA AMGX.^[27] LOBPCG реализован,^[28] но не входит в TensorFlow.

Приложения

Материаловедение

LOBPCG реализован в ABINIT^[29] (включая CUDA версия) и Осьминог.^[30] Он был использован для многомиллиардных матриц размером Приз Гордона Белла финалистов, на Симулятор Земли суперкомпьютер в Японии.^[31][32] Модель Хаббарда для сильно коррелированных электронных систем, чтобы понять механизм, лежащий в основе сверхпроводимость использует LOBPCG для расчета основное состояние из Гамильтониан на K компьютер.^[33] Есть MATLAB ^[34] и Юля^[35][36][37]версии LOBPCG для Кон-Шам уравнения и теория функционала плотности (DFT) с использованием простой волны. Недавние реализации включают TTPY,^[38] Утконос ‐ QM,^[39] МФДн,^[40] АПФ-молекула,^[41] LACONIC.^[42]

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

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

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

LOBPCG - один из основных решателей собственных значений в PYFEMax и высокопроизводительной мультифизике. заключительный элемент программное обеспечение Netgen / NGSolve. LOBPCG из гипре включен в Открытый исходный код легкий масштабируемый C ++ библиотека для заключительный элемент методы MFEM, который используется во многих проектах, в том числе ВЗРЫВ, XBraid, Посещение, xSDK, институт FASTMath в SciDAC, и Центр совместного проектирования по эффективной экзадачной дискретизации (CEED) в Exascale вычисления Проект.

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

Итеративная аппроксимация на основе LOBPCG фильтр нижних частот может использоваться для шумоподавление; видеть,^[44] например, чтобы ускорить полное изменение шумоподавления.

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

Сегментация изображения через спектральная кластеризация выполняет низкоразмерный встраивание используя близость матрица между пикселями с последующей кластеризацией компонентов собственных векторов в низкоразмерном пространстве. LOBPCG с многосеточный предварительная подготовка был впервые применен к сегментация изображения в ^[45] через спектральный разбиение графа с использованием граф лапласиан для двусторонний фильтр. Scikit-Learn использует LOBPCG из SciPy с алгебраическое многосеточное предварительное кондиционирование для решения задачи на собственные значения.^[46]

Сбор данных

Программные пакеты scikit-learn и Мегамен^[47] используйте LOBPCG для масштабирования спектральная кластеризация^[48] и многообразное обучение^[49] через Собственные карты лапласиана для больших наборов данных. NVIDIA реализовал^[50] LOBPCG в своей библиотеке nvGRAPH, представленной в CUDA 8.

Рекомендации

внешняя ссылка

• LOBPCG в MATLAB
• LOBPCG в Октава
• LOBPCG в SciPy
• LOBPCG в Ява в Код Google
• LOBPCG в блоке локально оптимальных предобусловленных Xolvers собственных значений (BLOPEX) в GitHub и в архиве в Код Google