Экономичные критерии останова итераций в методе сопряженных градиентов Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 20, № 2, 2015

Экономичные критерии останова итераций в методе сопряженных градиентов

И. В. Киреев

Институт вычислительного моделирования СО РАН, Красноярск, Россия Контактный e-mail: kiv@icm.krasn.ru

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

Ключевые слова: метод сопряженных градиентов, критерии останова итераций.

Введение

Идеальный алгоритм решения систем линейных алгебраических уравнений теоретически должен иметь конечное завершение, но при досрочном завершении он должен дать приемлемое приближение к искомому решению. Одним из таких методов является метод сопряженных градиентов (MCG).

История MCG началась с работы [1]. Этот стандартный на сегодня вычислительный алгоритм используется для решения больших систем уравнений с симметричными положительно определенными матрицами, для которых прямые методы непрактичны. Алгоритм MCG можно применять не только для решения систем линейных уравнений, но и для исследования спектра матрицы [2].

Журнал "Computing in Science and Engineering" назвал метод сопряженных градиентов одним из 10 лучших алгоритмов XX века [3]. И тем приятнее отметить, что алгоритм ортогонализации степенной последовательности [4], лежащий в основе теоретического обоснования MCG, предложил русский математик А.Н. Крылов. История создания MCG и его развитие (до 1977 г.) отражены в обзоре [5].

1. Постановка задачи

При численном решении системы линейных уравнений каким-либо методом возникает естественный вопрос: как ошибки округления в компьютере влияют на полученное приближенное решение? Существуют два подхода к оценке точности алгоритма [6]: прямой, когда осуществляется непосредственный анализ погрешности (данные точные, но сам алгоритм работает с погрешностью), и обратный, основанный на идее, что реально вычисленное решение рассматривается как точное для той же задачи, но с возмущенными входными данными. При этом само возмущение выбирается так, чтобы его действие оказалось эквивалентным совокупному влиянию всех ошибок округления. Обратный

© ИВТ СО РАН, 2015

анализ позволяет оценить совместное влияние ошибок округления и ошибок входных данных на точность результатов.

В 1992 г. впервые был проведен обратный анализ погрешности метода сопряженных градиентов [7]. При этом показано [8], что характер погрешности численного решения линейной системы уравнений МСО при конечной точности машинной арифметики имеет сходство с точными вычислениями по этому же алгоритму, но приложенными к иной матрице, у которой гораздо больше, нежели у первоначальной, собственных значений, распределенных, однако, вблизи спектральных чисел исходной матрицы.

Метод сопряженных градиентов при точных вычислениях приводит к ответу за конечное число шагов, но по сути является итерационным процессом. Его слабым местом является критерий останова — определение номера шага процесса, после которого точность приближения к решению системы линейных уравнений на данном компьютере не может быть существенно улучшена. Практическое применение результатов обратного анализа погрешности для выработки условий прекращения вычислений в МСО опирается на априорные оценки границ спектра матрицы системы, получение которых является очень трудоемкой задачей. Поэтому весьма актуально построение экономичных критериев останова для ЫОО, что является целью данной работы.

2. Численная реализация МСС

Рассмотрим совместную систему линейных алгебраических уравнений

Ах = Ь, (1)

где х, Ь Е И™ — т-мерное евклидово пространство, А — квадратная симметричная положительно определенная вещественная матрица порядка т. При сделанных ограничениях решение х системы (1) доставляет минимум квадратичному функционалу

а(у) = 2 (АУ, У)-(Ь, У), (2)

т

где (Ь, у) = Ьк■Ук — скалярное произведение векторов Ь и у из Ит; евклидову норму к=1 _

вектора у обозначим через ||у || = л/(у, у). Для тех же векторов Ь и у их энергетическое

скалярное произведение обозначим через (Ь, у) а = (Ь, Ау). Тогда ||у||А = л/(у, у)А —

энергетическая норма вектора у. Необходимо заметить [9], что из (1) и (2) следует

равенство ||х — у ||А = 2а(у) + ||х||А, т. е. величина а(у) характеризует близость вектора у

к решению х системы (1) в энергетической норме.

Систему (1) численно можно решить методом сопряженных градиентов [9]:

при п = 0 задан х0; вычисляем г0 = Ах0 — Ь и полагаем б1 = г0;

при п > 0 вычисляем ап и вп по одной из формул

||г"—1||2 (гП-1 бП) (гП бП)^ ||гП ||2

