Общая схема двухслойного итерационного метода. Сходимость.

Среди недостатков прямых методов решения систем ЛАУ следует отметить сравнительно большие вычислительные затраты, возрастающие в общем случае пропорционально , и склонность к накоплению вычислительной погрешности. В силу этого прямые методы используются, как правило, при решении систем сравнительно небольшой размерности .

В большинстве случаев при решении систем ЛАУ большой размерности приходится иметь дело с так называемыми разреженными матрицами, количество ненулевых элементов которых возрастает пропорционально при общем числе элементов . Нулевые элементы матрицы не несут полезной информации и могут быть проигнорированы в матричных вычислениях. В рамках прямых методов игнорирование нулевых элементов матрицы весьма нетривиальная задача. Основная проблема здесь в том, что в процессе матричных преобразований общее число ненулевых элементов не сохраняется – возрастает. Например, суммарное число ненулевых элементов матриц при разложении или при факторизации Холецкого может многократно превосходить число ненулевых элементов исходной матрицы. Предварительное упорядочение элементов матрицы при использовании алгоритмов факторизации позволяет в ряде случаев снизить остроту проблемы, однако полностью устранить рост числа ненулевых элементов в общем случае, по-видимому, не представляется возможным.

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

Общая схема большинства итерационных методов решения систем ЛАУ

(3.1)

с невырожденной матрицей и заданным вектором правой части имеет вид

, , (3.2)

где – матрица итерационного метода или матрицей перехода от k-ой итерации к (k+1)-ой, – произвольный вектор, называемый начальным приближением итерационного процесса. Последовательность , будем называть итерационными приближениями искомого решения.

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

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

1. последовательность векторов , , сходится.

2. предел данной последовательности является решением (3.1).

Из условия 2 следует, что матрица и вектор могут быть заданы в виде:

, , (3.3)

где – единичная матрица, – произвольная невырожденная матрица, выбор которой призван удовлетворить условие 1. Заметим, что для произвольных невырожденных матриц и существует единственное значение вектора такое, что и с учетом выбора (3.3)

.

Теорема 3.1. Пусть , – решение (3.1). Тогда итерационный процесс (3.2), (3.3) сходится со скоростью геометрической прогрессии к вектору и для погрешности итерационного метода выполняется оценка

. (3.4)

. Доказательство: Вычитая из (3.2) равенство для погрешности метода имеем:

.

Поскольку , то . Теорема доказана.

Разнообразие итерационных методов связано с выбором конкретного вида матрицы , которую принято называть переобусловливателем (preconditioner). Если матрица одинакова для всех итераций, то такой итерационный процесс называется стационарным. Среди нестационарных итерационных методов традиционно используются переобусловливатели вида , где итерационный параметр для каждой итерации выбирается из расчета наибольшей скорости сходимости. С точки зрения алгоритмической реализации итерационный процесс (3.2), (3.3) удобно представить в виде

. (3.5)

При использовании (3.5), в отличие от (3.2), (3.3), не требуется явный вид матрицы , и вычисление очередного итерационного приближения сводится к решению системы алгебраических уравнений вида

