1611688890-f641c9ec8276824e4686da772eb56520 (826652), страница 50
Текст из файла (страница 50)
Недостатком этого подхода являетсясущественная опора на LU-разложение матрицы, и потому желательноиметь более прямой способ нахождения разложения Холесского.Выпишем равенство A = CC ⊤ , определяющее множитель Холес-3063. Численные методы линейной алгебрыского, в развёрнутой форме с учётом симметричности A:a11 a21 .. .an1c11 c21 . . .cn1a22...an2..cn2···▽...... 0c22....cnnc11ann=c21c22······... · 0cik cjkпри i ≥ j.cn1cn2 ..
,. cnn(3.74)где «▽» означает симметричные относительно главной диагонали элементы матрицы, которые несущественны в последующих рассмотрениях. Можно рассматриваеть это равенство как систему уравнений относительно неизвестных переменных c11 , c21 , c22 , . . . , cnn — элементовнижнего треугольника множителя Холесского. Всего их 1+2+. . .+n =12 n(n + 1) штук. Для их определения имеем столько же соотношений,вытекающих в матричном равенстве (3.74) из выражений для элементов aij , i ≥ j, которые образуют нижний треугольник симметричнойматрицы A = (aij ).В поэлементной форме система уравнений (3.74) имеет вид, определяемый правилом умножения матриц и симметричностью A:aij =jXk=1(3.75)Выписанные соотношения образуют, фактически, двумерный массив, вкотором уравнения имеют двойные индексы — i и j, но их можно линейно упорядочить таким образом, что система уравнений (3.75) получитспециальный вид, очень напоминающий треугольные СЛАУ.
Далее этасистема может быть решена с помощью процесса, сходного с прямойподстановкой для треугольных СЛАУ (см. §3.6а).В самом деле, если выписывать выражения для элементов aij постолбцам матрицы A, начиная в каждом столбце с диагонального элемента ajj и идя сверху вниз до ajn , то все уравнения из (3.75) разбиваются на n следующих групп, которые удобно занумеровать столбцовым3073.6.
Прямые методы решения линейных системиндексом j = 1, 2 . . . , n:для j = 1(для j = 2(c221 + c222 = a22 ,для j = 3(c231 + c232 + c233 = a33 ,···c211 = a11 ,ci1 c11 = ai1 ,i = 2, 3, . . . , n,ci1 c21 + ci2 c22 = ai2 ,i = 3, 4, . . . , n,ci1 c31 + ci2 c32 + ci3 c33 = ai3 ,···i = 4, 5, . . . , n,.В краткой записи получающаяся система может быть записана следующим образом (c2j1 + c2j2 + .
. . + c2j,j−1 + c2jj = ajj ,ci1 cj1 + ci2 cj2 + . . . + cij cjj = aij ,j = 1, 2, . . . , n.i = j + 1, . . . , n,(3.76)где считается, что cji = 0 при j < i.Получается, что в уравнениях из (3.76) для j-го столбца множителяХолесского присутствуют все элементы j-го и предшествующих столбцов. Если последовательно рассматривать группы уравнений в порядкевозрастания номера j, то реально неизвестными к моменту обработкиj-го столбца (т. е.
решения j-ой группы уравнений) являются только(n − j + 1) элементов cij именно этого j-го столбца, которые к томуже выражаются несложным образом через известные элементы и другчерез друга.В целом выписанная система уравнений (3.76) действительно имееточень специальный вид, пользуясь которым можно находить элементыcij матрицы C последовательно друг за другом по столбцам в порядке,3083.
Численные методы линейной алгебрыC=↓↓↓↓...y↓...y0......···↓y×,Рис. 3.16. Схема определения элементов треугольногомножителя в методе Холесского.который наглядно изображён на Рис. 3.16. Более точно,(√c11 = a11 ,при j = 1ci1 = ai1 /c11 ,i = 2, 3, .
. . , n,при j = 2(при j = 3(pa22 − c221 ,= ai2 − ci1 c21 /c22 ,c22 =ci2c33 =pa33 − c231 − c232 ,i = 3, 4, . . . , n,ci3 = ai3 − ci1 c31 − ci2 c32 /c33 ,i = 4, 5, . . . , n,и так далее для остальных j. Псевдокод этого процесса приведён вТабл. 3.1, где считается, что если верхний предел суммирования превосходит нижний, то сумма «пуста» и суммирование не выполняется.Если A — симметричная положительно определённая матрица, то всилу Теоремы 3.6.4 система (3.76) обязана иметь решение, и наш алгоритм успешно прорабатывает до конца, находя его. Если же матрицаA не является положительно определённой, то алгоритм (3.77) аварийно прекращает работу при попытке извлечь корень из отрицательногочисла либо разделить на нуль.
Вообще, запуск алгоритма (3.77) — этосамый экономичный способ проверки положительной определённостисимметричной матрицы.Метод решения СЛАУ, основанный на разложении Холесского и использующий соотношения (3.73) и алгоритм (3.77), называют методомХолесского.
Он был предложен в 1910 году А.-Л. Холесским в неопубликованной рукописи, которая, тем не менее, сделалась широко извест-3093.6. Прямые методы решения линейных системТаблица 3.1. Алгоритм разложения Холесского(прямой ход метода Холесского)DO FOR j = 1 TO nvuj−1uXc ← ta −c2jjjjjkk=1DO FOR i = j + 1 TO ncij ←END DOaij −j−1Xk=1cik cjk!(3.77)cjjEND DOной во французской геодезической службе, где решались такие системы уравнений.
Позднее метод неоднократно переоткрывался, и потомуиногда в связи с ним используются также термины «метод квадратногокорня», «метод квадратных корней» или даже другие имена, данныеего позднейшими авторами.Метод Холесского можно рассматривать как специальную модификацию метода Гаусса, которая требует вдвое меньше времени и памяти ЭВМ, чем обычный метод Гаусса в общем случае. Замечательнымсвойством метода Холесского является также то, что обусловленностьмножителей Холесского, вообще говоря, является лучшей, чем у матрицы исходной СЛАУ: она равна корню квадратному из обусловленностиматрицы системы (это следует из самого разложения Холесского).
Тоесть, в отличие от обычного метода Гаусса, треугольные системы линейных уравнений из (3.73), к решению которых сводится задача, менеечувствительны к ошибкам, чем исходная линейная система. В следующем пункте мы увидим, что подобную ситуацию следует рассматриватькак весьма нетипичную.Если при реализации метода Холесского использовать комплекснуюарифметику, то извлечение квадратного корня можно выполнять все-3103. Численные методы линейной алгебрыгда, и потому такая модификация применима к симметричным неособенным матрицам, которые не являются положительно определёнными.
При этом множители Холесского становятся комплексными треугольными матрицами.Другой способ распространения идеи метода Холесского на произвольные симметричные матрицы состоит в том, чтобы ограничитьсяразложением (3.72), которое называется LDL⊤-разложением матрицы.Если исходная матрица не является положительно определённой, тодиагональные элементы в матричном множителе D могут быть отрицательными. Но LDL⊤-разложение столь же удобно для решения системлинейных алгебраических уравнений, как и рассмотренные ранее треугольные разложения. Детали этих построений читатель может найти,к примеру, в [11, 15, 43, 67].Отметим также, что существует возможность другой организациивычислений при решении системы уравнений (3.75), когда неизвестныеэлементы c11 , c21 , c22 , .
. . , cnn последовательно находятся по строкаммножителя Холесского, а не по столбцам, как в (3.77). Этот алгоритмназывается схемой окаймления [15], и он по своим свойствам примерноэквивалентен рассмотренному выше алгоритму (3.77).3.7Прямые методы на основеортогональных преобразований3.7аЧисло обусловленностии матричные преобразованияПусть A и B — неособенные квадратные матрицы, и матрица Aумножается на матрицу B. Как связано число обусловленности произведения AB с числами обусловленности сомножителей A и B?Справедливы соотношенияи поэтомуkABk ≤ kAk kBk, (AB)−1 = B −1 A−1 ≤ kA−1 k kB −1 k,cond(AB) = (AB)−1 kABk ≤ cond A · cond B.(3.78)3.7. Методы на основе ортогональных преобразований311С другой стороны, если C = AB, то A = CB −1 , и в силу доказанногонеравенстваcond(A) ≤ cond(C) · cond(B −1 ) = cond(AB) · cond(B),коль скоро cond(B −1 ) = cond(B).
Поэтомуcond(AB) ≥ cond(A)/cond(B).Аналогичным образом из B = CA−1 следуетcond(AB) ≥ cond(B)/cond(A).Объединяя полученные неравенства, в целом получаем оценкуcond(A) cond(B).(3.79),cond(AB) ≥ maxcond(B) cond(A)Ясно, что её правая часть не меньше 1.Неравенства (3.78)–(3.79) кажутся грубыми, но они достижимы. Всамом деле, пусть A — неособенная симметричная матрица с собственнными значениями λ1 , λ2 , . . . и спектральным числом обусловленности,равным (см.
стр. 276)cond2 (A) =maxi |λi (A)|.mini |λi (A)|У матрицы A2 собственные векторы, очевидно, совпадают с собственными векторами матрицы A, а собственные значения равны λ21 , λ22 , . . . .Как следствие, числом обусловленности матрицы A2 становится2maxi |λi (A)|2maxi |λi (A)|maxi (λi (A))2,==cond2 A2 =mini (λi (A))2mini |λi (A)|2mini |λi (A)|и в верхней оценке (3.78) получаем равенство. Совершенно сходным образом можно показать, что для спектрального числа обусловленностиоценка (3.78) достигается также на произведениях вида A⊤A.Нижняя оценка (3.79) достигается, к примеру, при B = A−1 длячисел обусловлености, порождённых подчинёнными матричными нормами.Практически наиболее важной является верхняя оценка (3.78), иона показывает, что при преобразованиях и разложениях матриц число обусловленности может существенно расти.