аП = — /гп— 1 чп\ „ = |йП|2 < 0 вп = 1кп||2 = IIгп-1II2 > 0

(г , б М ||б ||б ||г Ц

и находим

гп = гп—1 + апАБп, бп+1 = гп + впБп, хп = хп—1 + апБп

Если все вычисления произведены при отсутствии ошибок округления, то соотношения (3), (4) гарантируют получение решения х системы (1) не более чем за т шагов. При этом совокупность векторов (б"} будет А-ортогональной ((б", Бк)а = 0, п = &), а векторы (г"} образуют ортогональную, т.е. сопряженную систему ((г", гк) = 0, п = &), причем г" является градиентом функции а(у) из (2) в точке х":

da(y)

dy

= rn = Axn - b. (5)

y=xn

Последнее обстоятельство отражено и в названии метода — Method of Conjugate Gradients, MCG.

Вариант MCG, в котором на каждом шаге процесса выполнены условия сопряженности (rn, rn-1) = (sn+1, s"")a = 0, связан с выражениями

lli-n-1 ||2 /ч-n ап\ ,

а =__НЕ_L_ Q =_ (r , s )A (6)

an = (rn-1, sn)A , вп = ||sn|2 . (6)

На каждом шаге версии (4), (6) МСО необходимо вычислить одно произведение матрицы А на вектор б" и четыре скалярных произведения: ||г"||2, (г"-1, б")а, (г", б")а „ ||0"|| 2 и ||б уА.

Напомним, что в авторской [1] версии МСО параметры ап и в" находятся по формулам

а

а = - ^_^ в = -(г , б )а .м

||ч"|| 2 , в" ||„"|| 2 . У)

||б 11А ||Б 11 А

Численная реализация алгоритма (4), (7) МСО требует меньших вычислительных затрат по сравнению с предыдущей версией — необходимо вычислить Аб", ||б"||А, (г"-1, б") и (г",б")а.

Предложенная в [10] и алгоритмически реализованная в [2] версия МСО, в которой ап и в" определяются соотношениями

|г"-1 у 2 |г"12

а" = -(г"-1 8")а , в" = у г"-1 у 2 , (8)

считается наиболее экономичной, поскольку на очередном шаге алгоритма (4), (8) необходимо вычислять одно произведение матрицы на вектор и только два скалярных произведения.

В монографии [11] обсуждается версия МСО, в которой гарантируется релаксаци-онность процесса, т.е. выполнение неравенства а(х") < а(х"-1) на каждом шаге п, для чего предложено делать два вычисления произведения матрицы на вектор. Однако, если рассматривать МСО как метод генерации направления спуска [12] для численного решения задачи минимизации квадратичной функции а (у), то релаксационность можно получить и за одно вычисление произведения матрицы на вектор. Покажем, как этого можно достичь.

Зафиксируем векторы у, б и обозначим через а5 решение одномерной задачи безусловной минимизации [12] квадратичной сильно выпуклой функции а(у + аэ) по вещественному аргументу а:

= _ (Ау - Ь, б) = (х - у, б)а

II II2 11112 ,

||б||А ||б||А

a(y + ass) = a(y) - ^(as ■ ||s||a)2 = mina(y + as).

2 «GR

Очевидно, что если направление спуска s не ортогонально градиенту функции a(y) в "точке" y, то будет выполнено неравенство a(y — ass) < a(y), означающее релаксаци-онность перехода от "точки" y к "точке" y — ass посредством "спуска по направлению s". Поэтому, если an и вп найти по формулам

= _ (xn-1, sra)a — (b, sn) = _ (rn, sra)a (Q)

an = ||sn|| 2 , вп = ||sn|| 2 , (Q)

||s Ha ||s Ha

эквивалентным (3) при точных вычислениях, то релаксационность процесса (4), (Q) гарантирована с той точностью, с которой вычисляются скалярные произведения и произведение матрицы на вектор. В такой версии MCG на каждом шаге следует выполнить одно произведение матрицы A на вектор sn и четыре скалярных произведения: (xn-1, sn)A, (bn-1, sn), (rn, sn)A и ||snHA.

Необходимо отметить, что при минимизации выпуклых гладких функций [12] методом сопряженных градиентов возникают соотношения, аналогичные (Q), в которых роль матрицы A играет матрица Гессе целевой функции, а выражение для в«. из (Q) обеспечивает A-ортогональность векторов sn и sn+1 с той же точностью, с которой вычисляется параметр an.

