1.5. Решение системы методом Холецкого.
Если матрица системы является симметричной и положительно определенной, то для решения системы применяют метод Холецкого (метод квадратных корней).
В основе метода лежит алгоритм специального LU-разложения матрицы A, в результате чего она приводится к виду A=
.
Если разложение получено, то как и в методе LU-разложения, решение системы сводится к последовательному решению двух систем с треугольными матрицами:
и 
.
Для нахождения коэффициентов матрицы L неизвестные коэффициенты матрицы L приравнивают соответствующим элементам матрицы A. Затем последовательно находят требуемые коэффициенты по формулам:
Пример 5. Решение системы методом Холецкого.
A= 
b=
Находим элементы матрицы L:
Таким образом разложение матрицы A имеет вид:
Последовательно решаем системы 
и
. Решением 1-ой системы является вектор
,
а решение 2-ой системы вектор 
.
Ответ: 
Решение с использованием встроенной функции системы Mathcad.
Пример 6.
1.6.Метод прогонки.
Рассмотрим систему уравнений с трех диагональной матрицей:
Преобразуем первое уравнение системы к виду
,
где 
,
Подставим полученное выражение во второе уравнение системы и преобразуем его к виду
На i-ом шаге уравнение преобразуется к виду
,
где 
,
.
На m-ом шаге подстановка в последнее уравнение выражения
дает возможность определить значение 
:
.
Значения остальных неизвестных находятся по формулам:
Вычислим прогоночные коэффициенты:
, 
, 
, 
, 
, 
, 
Обратный ход прогонки.
Находим значения неизвестных:
,
,
,
Ответ: 
.
2. Контрольные задачи.
Задача 1 . Дана система уравнений Ax=b.
ПОРЯДОК РЕШЕНИЯ ЗАДАЧИ.
Задать матрицу системы A и вектор правой части b.
Найти решение системы Ax=b.
а) используя встроенную функцию lsolve пакета MATHCAD,
б) с помощью метода Гаусса с выбором ведущего элемента по столбцу.
Пример выполнения фрагмента задачи
x:=lsolve(A,b) x=
Разложение Холецкого (метод квадратного корня)
| Разложение Холецкого | |
| Последовательный алгоритм | |
| Последовательная сложность | [math]O(n^3)[/math] |
| Объём входных данных | [math]\frac |
| Объём выходных данных | [math]\frac |
| Параллельный алгоритм | |
| Высота ярусно-параллельной формы | [math]O(n)[/math] |
| Ширина ярусно-параллельной формы | [math]O(n^2)[/math] |
Содержание
1 Свойства и структура алгоритма
1.1 Общее описание алгоритма
Первоначально разложение Холецкого использовалось исключительно для плотных симметричных положительно определенных матриц. В настоящее время его использование гораздо шире. Оно может быть применено также, например, к эрмитовым матрицам. Для повышения производительности вычислений часто применяется блочная версия разложения.
Для разреженных матриц разложение Холецкого также широко применяется в качестве основного этапа прямого метода решения линейных систем. В этом случае используют специальные упорядочивания для уменьшения ширины профиля исключения, а следовательно и уменьшения количества арифметических операций. Другие упорядочивания используются для выделения независимых блоков вычислений при работе на системах с параллельной организацией.
На этой странице представлено исходное разложение Холецкого с новых позиций нашего суперкомпьютерного века. Приведено описание конкретной версии разложения Холецкого — для плотных вещественных симметричных положительно определённых матриц, но структура для ряда других версий, например, для комплексного случая, почти такая же, различия состоят в замене большинства вещественных операций на комплексные. Список остальных основных вариантов разложения можно посмотреть на странице Метод Холецкого (нахождение симметричного треугольного разложения).
[math] L = \begin
[math] U = \begin
1.1.1 Симметричность и положительная определённость матрицы
Симметричность матрицы позволяет хранить и вычислять только чуть больше половины её элементов, что почти вдвое экономит как необходимые для вычислений объёмы памяти, так и количество операций в сравнении с, например, разложением по методу Гаусса. При этом альтернативное (без вычисления квадратных корней) LU-разложение, использующее симметрию матрицы, всё же несколько быстрее метода Холецкого (не использует извлечение квадратных корней), но требует хранения всей матрицы. Благодаря тому, что разлагаемая матрица не только симметрична, но и положительно определена, её LU-разложения, в том числе и разложение методом Холецкого, имеют наименьшее эквивалентное возмущение из всех известных разложений матриц.
1.1.2 Режим накопления
Алгоритм позволяет использовать так называемый режим накопления, обусловленный тем, что значительную часть вычислений составляют вычисления скалярных произведений.
1.2 Математическое описание алгоритма
Исходные данные: положительно определённая симметрическая матрица [math]A[/math] (элементы [math]a_
Вычисляемые данные: нижняя треугольная матрица [math]L[/math] (элементы [math]l_
Существует также блочная версия метода, однако в данном описании разобран только точечный метод.
1.3 Вычислительное ядро алгоритма
в режиме накопления или без него, в зависимости от требований задачи. Во многих последовательных реализациях упомянутый способ представления не используется. Дело в том, что в них вычисление сумм типа
в которых и встречаются скалярные произведения, ведутся не в порядке «вычислили скалярное произведение, а потом вычли его из элемента», а путём вычитания из элемента покомпонентных произведений, являющихся частями скалярных произведений. Поэтому следует считать вычислительным ядром метода не вычисления скалярных произведений, а вычисления выражений
в режиме накопления или без него.
Тем не менее, в популярных зарубежных реализациях точечного метода Холецкого, в частности, в библиотеках LINPACK и LAPACK, основанных на BLAS, используются именно вычисления скалярных произведений в виде вызова соответствующих подпрограмм BLAS (конкретно — функции SDOT). На последовательном уровне это влечёт за собой дополнительную операцию суммирования на каждый из [math]\frac
1.4 Макроструктура алгоритма
в режиме накопления или без него.
1.5 Схема реализации последовательного алгоритма
Последовательность исполнения метода следующая:
2. [math]l_
Далее для всех [math]i[/math] от [math]2[/math] до [math]n[/math] по нарастанию выполняются
1.6 Последовательная сложность алгоритма
Для разложения матрицы порядка n методом Холецкого в последовательном (наиболее быстром) варианте требуется:
Умножения и сложения (вычитания) составляют основную часть алгоритма.
При этом использование режима накопления требует совершения умножений и вычитаний в режиме двойной точности (или использования функции вроде DPROD в Фортране), что ещё больше увеличивает долю умножений и сложений/вычитаний во времени, требуемом для выполнения метода Холецкого.
При классификации по последовательной сложности, таким образом, метод Холецкого относится к алгоритмам с кубической сложностью.
1.7 Информационный граф
Опишем граф алгоритма [9] [10] [11] как аналитически, так и в виде рисунка.
Граф алгоритма состоит из трёх групп вершин, расположенных в целочисленных узлах трёх областей разной размерности.
Аргумент этой функции
Аргументы операции следующие:
Аргументы операции следующие:
Результат срабатывания операции является промежуточным данным алгоритма.
Интерактивное изображение графа алгоритма без входных и выходных данных для n = 4.
1.8 Ресурс параллелизма алгоритма
Для разложения матрицы порядка [math]n[/math] методом Холецкого в параллельном варианте требуется последовательно выполнить следующие ярусы:
При этом использование режима накопления требует совершения умножений и вычитаний в режиме двойной точности, а в параллельном варианте это означает, что практически все промежуточные вычисления для выполнения метода Холецкого в режиме накопления должны быть двойной точности. В отличие от последовательного варианта это означает увеличение требуемой памяти почти в 2 раза.
1.9 Входные и выходные данные алгоритма
Входные данные: плотная матрица [math]A[/math] (элементы [math]a_
Объём входных данных: [math]\frac
Выходные данные: нижняя треугольная матрица [math]L[/math] (элементы [math]l_
Объём выходных данных: [math]\frac
1.10 Свойства алгоритма
Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов, как хорошо видно, является квадратичным (отношение кубической к линейной).
При этом вычислительная мощность алгоритма, как отношение числа операций к суммарному объему входных и выходных данных – всего лишь линейна.
При этом алгоритм почти полностью детерминирован, это гарантируется теоремой о единственности разложения. Использование другого порядка выполнения ассоциативных операций может привести к накоплению ошибок округления, однако это влияние в используемых вариантах алгоритма не так велико, как, скажем, отказ от использования режима накопления.
Дуги информационного графа, исходящие из вершин, соответствующих операциям квадратного корня и деления, образуют пучки т. н. рассылок линейной мощности (то есть степень исхода этих вершин и мощность работы с этими данными — линейная функция от порядка матрицы и координат этих вершин). При этом естественно наличие в этих пучках «длинных» дуг. Остальные дуги локальны.
Наиболее известной является компактная укладка графа — его проекция на треугольник матрицы, который перевычисляется укладываемыми операциями. При этом «длинные» дуги можно убрать, заменив более дальнюю пересылку комбинацией нескольких ближних (к соседям).
Это явление обусловлено положительной определённостью матрицы. Среди всех используемых разложений матриц это наименьшее из эквивалентных возмущений.
2 Программная реализация алгоритма
2.1 Особенности реализации последовательного алгоритма
В простейшем (без перестановок суммирования) варианте метод Холецкого на Фортране можно записать так:
При этом для реализации режима накопления переменная [math]S[/math] должна быть двойной точности.
Отдельно следует обратить внимание на используемую в такой реализации функцию DPROD. Её появление как раз связано с тем, как математики могли использовать режим накопления вычислений. Дело в том, что, как правило, при выполнении умножения двух чисел с плавающей запятой сначала результат получается с удвоенными длинами мантиссы и порядка, и только при выполнении присваивания переменной одинарной точности результат округляется. Эта возможность даёт выполнять умножение действительных чисел с двойной точностью без предварительного приведения их к типу двойной точности. На ранних этапах развития вычислительных библиотек снятие результата без округление реализовали вставками специального кода в фортран-программы, но уже в 70-х гг. XX века в ряде трансляторов Фортрана появилась функция DPROD, реализующая это без обращения программиста к машинным кодам. Она была закреплена среди стандартных функций в стандарте Фортран 77, и реализована во всех современных трансляторах Фортрана.
Для метода Холецкого существует блочная версия, которая отличается от точечной не тем, что операции над числами заменены на аналоги этих операций над блоками; её построение основано на том, что практически все циклы точечной версии имеют тип SchedDo в терминах методологии, основанной на исследовании информационного графа и, следовательно, могут быть развёрнуты. Тем не менее, обычно блочную версию метода Холецкого записывают не в виде программы с развёрнутыми и переставленными циклами, а в виде программы, подобной реализации точечного метода, в которой вместо соответствующих скалярных операций присутствуют операции над блоками.
Особенностью размещения массивов в Фортране является хранение их «по столбцам» (быстрее всего меняется первый индекс). Поэтому для обеспечения локальности работы с памятью представляется более эффективной такая схема метода Холецкого (полностью эквивалентная описанной), когда исходная матрица и её разложение хранятся не в нижнем, а в верхнем треугольнике. Тогда при вычислениях скалярных произведений программа будет «идти» по последовательным ячейкам памяти компьютера.
Есть и другой вариант точечной схемы: использовать вычисляемые элементы матрицы [math]L[/math] в качестве аргументов непосредственно «сразу после» их вычисления. Такая программа будет выглядеть так:
СОДЕРЖАНИЕ
Заявление
Разложение Холецкого эрмитовой положительно определенной матрицы A является разложением вида
Положительные полуопределенные матрицы
Если эрмитова матрица A является только положительно полуопределенной, а не положительно определенной, то она все еще имеет разложение формы A = LL *, где диагональные элементы матрицы L могут быть равны нулю. Разложение не обязательно должно быть уникальным, например:
Разложение ЛПНП
Близким вариантом классического разложения Холецкого является разложение ЛПНП,
Разложение LDL связано с классическим разложением Холецкого формы LL * следующим образом:
Пример
Вот разложение Холецкого симметричной вещественной матрицы:
А вот его разложение Т ЛПНП :
Приложения
Линейный метод наименьших квадратов
Нелинейная оптимизация
Моделирование Монте-Карло
Фильтры Калмана
Обращение матрицы
Неэрмитова матрица B также может быть обращена с помощью следующего тождества, где BB * всегда будет эрмитовой:
Вычисление
Какой из приведенных ниже алгоритмов работает быстрее, зависит от деталей реализации. Как правило, первый алгоритм будет немного медленнее, потому что он обращается к данным менее регулярно.
Алгоритм Холецкого
Рекурсивный алгоритм начинается с i : = 1 и
На шаге i матрица A ( i ) имеет следующий вид:
Если теперь мы определим матрицу L i как
тогда мы можем записать A ( i ) как
Алгоритмы Холецкого – Банахевича и Холецкого – Краута.
Если мы выпишем уравнение
и, следовательно, следующие формулы для элементов L :
Для комплексных и вещественных матриц разрешены несущественные произвольные изменения знаков диагональных и связанных недиагональных элементов. Выражение под квадратным корнем всегда положительно, если A вещественно и положительно определено.
Для комплексной эрмитовой матрицы применяется следующая формула:
Любой из вариантов доступа позволяет при желании выполнять все вычисления на месте.
Устойчивость расчета
Предположим, что мы хотим решить хорошо обусловленную систему линейных уравнений. Если используется разложение LU, то алгоритм будет нестабильным, если мы не используем какую-то стратегию поворота. В последнем случае ошибка зависит от так называемого фактора роста матрицы, который обычно (но не всегда) невелик.
Разложение ЛПНП
Альтернативной формой, устраняющей необходимость извлекать квадратные корни, когда A является симметричным, является симметричная неопределенная факторизация
Следующие рекурсивные отношения применяются для записей D и L :
Это работает до тех пор, пока сгенерированные диагональные элементы в D остаются ненулевыми. Тогда разложение единственное. D и L реальны, если A реально.
Для комплексной эрмитовой матрицы A применяется следующая формула:
Опять же, шаблон доступа позволяет при желании выполнять все вычисления на месте.
Вариант блока
где каждый элемент в матрицах выше является квадратной подматрицей. Отсюда следуют аналогичные рекурсивные отношения:
Это включает в себя матричные произведения и явную инверсию, что ограничивает практический размер блока.
Обновление декомпозиции
Обновление первого ранга
знак равно А + Икс Икс * <\ Displaystyle <\ тильда <\ mathbf >> = \ mathbf + \ mathbf 
Вот функция, написанная в синтаксисе Matlab, которая реализует обновление первого ранга:
Понижение первого ранга
Код для обновления ранга один, показанный выше, может быть легко адаптирован для понижения ранга один: нужно просто заменить два добавления в назначении на вычитания r и L((k+1):n, k) на вычитания.
Добавление и удаление строк и столбцов
Если у нас есть симметричная и положительно определенная матрица, представленная в блочной форме как А <\ displaystyle \ mathbf >
и его верхний фактор Холецкого
Эти формулы можно использовать для определения фактора Холецкого после вставки строк или столбцов в любую позицию, если мы соответствующим образом установим размеры строки и столбца (в том числе равными нулю). Обратная задача, когда у нас есть
с известным разложением Холецкого
и желаете определить фактор Холецкого
матрицы с удаленными строками и столбцами, А <\ displaystyle \ mathbf >
дает следующие правила:
знак равно А ± Икс Икс * <\ Displaystyle <\ тильда <\ mathbf >> = \ mathbf \ pm \ mathbf
Доказательство для положительных полуопределенных матриц.
Доказательство ограничивающим аргументом
Приведенные выше алгоритмы показывают, что каждая положительно определенная матрица имеет разложение Холецкого. Этот результат можно распространить на положительный полуопределенный случай ограничивающим аргументом. Аргумент не является полностью конструктивным, т. Е. Не дает явных численных алгоритмов для вычисления факторов Холецкого. А <\ displaystyle \ mathbf >
Доказательство с помощью QR-разложения.
Обобщение
действуя на прямую сумму
