Численные методы линейной алгебры

Белорусский государственный университет

 

В.М.Волков

 

 

Численные методы линейной алгебры

 

Теоретический минимум для студентов 5-го курса

Механико-математического факультета

(заочное отделение)

 

Минск 2009

В.М.Волков. Численные методы линейной алгебры.

Теоретический минимум для студентов 5-го курса механико-математического факультета

 

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

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

 

Введение

 

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

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

 

Задачи линейной алгебры.

Основные понятия

 

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

– сложение;

– умножение на скаляр;

– скалярное умножение[*].

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

 

Матрицу можно рассматривать также как упорядоченную совокупность векторов, образующих ее строки или столбцы. Индексы элемента определяют его позицию в k-той строке и m-том столбце (изменениям первого индекса соответствует изменение положения в столбце, а второго – положения в строке). Размерность матрицы указывает на то, что матрица имеет строк и столбцов соответственно. Если , то матрица называется квадратной и определяет оператор в N-мерном векторном пространстве.

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

1. диагональные матрицы: , ;

2. нижние (левые) треугольные матрицы: , ;

3. верхние (правые) треугольные матрицы: , .;

4. ленточные матрицы: , ;

5. матрицы перестановок, вращений и отображений.

 

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

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

, , .

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

.

 

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

,

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

 

Нормы векторов и матриц.

1. – положительная определенность; 2. , – линейность при умножении на скаляр; 3. – неравенство треугольника.

Системы линейных алгебраических уравнений. Разрешимость и устойчивость.

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

Прямые методы решения систем ЛАУ

 

Метод Гаусса

Большинство прямых методов решения систем ЛАУ в той или иной мере наследуют идею алгоритма последовательного исключения неизвестных – метода Гаусса.… , . (2.1) Последнее уравнение системы содержит всего лишь одно неизвестное ,.и его значение вычисляется первым:

Метод Гаусса с выбором ведущего элемента

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

LU-факторизация.

Основной теоретический результат, касающийся существования и единственности представления матрицы вида заключается в следующем Теорема. Если , то существует матрица перестановок такая, что имеет место… , (2.8)

Разложение Холецкого (метод квадратного корня).

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

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

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

Геометрический смысл и примеры плохой обусловленности матриц.

 

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

(2.26)

Каждое их уравнений системы (2.26) можно трактовать как уравнение прямой на плоскости, и тогда решение рассмотренной системы определяется координатами точки пересечения данных прямых.

Тангенсы углов, образуемых прямыми с осью , соответственно равны

, .

Вычислим определитель матрицы системы (2.26) и покажем, что условие эквивалентно условию :

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

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

. (2.27)

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

Рассмотрим также возмущенную систему с вектором правой части , , . Первое уравнение и определенная им прямая остаются неизменными, а график второй прямой смещается вдоль оси вверх на расстояние . (см. Рис.2.1). Несмотря на малость смещения второй прямой, точка пересечения прямых перемешается с исходной позиции на весьма существенное расстояние: , . Учитывая неравенство (2.21) , мы можем оценить снизу число обусловленности матрицы задачи:

. (2.28)

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

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

Примером плохо обусловленной матрицы может служить матрица Гильберта

 

. (2.29)

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

Итерационные методы решения систем ЛАУ

 

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

 

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

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

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

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

(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) достичь существенного уменьшения числа обусловленности матрицы по сравнению с числом обусловленности исходной матрицы .

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

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

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

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

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

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

:

, , .

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

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

 

Проблема собственных значений и собственных векторов.

 

Свойства собственных значений и собственных векторов матриц.

Число называется собственным значением матрицы , если существует такой ненулевой вектор , удовлетворяющий равенству . (4.1) Вектор , удовлетворяющий равенству (4.1), называется собственным вектором матрицы . Совокупность всех собственных…

Степенной метод.

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

Метод вращений.

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

Литература

1. А.А.Самарский, А.В.Гулин. Численные методы. М., Наука, 1989.

2. В.В.Воеводин, Ю.А.Кузнецов. Матрицы и вычисления. М., Наука, 1984.

3. Д.К.Фадеев, В.Н.Фадеева. Вычислительные методы линейной алгебры. СПт. 2002.

 

П.1. Реализация итерационного метода вращений для расчета собственных значений симметричной матрицы

% Проблема собственных значений % Метод Вращений N=22; % Задание Размерности матрицы

Содержание

Введение………………………………………………………………………………………..3

1. Задачи линейной алгебры…………………………………………………………………4

1.1 Основные понятия ………………………………………………………………..4

1.2 Нормы векторов и матриц………………………………………………………...5

1.3. Системы линейных алгебраических уравнений. Разрешимость

и устойчивость………………………………………………………………………….6

2. Прямые методы решения систем ЛАУ……………………………………………………..7

2.1 Метод Гаусса………………………………………………………………………..8

2.2 Метод Гаусса с выбором ведущего элемента…………………………………….9

2.3 LU-факторизация………………………………………………………………….10

2.4 Разложение Холецкого (метод квадратного корня)…………………………….12

2.5 Число обусловленности матрицы и оценки погрешности решения систем линейных алгебраических уравнений………………………………………………..13

2.6 Геометрический смысл и примеры плохой обусловленности матриц………...14

3. Итерационные методы решения систем ЛАУ

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

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

3.3 Метод простой итерации………………………………………………………….19

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

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

3.6 Методы сопряженных градиентов………………………………………….…….23

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

4. Проблема собственных значений и собственных векторов.

4.1 Свойства собственных значений и собственных векторов матриц……………..27

4.2 Степенной метод……………………………………………………………………29

4.3 Метод вращений……………………………………………………………………30

Литература ……………………………………………………………………………...31

Приложение

 

[*] Скалярное произведение может быть определено также другими способами в соответствии с аксиомами скалярного умножения.

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