О реальной эффективности приведенных выше версий MCG можно судить по значению энергетической нормы невязки ||x — xn||a на каждом шаге процесса. Для того чтобы отслеживать накопление ошибки, присущее только версии алгоритма MCG, была проведена серия численных расчетов для диагональных матриц A; подобный выбор матрицы системы (1) гарантирует, что "кососимметричная составляющая" погрешности вычисления произведения матрицы A на вектор sn будет пренебрежимо мала. На рис. 1 в графическом виде приведены результаты вычислительных экспериментов для систем

Рис. 1. Изменение по шагам МСО энергетической нормы невязки х — хп; (6)—(9) — версии МОО (то же на рис. 2)

с матрицами вида Лц = i-M(i-1), i = 1, 2,... ,m; здесь M — натуральное число. Число обусловленности такой матрицы condA = A11/Amm = (m — 1)M(m-1) очень быстро увеличивается с ростом m и M: если m = 256, то при M =2 cond A ~ 6.6 • 104, а при M = 2 cond A ~ 4.3 • 109. Вычисления проводились в режиме "double", при котором относительная погрешность представления чисел в компьютере не превышает величину £с = 2-52 ^ 2.22045 • 10-16, а наименьшее, отличное от нуля, положительное число равно 2-1022 « 2.22507 • 10 - 308 (машинный ноль). Эти характеристики вычислительного процесса можно легко вычислить [13].

Отметим монотонность приведенных зависимостей, характер которых в основном определяется числом обусловленности матрицы системы и величиной £с. Это обстоятельство говорит о том, что целесообразно не прерывать итерации после m шагов [2]. Кроме того, "экспериментально" установлено, что наименее эффективной оказывается релаксационная версия MCG (4), (9), в то время как алгоритм (4), (6), обеспечивающий на каждом шаге MCG условия сопряженности (rn, rn-1) = (sn+1, s"")a = 0, оказался самым точным из рассматриваемых.

3. Критерии останова итераций MCG

Весьма сложным вопросом численной реализации MCG является определение рационального числа итераций, т. е. момента, начиная с которого полученное приближенное решение xn не может быть существенно улучшено и итерационный процесс целесообразно прекратить.

Наиболее естественным было бы остановить процесс в тот момент, когда величина ||x — хп||а становится соизмеримой с погрешностью ее вычисления. Линейная теория накопления ошибок вычислений [7, 13, 14] позволяет, отправляясь от точностной характеристики вычислений £с, получить верхнюю оценку |i||x — хп||а| погрешности вычисления ||х — хп||а. Рисунок 2 иллюстрирует типовой характер зависимости |$||x — хп||а| от относительного номера итерации n/m.

Рис. 2. Верхняя оценка абсолютной погрешности вычисления энергетической нормы невязки

Анализ численных экспериментов и "затратность" вычисления как энергетической нормы, так и верхней (как правило, весьма завышенной) оценки погрешности ее вычисления приводят к выводу, что использовать критерии останова итерационного процесса MCG на базе линейной теории накопления ошибок нецелесообразно.

В работе [10] предложено оценивать накопленную погрешность вычислений, анализируя векторы rn и Axn—b. Численные эксперименты показали, что величины log10 ||rn|| и log10 ||Axn — b — rn|| для MCG-версий (4), (6)-(4), (8) меняются с ростом n схожим образом (рис. 3, а), тогда как соотношения (4), (Q) дают несколько иную картину (рис. 3, б).

0 -5 -10 -15

2 4

m = 256, M = 2

n/m

-5

-10

-15

log10 ||Axn — b — rn|| — ~

_I_I_I_I_I_I_I_I_l_

m = 256, M = 4

n/m

-10

-15-

2 3

m = 256, M = 2

n/m

-5

-10

-15

log10 ||Axn — b — rn|| ыгут^

L-l_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_l_

12 3 4

m = 256, M = 4

n/m

Рис. 3. Изменение ||rn|| и ||Axn — b — rn|| по шагам MCG на основе соотношений (4), (7) (а) и (4), (9) (б)

а

0

2

4

0

0

1

4

0

Общим для рассматриваемых версий МСО является то, что для систем, числа обусловленности которых не очень велики, евклидова норма ||Ах" — Ь — г"|| с ростом числа итераций меняется незначительно, тогда как величина ||г"|| уменьшается начиная с некоторого шага. Поэтому в [2] в качестве критерия останова предлагается проверять неравенство

||г"|| < е("/т)2||Ах" — Ь — г"||, (10)

