В математике обобщенный метод минимальный остаточный (GMRES) представляет собой итерационный метод для численного решения неопределенной несимметричной системы линейных уравнений . Метод аппроксимирует решение вектором из подпространства Крылова с минимальной невязкой . Для нахождения этого вектора используется итерация Арнольди .
Метод GMRES был разработан Юсефом Саадом и Мартином Х. Шульцем в 1986 году. [1] GMRES является обобщением и улучшением метода MINRES , в то время как последний также требует, чтобы матрица была симметричной, но имеет то преимущество, что требует только обработки трех векторов. [2] [3] Он был разработан Крисом Пейдж и Майкл Сондерс в 1975 г. GMRES также является частным случай ДИМЕЙ метод , разработанных Питером Пули в 1980 году DIIS также применим к нелинейным системам.
Метод
Обозначим евклидову норму любого вектора v через. Обозначим (квадратную) систему линейных уравнений, которую необходимо решить, через
Матрица считаются обратимым размером м матрица с размерностью м . Кроме того, предполагается, что b нормировано, т. Е. Что.
В н -го Крылов подпространство для этой задачи
где это начальная ошибка при первоначальном предположении . Четко если .
GMRES приближает точное решение по вектору минимизирующий евклидову норму невязки .
Векторы может быть близка к линейно зависимой , поэтому вместо этого базиса используется итерация Арнольди для поиска ортонормированных векторов которые составляют основу . В частности,.
Следовательно, вектор можно записать как с участием , где представляет собой M матрицу с размерностью п матрица , составленная.
Процесс Арнольди также дает ()-от-верхняя матрица Гессенберга с участием
Для симметричных матриц фактически получается симметричная трехдиагональная матрица, что приводит к методу Минреса .
Поскольку столбцы ортонормированы, мы имеем
где
первый вектор в стандартном базисе из, а также
являющийся первым пробным вектором (обычно нулевым). Следовательно, можно найти, минимизируя евклидову норму невязки
Это линейная задача наименьших квадратов размера n .
Это дает метод GMRES. На-я итерация:
- вычислить по методу Арнольди;
- Найти что сводит к минимуму ;
- вычислить ;
- повторите, если остаток еще недостаточно мал.
На каждой итерации произведение матрица-вектор должны быть вычислены. Это стоит около операции с плавающей запятой для общих плотных матриц размера, но стоимость может снизиться до для разреженных матриц . В дополнение к произведению матрица-вектор,операции с плавающей точкой должны вычисляться на n -й итерации.
Конвергенция
П - й итерации минимизирует невязку в подпространстве Крылова. Поскольку каждое подпространство содержится в следующем подпространстве, невязка не увеличивается. После m итераций, где m - размер матрицы A , пространство Крылова K m - это все R m, и, следовательно, метод GMRES приходит к точному решению. Однако идея состоит в том, что после небольшого количества итераций (относительно m ) вектор x n уже является хорошим приближением к точному решению.
Такого не бывает вообще. Действительно, теорема Гринбаума, Патака и Стракоша утверждает, что для любой невозрастающей последовательности a 1 ,…, a m −1 , a m = 0 можно найти такую матрицу A , что || r n || = a n для всех n , где r n - невязка, определенная выше. В частности, можно найти матрицу, для которой невязка остается постоянной в течение m - 1 итераций и падает до нуля только на последней итерации.
Однако на практике GMRES часто работает хорошо. Это можно доказать в конкретных ситуациях. Если симметричная часть A , то есть, положительно определена , то
где а также обозначают наименьшее и наибольшее собственное значение матрицы, соответственно. [4]
Если является симметричной и положительно определена, то у нас даже есть
где обозначает число обусловленности из А в евклидовой норме.
В общем случае, когда A не является положительно определенным, имеем
где Р п обозначает множество многочленов степени не выше п с р (0) = 1, V является матрицей , появляющаяся в спектральном разложении из A , и σ ( ) является спектр из A . Грубо говоря, это говорит о том, что быстрая сходимость происходит, когда собственные значения A сгруппированы далеко от начала координат, а A не слишком далеко от нормальности . [5]
Все эти неравенства ограничивают только остатки, а не фактическую ошибку, то есть расстояние между текущей итерацией x n и точным решением.
Расширения метода
Как и другие итерационные методы, GMRES обычно сочетается с методом предварительной обработки , чтобы ускорить сходимость.
Стоимость итераций растет как O ( n 2 ), где n - номер итерации. Поэтому метод иногда перезапускается после некоторого количества, скажем k , итераций с x k в качестве начального предположения. Полученный метод называется GMRES ( k ) или Restarted GMRES. Для неположительно определенных матриц этот метод может страдать от стагнации сходимости, поскольку перезапущенное подпространство часто близко к предыдущему подпространству.
Недостатки GMRES и перезапущенного GMRES устраняются повторным использованием подпространства Крылова в методах типа GCRO, таких как GCROT и GCRODR. [6] Повторное использование подпространств Крылова в GMRES также может ускорить сходимость, когда необходимо решить последовательности линейных систем. [7]
Сравнение с другими решателями
Итерация Арнольди сводится к итерации Ланцоша для симметричных матриц. Соответствующий метод подпространств Крылова является методом минимальных невязок (MinRes) Пейдж и Сондерс. В отличие от несимметричного случая, метод MinRes задается трехчленным рекуррентным соотношением . Можно показать, что не существует метода подпространств Крылова для общих матриц, который задается коротким рекуррентным соотношением и все же минимизирует нормы остатков, как это делает GMRES.
Другой класс методов основан на несимметричной итерации Ланцоша , в частности метод BiCG . Они используют трехчленное рекуррентное соотношение, но они не достигают минимальной невязки, и, следовательно, невязка не уменьшается монотонно для этих методов. Сходимость даже не гарантируется.
Третий класс образован такими методами, как CGS и BiCGSTAB . Они также работают с трехчленным рекуррентным отношением (следовательно, без оптимальности) и могут даже преждевременно завершиться без достижения сходимости. Идея этих методов заключается в подходящем выборе порождающих полиномов итерационной последовательности.
Ни один из этих трех классов не является лучшим для всех матриц; всегда есть примеры, в которых один класс превосходит другой. Поэтому на практике используются несколько решателей, чтобы увидеть, какой из них лучше всего подходит для данной проблемы.
Решение задачи наименьших квадратов
Одна часть метода GMRES - найти вектор что сводит к минимуму
Обратите внимание, что является матрицей размером ( n + 1) на n , следовательно, она дает линейную систему из n +1 уравнений с избыточными ограничениями для n неизвестных.
Минимум может быть вычислен с помощью QR-разложения : найти ( n + 1) -by- ( n + 1) ортогональную матрицу Ω n и ( n + 1) -by- n матрицу верхнего треугольника такой, что
Треугольная матрица имеет на одну строку больше, чем столбцов, поэтому ее нижняя строка состоит из нуля. Следовательно, его можно разложить как
где является н - матрицу с размерностью п (таким образом , квадрат) треугольная матрица.
QR-разложение можно дешево обновлять от одной итерации к следующей, потому что матрицы Хессенберга различаются только строкой нулей и столбцом:
где ч п + 1 = ( ч 1, п + 1 , ..., ч п + 1, п + 1 ) Т . Это означает, что предварительное умножение матрицы Хессенберга на Ω n , дополненное нулями и строкой с мультипликативной единицей, дает почти треугольную матрицу:
Это было бы треугольником, если σ равно нулю. Чтобы исправить это, нужно вращение Гивенса.
где
С помощью этого вращения Гивенса мы формируем
Действительно,
представляет собой треугольную матрицу.
Учитывая QR-разложение, проблема минимизации легко решается, если заметить, что
Обозначая вектор от
при g n ∈ R n и γ n ∈ R это
Вектор y, который минимизирует это выражение, задается следующим образом:
И снова векторы легко обновляются. [8]
Пример кода
Обычный GMRES (MATLAB / GNU Octave)
функция [x, e] = gmres ( A, b, x, max_iterations, threshold ) n = длина ( A ); m = max_iterations ; % использовать x как начальный вектор г = Ь - А * х ; b_norm = норма ( б ); ошибка = norm ( r ) / b_norm ; % инициализировать одномерные векторы sn = нули ( m , 1 ); cs = нули ( m , 1 ); % e1 = нули (n, 1); e1 = нули ( m + 1 , 1 ); е1 ( 1 ) = 1 ; е = [ ошибка ]; r_norm = норма ( г ); Q (:, 1 ) = r / r_norm ; бета = r_norm * e1 ; % Примечание: это не бета-скаляр из раздела «Метод» выше, а бета-скаляр, умноженный на e1 для k = 1 : m % беги арнольди [ H ( 1 : k + 1 , k ) Q (:, k + 1 )] = arnoldi ( A , Q , k ); % удалить последний элемент в строке H и обновить матрицу вращения [ H ( 1 : k + 1 , k ) cs ( k ) sn ( k )] = apply_givens_rotation ( H ( 1 : k + 1 , k ), cs , sn , k ); % обновить остаточный вектор бета ( к + 1 ) = - sn ( к ) * бета ( к ); бета ( к ) = сс ( к ) * бета ( к ); error = abs ( beta ( k + 1 )) / b_norm ; % сохранить ошибку е = [ е ; ошибка ]; если ( ошибка <= порог ) перерыв ; конец конец %, если порог не достигнут, в этой точке k = m (а не m + 1) % рассчитать результат у = Н ( 1 : к , 1 : к ) \ бета ( 1 : к ); х = х + Q (:, 1 : к ) * у ; конец% ------------------------------------------------- ---%% Функция Арнольди%% ------------------------------------------------- ---%функция [h, q] = arnoldi ( A, Q, k ) q = A * Q (:, k ); % Крылова Вектор для i = 1 : k % модифицированного Грама-Шмидта, сохраняя матрицу Гессенберга h ( i ) = q ' * Q (:, i ); q = q - h ( i ) * Q (:, я ); конец h ( k + 1 ) = norm ( q ); q = q / h ( k + 1 ); конец% ------------------------------------------------- --------------------%% Применение вращения Гивенса к H col%% ------------------------------------------------- --------------------%функция [h, cs_k, sn_k] = apply_givens_rotation ( h, cs, sn, k ) % подать заявку на i-й столбец для i = 1 : k - 1 темп = cs ( i ) * h ( i ) + sn ( i ) * h ( i + 1 ); h ( i + 1 ) = - sn ( i ) * h ( i ) + cs ( i ) * h ( i + 1 ); h ( i ) = темп ; конец % обновить следующие значения sin cos для вращения [ cs_k sn_k ] = givens_rotation ( h ( k ), h ( k + 1 )); % удалить H (i + 1, i) h ( k ) = cs_k * h ( k ) + sn_k * h ( k + 1 ); h ( k + 1 ) = 0,0 ; конец%% ---- Вычислить заданную матрицу вращения ---- %%функция [cs, sn] = givens_rotation ( v1, v2 ) % если (v1 == 0)% cs = 0;% sn = 1;% еще т = sqrt ( v1 ^ 2 + v2 ^ 2 ); % cs = abs (v1) / t;% sn = cs * v2 / v1; cs = v1 / t ; % см. http://www.netlib.org/eispack/comqr.f sn = v2 / t ; % конецконец
Смотрите также
- Метод двусопряженного градиента
Рекомендации
- ^ Y. Саад и М. Х. Шульц
- ^ Пейдж и Сондерс, "Решение разреженных неопределенных систем линейных уравнений", SIAM J. Numer. Anal., Том 12, стр. 617 (1975) https://doi.org/10.1137/0712047
- ^ Н.Нифа. «Докторская диссертация» (PDF) .
- ^ Eisenstat, Elman & Schultz, Thm 3.3. NB: все результаты для GCR справедливы и для GMRES, см. Саад и Шульц
- ^ Trefethen & Bau, Thm 35.2
- ^ Амриткар, Амит; де Стерлер, Эрик; Свиридович, Катажина; Тафти, Данеш; Ахуджа, Капил (2015). «Переработка подпространств Крылова для приложений CFD и новый гибридный решатель рециклинга». Журнал вычислительной физики . 303 : 222. arXiv : 1501.03358 . Bibcode : 2015JCoPh.303..222A . DOI : 10.1016 / j.jcp.2015.09.040 .
- ^ Галлия, Андре (2014). Методы рециклирования подпространств Крылова для последовательностей линейных систем . ТУ Берлин. DOI : 10,14279 / depositonce-4147 .
- ^ Stoer и Bulirsch, §8.7.2
Заметки
- А. Майстер, Numerik linearer Gleichungssysteme , 2-е издание, Vieweg 2005, ISBN 978-3-528-13135-7 .
- Ю. Саад, Итерационные методы для разреженных линейных систем , 2-е издание, Общество промышленной и прикладной математики , 2003. ISBN 978-0-89871-534-7 .
- Ю. Саад и М. Х. Шульц, "GMRES: обобщенный алгоритм минимальной невязки для решения несимметричных линейных систем", SIAM J. Sci. Стат. Comput. , 7 : 856-869, 1986. DOI : 10,1137 / 0907058 .
- SC Eisenstat, HC Elman и MH Schultz, "Вариационные итерационные методы для несимметричных систем линейных уравнений", Журнал SIAM по численному анализу , 20 (2), 345–357, 1983.
- Дж. Стоер, Р. Булирш, Введение в численный анализ , 3-е издание, Спрингер, Нью-Йорк, 2002. ISBN 978-0-387-95452-3 .
- Ллойд Н. Трефетен и Дэвид Бау, III, Численная линейная алгебра , Общество промышленной и прикладной математики, 1997. ISBN 978-0-89871-361-9 .
- Донгарра и др. , Шаблоны для решения линейных систем: строительные блоки для итерационных методов , 2-е издание, SIAM, Филадельфия, 1994
- Амриткар, Амит; де Стерлер, Эрик; Свиридович, Катажина; Тафти, Данеш; Ахуджа, Капил (2015). «Переработка подпространств Крылова для приложений CFD и новый гибридный решатель рециклинга». Журнал вычислительной физики 303: 222. DOI: 10.1016 / j.jcp.2015.09.040