. (3.5')

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

• Система уравнений должна допускать существенно более эффективный способ решения, нежели исходная система (3.1).
• Матрица итерационного метода по норме должна иметь как можно меньшее значение .

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

. 3 .2 Вычислительная сложность итерационных методов. Число итераций.

Вычислительные затраты при решении систем ЛАУ с использованием итерационных методов определяются двумя основными факторами. Во-первых, это число итераций, требуемых для достижения заданной точности, а во-вторых – вычислительные затраты на реализацию одной итерации. Именно эти две характеристики позволяет определять и сравнивать эффективность различных итерационных методов. В ряде случаев удается получить вполне реалистичные теоретические оценки скорости сходимости.

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

.

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

, (3.6)

где – некоторая малая величина, определяющая желаемую погрешность искомого приближенного решения. Количество итераций для достижения заданной точности можно оценить, зная норму матрицы итерационного процесса. Если положить, что при нулевом начальном приближении относительная погрешность равна 100%, тогда, например, для достижения точности с учетом оценки (3.4) потребуется число итераций такое, что

. (3.7)

Норма матрицы итерационного процесса характеризует скорость его сходимости. Максимальная скорость сходимости достигается при минимальном значении . Оценка (3.7) показывает, что минимальное число итераций для достижения заданной точности может неограниченно возрастать при .

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

Положим для начала, что . В этом случае одна итерация согласно (3.5) включает одно умножение матрицы на вектор, , и сложение трех векторов, . Если , то каждая итерация (3.5) требует порядка арифметических операций. При количестве итераций вычислительная сложность итерационного процесса оценивается числом арифметических операций порядка , т.е. сопоставима с вычислительной сложностью метода Гаусса.

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

С другой стороны, целесообразность использования переобусловливателя связана с возможностью уменьшения числа итераций для достижения заданной точности решения (переобусловливатель – единственный свободный параметр итерационного процесса (3.5), выбор которого способен влиять на норму матрицы итерационного процесса). Вычислительная сложность отдельной итерации при этом возрастает ровно настолько, сколько требуется для решения системы уравнений вида . Таким образом, использование переобусловливателя способно реально снизить вычислительные затраты на поиск приближенного решения, если обеспеченное им относительное уменьшение числа итераций превосходит относительный рост вычислительной сложности отдельной итерации.

3 .3 Метод простой итерации.

Стационарный одношаговый итерационный метод вида

(3.8)

называется методом простой итерации. В трактовке общей схемы одношаговых методов (3.2) матрица метода простой итерации имеет вид . В данном методе фактически отсутствует переобусловливатель. Его роль выполняет единичная матрица и итерационный параметр . Согласно представлению (3.5), для метода простой итерации .

Теорема 3.2. Пусть – симметричная положительно определенная матрица , тогда итерационный метод (3.8) сходится при .

Доказательство: Спектральная норма симметричной матрицы определяется модулем ее наибольшего собственного значения: . Если – симметричная матрица, то матрица также будет симметричной и . Тогда . Из положительной определенности матрицы следует, что при выполняется оценка , из которой следует, что .

Определим условия, при которых скорость сходимости итерационного метода (3.8) максимальна.

Теорема 3.3. Пусть – симметричная положительно определенная матрица: , , где положительные постоянные , – соответственно минимальное и максимальное собственные значения матрицы . Тогда максимальная скорость сходимости итерационного процесса (3.8) достигается при , при этом

. (3.9)

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

.

Несложно заметить, что и , следовательно, в интервале значений функция принимает минимальное значение. Поскольку функция определяется максимальным значением модулей двух линейных функций, то минимум такой функции может достигаться только в точке равенства модулей данных линейных функций. Уравнение имеет единственный корень на интервале : . При этом

.

Теорема доказана.

Заметим, что скорость сходимости метода простой итерации зависит от отношения даже в случае оптимального выбора итерационного параметра. Для симметричных положительно определенных матриц , . Следовательно , где – число обусловленности матрицы . В случае плохо обусловленных матриц значение велико, и тогда согласно оценке (3.9)

.

Таким образом, эффективность метода простой итерации может катастрофически ухудшаться при .

3 .4 Оптимальный выбор параметра нестационарного итерационного метода.

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

Для решения системы ЛАУ с положительно определенной матрицей используем итерационный метод

. (3.10)

В отличие от метода простой итерации (3.8) в итерационном процессе (3.10) используется переменный итерационный параметр. Определим, каким должно быть значение итерационного параметра , чтобы норма погрешности для очередного итерационного приближения имела минимальное значение.

Для погрешности итерационного метода (3.10) имеем

. (3.11)

Умножим скалярно уравнение (3.11) само на себя, предварительно умножив его слева на матрицу :

.

Условие минимума нормы погрешности определим из равенства нулю ее производной:

.

Из последнего равенства имеем

.

Учитывая, что , получаем выражение для оптимального значения итерационного параметра

, . (3.12)

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

. (3.13)

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

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

,

то скорость сходимости итерационного метода (3.10) будет определяться показателем геометрической прогрессии

. (3.14)

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

3 .5 Градиентные методы. Методы наискорейшего спуска.

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

(3.15)

эквивалентна задаче минимизации функционала (квадратичной формы)

(3.16)

Использование данной эквивалентности лежит в основе так называемых градиентных итерационных методов для систем ЛАУ с положительно определенной симметричной матрицей. Наиболее распространенные из градиентных методов являются методы наискорейшего спуска и метод сопряженных градиентов.

Несложно вычислить градиент квадратичной формы и убедиться, что для любого вектора он совпадает с невязкой системы:

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

(3.17)

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

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

Для вычисления величины шага рассмотрим условие минимума функционала в направлении вектора невязки

Таким образом, величина шага должна быть такой, чтобы новая невязка была ортогональна невязке на предыдущем итерационном приближении. Для вычисления умножим слева уравнение (3.17) на A и вычтем из полученного равенства .

Последнее равенство умножим скалярно на , в результате чего, учитывая требование ортогональности и , имеем

. (3.18)

Метод наискорейшего спуска (3.17), (3.18) и по структуре и по скорости сходимости аналогичен методу минимальных невязок (3.13).

3 .6 Методы сопряженных градиентов.

Выбор направления для минимизации функционала в методе наискорейшего спуска представляется вполне разумным, но вряд ли самым оптимальным. Наблюдая за траекторией сходимости решения системы двух уравнений (см. Рис. 3.1.) можно заметить, что каждое новое приближение (за исключением некоторых частных случаев выбора начального приближения) находится в направлении, которое не совпадает с направлением на точку глобального минимума функционала. Более того, в методе наискорейшего спуска одно и то же направление поиска может использоваться на нескольких итерациях.

Среди градиентных итерационных методов наибольшее распространение получил метод сопряженных градиентов (Conjugate Gradients). Его отличительная особенность состоит в том, что теоретически он позволяет получить точное решение за фиксированное число итераций, не превосходящее размерности матрицы (т.е. он аналогичен прямым методам решения СЛАУ).

Структура метода сопряженных градиентов (CG) во многом напоминает метод наискорейшего спуска. В частности, первая итерация CG метода выполняется по тому же алгоритму, что и для метода скорейшего спуска

(3.19)

(3.20)

(3.21)

(3.22)

Далее, стратегия выбора направления для минимизации квадратичной формы меняется, и последующие итерации вычисляются по следующему алгоритму (вначале приведем структуру алгоритма, а затем рассмотрим идеи, положенные в основу данного метода):

(3.23)

(3.24)

; (3.25)

; (3.26)

(3.27)

(3.28)

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

Как работает метод сопряженных градиентов?

Заметим, что направление поиска глобального минимума квадратичной формы (решения задачи) определяется вектором d, который конструируется по рекуррентным формулам с учетом невязок задачи на всех предыдущих итерациях. Идея, приводящая к методу сопряженных градиентов, состоит в нахождении такого набора линейно независимых направлений , чтобы в каждом из этих направлений делать лишь один шаг. Это возможно лишь в том случае, если на каждом итерационном приближении происходит полное исключение погрешности вдоль текущего направления поиска .

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

. (3.29)

Пусть , , т.е. направления поиска А -ортогональны, где А – симметричная положительно определенная матрица рассматриваемой системы ЛАУ. Если потребовать, чтобы на текущей итерации полностью устранялась k -я компонента погрешности, то каждое новое приближение должно находиться по следующему алгоритму

,

и тогда

. (3.30)

Умножая уравнение (3.30) слева на матрицу A, а затем, умножая полученное равенство скалярно на , с учетом имеем:

,

. (3.31)

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

С учетом А -ортогональности направлений поиска, для этого достаточно взять любой набор линейно независимых векторов и применить к нему алгоритм ортогонализации. В методе сопряженных градиентов для построение требуемой системы А-ортогональных векторов используется тот факт, что погрешность на -й итерации A-ортогональна вектору поиска (поскольку -я итерация исключает погрешность вдоль вектора поиска , а остальные вектора в разложении (3.30) А-ортогональны направлению ):

.

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

В качестве первого вектора используем вектор невязки нулевого приближения

.

Для получения следующего направления поиска следует взять невязку и вычесть из неё составляющие, которые не А-ортогональны вектору

, , .

Повторяя данную процедуру, приходим к формулам построения A-ортогональной системы векторов .

,

,

Выражение для вычисления коэффициентов можно упростить, используя тот факт, что

,

Из последнего равенства следует , откуда

С учетом того, что , , ,

Таким образом, мы получили формулы градиентного метода (3.20)-(3.28), использующего в качестве направлений поиска минимума функционала полную системы A-ортогональных векторов. При этом на каждом шаге итерационного метода невязка приближенного решения ортогональна направлению поиска, т.е. направления поиска являются сопряженными с вектором градиента в точке текущего итерационного приближения. Отсюда происходит название данного метода – метод сопряженных градиентов.

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

3 .7 Неявные итерационные методы и выбор переобусловливателя.

В случае плохо обусловленных матриц сходимость итерационных методов вида (3.2), (3, 3) с оператором может оказаться очень медленной (см. оценки (3.9), (3.14)). В данной ситуации даже выбор оптимальных итерационных параметров не способен решить проблему достижения желаемой скорости сходимости. Решение данной проблемы в большинстве случаев может быть получено с использованием неявных итерационных методов или итерационных методов с переобусловливателем.

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

, (3.15)

эквивалентен явному итерационному методу

, (3.16)

где . Основное функциональное назначение матрицы в том, чтобы в итерационных процессах (3.15) и (3.16) достичь существенного уменьшения числа обусловленности матрицы по сравнению с числом обусловленности исходной матрицы .

Второй важный момент при выборе переобусловливателя состоит в вычислительной эффективности – возможности вычисления матрицы (или решения системы ) намного эффективнее, нежели обращение матрицы (или решение системы ).

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

. Рассмотрим наиболее простые неявные итерационные методы.

Метод Якоби. . В качестве переобусловливателя используется диагональная матрица, элементы которой совпадают с диагональными элементами матрицы . Выбор диагонального переобусловливателя практически не увеличивает вычислительную сложность отдельной итерации по сравнению с явным методом. Данный метод может оказаться полезным для разреженных матриц с диагональным преобладанием в случае, когда диагональные элементы матрицы существенно отличаются друг от друга. Такие матрицы возникают, например, при дискретизации многомерных уравнений математической физики с сильно неоднородными коэффициентами. Если диагональные элементы матрицы одинаковы (или почти совпадают), то метод Якоби не имеет преимуществ по сравнению с явным методом.

Метод Зейделя (Гаусса-Зейделя). Матрица имеет треугольный вид и строится непосредственно из соответствующих элементов матрицы . В силу треугольности матрицы данный метод имеет небольшой рост вычислительных затрат на одну итерацию и примерно вдвое (в отдельных случаях и более чем вдвое) сокращает число итераций для достижения заданной точности по сравнению с явным методом. Использование данного метода практически всегда имеет положительный эффект по сравнению с явным методом и методом Якоби.

Метод последовательной верхней релаксации. В некотором роде является обобщением метода Зейделя и метода Якоби. Переобусловливатель строится из верхней треугольной части матрицы без главной диагонали и диагональных элементов матрицы

:

, , .

При метод последовательной верхней релаксации совпадает с методом Гаусса-Зейделя, при – метод совпадает с методом Якоби. Оптимальное значение параметра обычно лежит в интервале . В большинстве случаев метод последовательной верхней релаксации превосходит по эффективности методы Якоби и Зейделя. Метод весьма популярен для многомерных задач математической физики. Существуют модификации данного метода, основанные на чередовании верхней и нижней треугольных матриц в переобусловливателе .

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