где п — номер текущей итерации, т — порядок исходной системы уравнений (1); при нарушении (10) итерации прекращаются. Экспоненциальный множитель е("/т) призван, наряду с ограничением числа шагов сверху [2] величиной 5т, бороться с плохой обусловленностью системы уравнений.

Каждая из приведенных выше версий МСО по своей природе есть разновидность метода спуска решения задачи минимизации квадратичного функционала (2), а характерной особенностью таких алгоритмов [16] является высокая производительность первых итераций. Поэтому можно отследить выполнение неравенства (10), рассматривая лишь в"+1-составляющие соответствующих векторов.

Если ввести е" = в"/||в"||, то

(е"+1, Ах"— Ь — г") = ((Ав"+1, х") — (в"+1, Ь) — (в"+1, гга)) /||в"+1'

и для определения момента окончания итерационного процесса вместо (10) проверяем неравенство

|(в"+1, г")| < е("/т)2 |(Ав"+1, х") — (в"+1, Ь) — (в"+1, г")|. (11)

Так как величины (в", г"-1) и Ав" определяются на каждом шаге МСО при построении коэффициентов ага и в", то на очередной итерации можно не вычислять произведение матрицы А на вектор х" [17], а достаточно найти лишь скалярные произведения (Ав"+1, х") и (в"+1, Ь), что значительно уменьшает вычислительные затраты.

Совместное изменение величин |(е"+1, Ь)| и |(е"+1, Ах" — Ь —г")| отражено на рис. 4, а для версии (4), (7) и на рис. 4, б для версии (4), (9). Сравнение последних графиков с рис. 3 показывает хорошую коррелированность последовательностей ||(е"+1, г")|} и {||г"||}, {|(е"+1, Ах" — Ь — г")|} и {||Ах" — Ь — г"||} для рассматриваемых версий МСО.

Анализ вычислительных экспериментов [18] показывает, что рассмотренные критерии останова наименее эффективно работают в релаксационной (4), (9) версии МСО, для которой, однако, возможен и иной критерий остановки вычислений.

Рассмотрим последовательность {а"} числовых значений

= 2а(х"), (12)

члены которой, согласно (9), можно вычислить по рекуррентной формуле

а" = а"-1 — в|")"А||— ^ 2 , ао = (г0, х0) — (Ь, х0). (13)

С другой стороны, исходя из определения (2) квадратичного функционала а(у), чис-

и х":

а" = (г", х") —(Ь, х"). (14)

ло а" можно выразить через г" и х"

en+1 r«,)|

2

m = 256, M = 2

n/m

-5

-10

-15

-20

log10 |(en+1, rn)|

log10 |(en+1Axn— b — rn)|

m = 256, M = 4

n/m

log10 |(en+1, rn)|

-10

-15 t

-20

m = 256, M = 2

n/m

-5

-10

-15

-20

log10 |(en+1, rn)

log10 |(en+1Axn— b — rn)|

2

m = 256, M = 4

n/m

Рис. 4. Изменение модулей проекций векторов rn и Axn — b — rn на направление спуска sn+1 по шагам MCG на основе соотношений (4), (7) (а) и (4), (9) (б)

а

0

4

0

4

Замечание к формуле (2) наводит на мысль, что разность между вычисленными по формулам (13) и (14) значениями ага может служить оценкой накопленной в процессе вычислений погрешности определения хп, и требование, чтобы на каждой итерации величина |ага — ага-1| не превышала этой погрешности, является вполне естественным. Введем

= (Ъ, хп), с„ = (Л хп). (15)

Тогда из (4) следует

= (Ъ, хп) = (Ъ, хп-1 + а„вп) = 6га_1 + ага(Ъ, вга), Ьо = (Ъ, х0).

Аналогично для чисел С" получаем

(г", х") = (г", х"-1 + «"в") = (г"-1 + «"Ав", х"-1) + а"(г", в")

С"-1 + а"(х"-1, в")а + а"(г", в"), Со = (г0, х0).

Теперь можем сформулировать для МСО в форме (4), (9) критерий прекращения вычислений: считая х0 заданным,

при п = 0, вычисляем г0 = Ах0 — Ь и полагаем в1 = г0;

а0 = (г0, х0) —(Ь, х0), Ь0 = (Ь, х0), С0 = (г0, х0);

