С формальной точки зрения каждый шаг исключения представляет собой умножение слева матрицы системы на специфическую элементарную матрицу (разд, 2.2,1,5), Например, элементарная матрица, осуществляющая первый шаг исключения для системы (2.36). такова:

М, = 1 -

1/2 \2 /

1 О 0) = ( -1/2

0 0\

о ij

Заметим, что Mi - единичная нижняя треугольная матрица, которая отличается от диагональной только первым столбцом (столбцом с номером шага исключения), где есть ненулевые субдиагональные элементы.

Итак, если в случае с п неизвестными описанный процесс гауссовых исключений не прервется из-за того, что какой-то ведущий элемент окажется нулем, то этот процесс можно будет проинтерпретировать как поэтапное умножение А слева на пакет [М} из п-1 единичных нижних треугольных матриц с целью получить в результате верхнюю треугольную матрицу U:

жл = ж„ ,. . .ММА = и. (2.38)

Каждая Mi преобразует 1-й столбец частично триангуляризованной матрицы к желаемому виду Следует подчеркнуть, что формализм описания процедуры исключений с помощью матриц {Mi) вовсе не предполагает запоминания их целиком при реализации алгоритма на машине; 1-я матрица полностью задается п - i множителями, так что исчерпывающие данные о преобразовании поместятся, например, в освобождающийся поддиатональный блок записи А.

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

MAx=Ux=Mb,

(2.39)

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

Поскольку матрица М невырождена, равенство (2.38) можно переписать так:

АМи=М:ГМ7\. .М-1,и. (2.40)

В силу специальной структуры матриц М обратные к ним матрицы также будут единичными нижними треугольными. Читатель без труда убедится, что Му отличается от только знаками элементов /-го столбца. Соответственно, если переобозначить M~ через L, равенство (2.40) принимает вид представления А произведением единичной нижней треугольной матрицы L на

верхнюю треугольную матрицу U:

A=LU.

Это представление называется LU-разложением А.

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

(О вычислить L и и, такие, что A=LU (ведущие элементы предполагаются ненулевыми);

(ii) решить систему Ly=b (что эквивалентно формированию У=МЬУ,

(iii) решить систему Ux=y.

Если не принимать во внимание малых п, то можно сказать, что основной объем вычислений по методу гауссовых исключений приходится на этап формирования L и U. Для этого требуется выполнить Vs" умножений, в то время как сложность второго и третьего этапов характеризуется умножениями.

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

/о 1 1\

Л=[ 2 3 4 . Vl о 7j

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

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

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

в том, что при выборе очередной, k-vt ведущей строки предпочтение отдается той из еще не побывавших ведущими строк, в которой модуль элемента из fe-ro столбца максимален. Она и становится fe-й, «поменявшись местами» со строкой, которая занимала это место ранее (а указанный элемент становится ведущим). Например, если бы частичное упорядочение использовалось на первом шаге триангуляризации матрицы (2.37), то строки 1 и 3 следовало бы переставить, получив тем самым матрицу

/4 -2 3\ Л = ( I 3 2

Подобные перестановки возможны на каждом шаге. Формально они описываются умножением на перестановочную матрицу, элементы которой могут быть только нулями или единицами, причем в каждом столбце и каждой строке есть в точности по одной единице. Умножение некоторой матрицы на перестановочную слева меняет местами строки, справа - столбцы. Например, перестановка строк 1 и 3 матрицы А в (2.37) достигается умножением А слева на

/о О Г

1 О О О

Дополнение гауссовых исключений описанной выше стратегией выбора ведущей строки означает, что всякий раз, когда потребуется поменять строки местами, частично триангуляризован-ная матрица будет умножена слева на соответствующую перестановочную. В результате это приведет к /,г7-разложению матрицы, которая получается из исходной перестановкой строк. Таким образом, будет справедливо равенство РА = Ш, где Р - некоторая перестановочная матрица. При этом модули всех поддиагональных элементов матрицы L будут лежать между нулем и единицей. Последнее вытекает из способа выбора ведущего элемента; его модуль всегда максимален в рассматриваемой части очередного столбца, и, стало быть, модули всех множителей (отношений прочих элементов этой части столбца к ведущему) меньше или равны единице.

Обратный анализ ошибок для гауссовых исключений с частичным упорядочением - довольно громоздкая процедура. Поэтому ограничимся сводкой результатов; численно найденные L и U будут точными треугольными сомножителями разложения некоторой матрицы Р(Л+£). Введем величину

S max I Oij I

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

l£ll.<n"VEMgMIL, где Y - константа порядка единицы, не зависящая ни от Л, ни от п, а E.i - машинная точность.

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

2.2.5.2. LDL-разложение и разложение по Холесскому. Любая положительно определенная матрица Л может быть представлена произведением вида

ALDL. (2.41)

где L - единичная нижняя треугольная матрица, а D - диагональная матрица со строго положительными диагональными элементами. Это представление называется LDL-разложением.

Поскольку диагональ D положительна, равенство (2.41) можно переписать так;

А-ШЧЮЧг1 = Ш = НтН. (2.42)

Здесь L - невырожденная нижняя треугольная матрица общего вида, а R - невырожденная верхняя треугольная матрица общего вида. Это представление принято называть разложением по Холесскому, а матрицу R - фактором Холесского или «квадратным корнем») из Л (по аналогии со скалярным случаем). Впредь, учитывая близость (2.41) и (2.40), мы будем называть разложением по Холесскому и то и другое представление.

Факторы Холесского могут быть вычислены путем прямого поэлементного сопоставления. По определению имеем

Сц 012 013

aia оая оаз 013 Озз Сзз .

П1„ > Qan

(2.43)

Соответственно, приравнивая (1, 1)-й элемент матрицы слева и отвечающий ему элемент произведения справа, получим Cn/fi, fxiVau. После этого, сравнивая первые строки левой и правой

частей (2.43), найдем все оставшиеся элементы первой строки R. В самом деле, это сравнение показывает, что

ain=rn>4„.

Далее, элементы (2, 2) матриц R и А удовлетворяют равенству аг1=г\+4. Отсюда, зная г,, находим rg, затем вторую строку R и т. д. Приведенный алгоритм иногда называют пестренной факторизацией Холесского. Элементы матрицы R можно вычислять и в таком порядке: /„, ги, Л22, п,, rs, г,,, .... По понятным причинам соответствующую процедуру принято называть столбцовой факторизацией.

Для построения разложения по Холесскому пхп-матрицы требуется выполнить около /п операций умножения (сложения) и вычислить п квадратных корней (последнее - только для формы (2.42)).

Замечательная особенность алгоритма Холесского состоит в том, что в отличие от гауссовых исключений он численно устойчив без каких-либо перестановок. Это свойство определяется следующим соотношением между элементами А и R:

*-1--й+...---м = ам, k=\.....п.

В силу него существует априорное ограничение сверху на модули элементов R: каждый из них не превосходит максимальной величины Gftfe. Соответственно «лавинообразный рост» элементов R невозможен независимо от того, будут ведущие элементы (в данном случае естественно назвать так элементы г„) малыми или нет. Обратный анализ ошибок дает для численной реализации алгоритма Холесского формуЬу RR=A+E с ограничением на 1£ сверху, аналогичным полученному для гауссовых исключений, но не содержащим фактора роста.

Было бы заманчиво предположить, что любую симметричную матрицу можно представить в виде (2.41), если не оговаривать знаки диагональных элементов D. К сожалению, что следует подчеркнуть особо, данная гипотеза ошибочна. Хля произвольной наперед взятой знакоиеопределенной симметричной матрицы LDL-представление может оказаться невозможным; например, его не существует для матрицы

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

разложения

/е 1\П Оув О yi

становятся очень большими.

2.2.5.3. ОЛ-разложение. В разд. 2.2.5.1 был изложен способ преобразования матриць! в верхнюю треугольную применением последовательности элементарных нижних треугольных матриц. Аналогичного результата можно достичь и по-другому, используя набор элементарных ортогональных матриц.

Для любого ненулевого вектора w определено так называемое преобразование Хаусхолдера с элементарной симметричной матрицей

где р=/211и). Эта матрица ортонормальна и, следовательно, преобразует векторы с сохранением евклидовой нормы. Для двух произвольных несовпадающих векторов с и Ь одинаковой евклидовой длины всегда можно подобрать матрицу Хаусхолдера, которая переведет один вектор в другой, т. е. обеспечит равенство

На = I-,

)a = a-w[-j-

(2.44)

Легко проверить, что оно справедливо при любом w, коллинеарном разности а - Ь.

Заметим также, что воздействие преобразования Н на вектор а сводится к вычитанию из а вектора w с некоторым множителем. Соответственно те компоненты исходного вектора, которым отвечают нулевые компоненты вектора Хаусхолдера w, при этом преобразовании не изменяются. Более того, преобразование Хаусхолдера вообще не изменит вектора а, если окажется, что wa=0.

Среди всевозможных ортогональных преобразований выделяют так называемые плоские повороты. Последаше часто используются для обнуления какой-нибудь одной компоненты вектора. Плоский поворот задается матрицей Хаусхолдера с вектором w, у которого есть только две ненулевые компоненты. Если они стоят в j-й и /-Й позициях, эта матрица будет отличаться от единичной только 1-й и j-й строками и столбцами. При i, меньшем /, ее элементы (i, i), (, Л. О, О и (/, /) образуют конфигурацию вида

где c+s=\ (так что c=cos 6 и s=sin 6 для некоторого в).

Умножение вектора слева на рассматриваемую матрицу плоского поворота изменит только i-ю и /-ю компоненты вектора.