при п > 0 вычисляем (х"-1, в")а, (Ь, в"), 11в"|А и (г", в")а и определяем

а" = а" 1

(в", х"-1)А — (в", Ь) ||в" ||а

Ь" = 6Т(-1 + а"(Ь, в"), С" = С"-1 + а" ((х"-1, в")а + (г", в"))

1С" Ь" а"1;

тогда, если на шаге п ИОС (4), (9) окажется, что

(в", х"-1)А — (в", Ь)

А_ Й|а

< е("/т)2£",

16)

то вычисления прекращаем.

Характер совместного изменения величин |2а(х") — 2а(х"-1 )| и |2а(х") — а" отражен на рис. 5. Из сравнения этих данных с рис. 2 можно сделать вывод, что точка пересечения графиков величин |2а(х") — 2а(х"-1)| и |2а(х") — а" соответствует тому значению п, при котором начинает устойчиво возрастать верхняя оценка |$||х — х"||а| погрешности вычислений энергетической нормы невязки ||х — х"||а.

2

3

т = 256, М = 2

п/т

0

-5

-10

1о§10 |2а(х") — аГ1

_|_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_I_1_

2

3

4

т = 256, М = 4

п/т

Рис. 5. Пошаговое для МСО на основе соотношений (4), (9) поведение величин |2а(х") — а" и |2а(х") — 2а(х"-1)| для различных чисел обусловленности

С

"

0

1

о

4

Заключение

В работе для четырех наиболее распространенных программных реализаций метода сопряженных градиентов решения систем линейных алгебраических уравнений численно исследовано пошаговое поведение энергетических норм невязок с точным решением при наличии погрешностей округления. Численные эксперименты показали, что все рассмотренные версии метода приближают искомое решение примерно с одинаковой погрешностью при числе итераций, не большем числа неизвестных. Однако картина существенно меняется для итераций с номерами, превышающими размерность системы в 3-5 раз. В этом диапазоне шагов наиболее точное приближение дает версия метода, в которой на каждом шаге обеспечиваются ортогональность в энергетической метрике направлений спуска и евклидова перпендикулярность векторов невязок, вычисляемых во время процесса. Наиболее медленным оказался алгоритм, в котором параметры итерационного процесса выбирались из условия его релаксационности. Отношение энергетических норм невязок в этих двух случаях достигало трех десятичных порядков.

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

Благодарности. Работа выполнена при финансовой поддержке РФФИ (грант № 1401-00130).

Список литературы / References

[1] Hestenes, M.R. and Stiefel, E. Methods of conjugate gradients for solving linear systems // J. of Res. of National Bureau of Standards. 1952. Vol. 49, No. 5. P. 409-435.

[2] Wilkinson, J.H. and Reinsch, С. Handbook for Automatic Computation. Vol. II. Linear Algebra. N.Y.: Springer-Verlag, 1971. 450 p.

[3] Cipra, B.A. The best of the 20th century: Editors name top 10 algorithms // Comput. Sci. Eng. 2000. Vol. 332, No. 4. P. 291-293.

[4] Крылов А.Н. О численном решении уравнения, которым в технических вопросах определяются частоты малых колебаний материальных систем // Изв. АН СССР. Отделение математических и естественных наук. 1931. № 4. С. 491-539.

Krylov, A.N. On the numerical solution of the equation by which in technical questions frequencies of small oscillations of material systems are determined // Izv. AN SSSR (News of Academy of Sciences of the USSR). Otdel. Matem. i Estest. Nauk. 1931. No. 4. P. 491-539. (in Russ.)

[5] Golub, G.H. and O'Leary, D.P. Some history of the conjugate gradient and Lanczos algorithms: 1948-1976 // SIAM Rev. 1989. Vol. 31, No. 1. P. 50-102.

[6] Воеводин В.В. Вычислительные основы линейной алгебры. М.: Наука, 1977. 304 c. Voevodin, V.V. Numerical Foundations of Linear Algebra. Moscow: Nauka, 1977. 304 p. (in Russ.)

[7] Greenbaum, A. Estimating the attainable accuracy of recursively computed residual methods. Techn. Rep. TR-95-1515. Department of Comput. Sci., Cornell Univ., Ithaca. N.Y. (May 1995).

[8] Greenbaum, A., Strakos, Z. Predicting the Behavior of Finite Precision Lanczos and Conjugate Gradient Computations. Techn. Rep. TR-538, Department of Comput. Sci., Cornell Univ., Ithaca. N.Y. (January 1991).

[9] Воеводин В.В. Линейная алгебра. М.: Наука, 1980. 400 с. Voevodin, V.V. Linear Algebra. Moscow: Nauka, 1980. 400 p. (in Russ.)

[10] Engeli, M., Ginsburg, Th., Rutishauser, H., Stiefel, E. Refined Iterative Methods for Computation of the Solution and the Eigenvalues of Self-Adjoint Boundary Value Problems. Mitteilungen aus dem Institut fur angewandte Mathematik der ETH Zurich, No. 8. Basel, Stuttgart: Birkhaiiser Verlag, 1959. 107 p.

[11] Воеводин В.В. Численные методы алгебры. Теория и алгорифмы. М.: Наука, 1966. 248 с. Voevodin, V.V. Computational Methods of Algebra. Theory and Algorithms. Moscow: Nauka, 1966. 248 p. (in Russ.)

[12] Карманов В.Г. Математическое программирование: Учеб. пособие. М.: Физматлит, 2004. 264 с.

Karmanov, V.G. Mathematical Programming: Textbook. Moscow: Fizmatlit, 2004. 264 p. (in Russ.)

[13] Малышев А.Н. Введение в вычислительную линейную алгебру. Новосибирск: Наука, 1991. 229 с.

Malyshev, A.N. Introduction to Computational Linear Algebra. Novosibirsk: Nauka, 1991. 229 p. (in Russ.)

[14] Годунов С.К. Решение систем линейных уравнений. Новосибирск: Наука, 1980. 178 с. Godunov, S.K. Solving Systems of Linear Equations. Novosibirsk: Nauka, 1980. 178 p. (in Russ.)

[15] Meurant, G. and Strakos, Z. The Lanczos and conjugate gradient algorithms in finite precision arithmetic // Acta Numerica. 2006. Vol. 15. P. 471-542.

[16] Фаддеев Д.К., Фаддеева В.Н. Вычислительные методы линейной алгебры. М.: Физ-матгиз, 1963. 229 с.

Faddeev, D.K., Faddeeva, V.N. Computational Methods in Linear Algebra. Moscow: Fizmatgiz, 1963. 229 p. (in Russ.)

[17] Чернышёва А.А., Киреев И.В. Модификация критерия Уилкинсона остановки итераций в методе сопряженных градиентов // Вест. КрасГУ. 2005. № 4. C. 173-177. Chernysheva, A.A., Kireev, I.V. Modification Wilkinson's criteria stopping iterations in the conjugate gradient method // Bulletin of Krasnoyarsk State Univ. 2005. No. 4. P. 173-177. (in Russ.)

[18] Киреев И.В. Численная реализация метода сопряженных градиентов. Красноярск, 2013. (Препр. ИВМ СО РАН; № 13-1. 26 с.)

Kireev, I.V. The Computational Implementation of the Conjugate Gradient Method. Krasnoyarsk, 2013. (Preprint. ICM SB RAS; № 13-1. 26 p.) (in Russ.)

Поступила в 'редакцию 1 октября 2014 г., с доработки — 9 февраля 2015 г.

Inexpensive stopping criteria in the conjugate gradient method

Kireev, Igor' V.

Institute of Computational Modelling SB RAS, Krasnoyarsk, 660036, Russia Corresponding author: Kireev, Igor' V., e-mail: kiv@icm.krasn.ru

In the paper, some aspects of the numerical implementation of the conjugate gradient method (CGM) for systems of linear algebraic equations with symmetric positive definite matrix in the presence of round-off errors are discussed. With exact calculations, CGM provides an exact solution in a finite number of iteration steps. But in fact CGM is an iterative process and the weak point in an iterative process is in a stopping criterion. It is required to determine the number of the iteration step, after which the accuracy of an approximation to a solution of a system of linear equations may not be considerably improved with a particular computer. Hence, the construction of inexpensive stopping criteria for CGM being the aim of this paper is an urgent problem. For four popular versions of CGM, the step-by-step behavior as well as stopping criteria for an iterative process are considered. Numerical results show that the most accurate approximation is achieved by the CGM-version where descent directions and residual vectors are orthogonal in the energy and Euclidean metrics, respectively, at each iteration step. A practical stopping criteria for CGM is proposed as a formula that enables one to determine the number of the CGM iteration step, starting with which the progress is no longer being made. The application of the constructed criteria to the solution of specific systems of linear algebraic equations with ill-conditioned matrices is demonstrated. Keywords: conjugate gradient method, stopping criteria.

Acknowledgements. This work has been supported by the Russian Foundation for Basic Research (RFBR) grant No. 14-01-00130.

Received 1 October 2014 Received in revised form 9 February 2015

© ICT SB RAS, 2015