Дж. Деммель - Вычислительная линейная алгебра, страница 13

Заметим, что к(А) = ЙА(! '2А 4)~ 4; таким образом, матрица А хорошо обусловлена и мы могли бы ожидать, что системы Ах = Ь решаются с хорошей точностью. 4) , число Я(1/10 ) равно 10, 01 -4, 4 Ю вЂ” 4 104 1), число Я(1 — 10 1) Равно — 10, Ю-4 1 10-4 таким образом, ИУ = 104 2.4. Анализ ошибок Заметим, что исходный элемент аз» был полностью «утрачен» для вычислений, когда из него пришлось вычесть число 104. Мы получили бы те же самые множители 1 и сГ, будь азз равен 1, О или — 2; вообще, любому числу, такому, что Я(азз — 104) = — 104. Поскольку далее алгоритм работает только с 1 и У, будет получен один и тот же ответ для всех этих различных значений азз, которые соответствуют совершенно различным матрицам А, а потому совершенно различным векторам х = А 'Ь; мы не можем, следовательно, гарантировать хорошую точность результата.

Это явление называется численной неустойчивостью, поскольку Ь и У не являются точными треугольными множителями матрицы, близкой к А. (По-другому, можно сказать, что число ОА — Х Ц! сравнимо с ОАО вместо того, чтобы быть величиной порядка е!)АО.) Посмотрим, что случится, если решать систему Ах = (1, 2)т, пользуясь вычисленным разложением. Правильный ответ: х (1,1] . Вместо него мы пот лучим следующее. Решая систему 1у = [1,2)~, находим у1 —— Я(1/1) = 1 и уз = Я(2 — 104 1) = — 104; отметим, что значение 2 было утеряно при вычитании из него числа 104.

Решение системы Ух = у дает хт — — Я(( — 104)/ * — 104)) = 1 и х1 — — Я((1 — 1)/10 4) = О, т. е. совершенно ошибочный результат. Другое предупреждение о потере точности можно получить, сравнивая число обусловленности матрицы А с числами обусловленности / и У. Напомним, что задача решения системы Ах = Ь была преобразована в решение двух других систем с матрицами 1 и У, поэтому мы не хотим, чтобы Е и У были обусловлены много хуже, чем А. Между тем, здесь число обусловленности матрицы А близко к 4, тогда как числа обусловленности для 1 и У равны 10э. В следующем разделе мы покажем, что использование метода СЕРР почти всегда исключает неустойчивость, иллюстрируемую данным примером.

Если бы здесь был применен этот метод, то прежде всего был бы изменен порядок двух уравнений. Предлагаем читателю проверить, что тогда были бы получены матрицы 01 ~ 1 01 Я(.0001/Ц 1 ~ = ~ .0001 0 Я(1 †.0001 1) 0 1 произведение которых весьма близко к А. Обе матрицы 1 и У (как и А) обусловлены очень хорошо. Очень хорошую точность имеет и вычисленное решение системы.

2.4.2. Формальный анализ ошибок гауссова исключения Изложим интуитивные соображения, стоящие за нашим анализом ошибок 1Л1- разложения. Если промежуточные величины, возникающие в произведении 1 У, очень велики по сравнению с ОАО, то информация, содержащаяся в элементах матрицы А, будет «утрачена», когда из них будут вычтены эти большие величины. Именно это произошло с элементом азз в примере из разд. 2.4.1. Если бы, напротив, промежуточные величины в произведении 1 У были сравнимы с элементами матрицы А, то можно было бы рассчитывать на получение Глава 2.

Решение линейных уравнений разложения с малой обратной ошибкой А — ЬП. Следовательно, мы хотим оценить наибольшие промежуточные величины в произведении Ь. 11. Мы сделаем это, оценивая элементы матрицы Д ~Ц (по поводу обозначений см. равд. 1.1). Наш анализ будет схож с тем, что мы проделали в равд. 1.б для задачи вычисления многочлена. Рассматривая там р = 2 '1 агх', мы показали, что если число (р) сравнимо с суммой 2 '„)агхз), то р будет вычислено с хорошей точностью. После того как будет дан общий анализ гауссова исключения, мы используем его для демонстрации того, что метод СЕРР (или, при ббльших затратах, метод СЕСР) обеспечивает сравнимость элементов матрицы Д Щ с 9А9 почти в любой реалистической ситуации.

К сожалению, наилучшие оценки для йоА9, которые мы сможем доказать, все же, в общем случае, много болыпе ошибок, наблюдаемых в реальных задачах. Поэтому оценки ошибок, используемые на практике, основываются на вычисленной невязке г и оценке (2.5) (или оценке (2.9)), а не на строгой, но пессимистичной оценке, выводимой в данном разделе. Чтобы упростить обозначения, предположим, что выбор главных элементов в матрице А проделан заранее.

Упростим алгоритм 2.2 до двух соотношений; одно из них относится к элементам ауе с у < Й, а другое — к элементам а а с у > й. Проследим вначале, что происходит в алгоритме 2.2 с элементом первого типа: этот элемент перевычисляется для каждого значения 1 от 1 до у — 1 путем вычитания 1;иге, после чего его значение присваивается элементу и а; таким образом, У вЂ” 1 и.е = ауа — ~~1 1зчизе.

Если у > и, то и здесь из а а вычитается 1чиге для 1 от 1 до а — 1, затем полученная сумма делится на иеа и результат становится значением элемента т Š— 1 а а — ~г .1 1 1 зигй иее Чтобы провести анализ ошибок округлений для этих двух формул, воспользуемся результатом, установленным в вопросе 1.1бя скалярное произведение, вычисленное в арифметике с плавающей точкой, удовлетворяет соотношению Я (~~ хгрз ~ = ~~ хзуг(1+61), где ~б1~ < 11е.

/ 1=1 Применим этот результат к формуле для и е, что даст1 1-1 и.ь = а е — ~~ 1чиге(1+ 61) (1+ б'), 1=1 СтрогО говоря, наш анализ формулы предполагает, что сначала вычисляется сумма, а затем накопленный результат вычитается из ауе. Однако окончательная оценка не зависит от порядка суммирования. 57 2А. Анализ ошибок где (61! < Ц вЂ” 1)е и (б'~ < е. Разрешая это равенство относительно а1зо получаем ,а-1 1 а ь =,и ь .1 . + ~~! 1,иаа(1+ 61), так как 1. = 1, 1 1 зися 1 + ~~ 1,чи ьб! а=1 1 и 1+61 =,) (где ~б,) < (1 — 1)е Величина Езь оценивается так: ~Е1а( = ~111 иаь б, < ~~ Яа~.

)иаь(.пе = пе()Ь). )У!)зр,. Проводя аналогичный анализ формулы для 1 ь, находим (1+ б')(абь — ~'„1 11,иаав(1+ ба)) 1,1 =(1+бн) где ~61) < (Ь вЂ” 1)е, (б'( < е и ~б"! < е. Разрешая это соотношение относительно а;и, имеем 1 1 †иьь111 + ~~ 11;иа (1+ 61) '1 + б!) '1 + ба!) а=1 1!1и,1а + ~~ 11,нагиба а=1 а=1 11и11 + Езь, 1 (1 + б') (1 + б" ) ' где (61! < ие. Итак, как и выше, ~Е ь) < пеЩ . )У( ы В целом, итоги этого анализа ошибок можно суммировать простой формулой А = 111+ Е, где ~Е! < пеф )Ц. Беря нормы, получаем (((ЕЯ < пе(ЩШ' )!)Ц)!. Если норма зависит только от абсолютных величин элементов матрицы (зто верно для шах-нормы, 1-нормы, нормы айробениуса, но не для 2-нормы), то это неравенство можно упростить к виду ()Е~! < ис)Щ ~Щ!.

Проанализируем теперь оставшуюся часть задачи: решить систему Шх = Ь, решая системы Еу = Ь и Ух = р. Согласно результату, полученному в вопросе 1.11, решение системы Ьу = Ь посредством прямой подстановки дает приближенное решение у, удовлетворяющее уравнению (А+ 61 ) р = Ь, где (Ы~ < ив ~Ц. Аналогично, при решении системы Ух = у мы получим вектор х, удовлетворяющий уравнению (1а'+ 61а')х = у с (бУ( < пе)У!.

58 Глава 2. Решение линейных уравнений Обьединяя эти соотношения, имеем Ь = (Х+6Х)д = (Ь+ 6Х)(У+ 6У)х = (ХХХ+ Ь6У+ 6Х1Х+ 6Х6У)х = (А — Е+ Х03+ 6ЬУ+ 6Х6ХХ)х = (А+6А)х, где 6А = — Е+ Ь6У+ 6ХУ+6Ь6У. Оценим теперь 6А, используя неравенство треугольника и наши оценки для Е, 6Х и6У: (6А~ = ) — Е+ Х6У+ 6Ы1 + 6Ь6Ц < )Е( + )Ь6Ц + )6Х Ц + )6Х 6Ц < (Е( + Д (6Ц + )6Х ) (Ц + (6Х! (6Ц < пг(Х) )Ц+пгД.

(Ц+пг~Х~. (Ц+и в '1Ц (Ц Зива )Ц. Беря здесь нормы и предполагая, что ЯХ)!) = '8Х(! (как и прежде, это верно для шах-нормы, 1-нормы, нормы Фробениуса, но не для 2-нормы), получаем '86А8 < ЗпгЗЩ ))Ц). Таким образом, чтобы понять, в каких случаях гауссово исключение является обратно устойчивым алгоритмом, нужно выяснить, когда верно соотношение ЗпвбЩ Щ! = О(в)0А0. Именно тогда величина "'(Гл)' в оценках теоРии возмущений имеет порядок О(в), как мы того желали (заметим, что 66 = 0), Основное эмпирическое наблюдение, поддерживаемое десятилетиями вычислительной практики, заключается в том, что метод СЕРР почти всегда обеспечивает соотношение Щ~ '0Щ 'ЗА8. Метод гарантирует, что каждый элемент в Л ограничен по абсолютной величине единицей, поэтому достаточно рассмотреть ЗЦ~.

Определим коэффициент роста в методе СЕРР' как дрр = !1П ак/!~А!! ак, где !)Ац„, = шахгу)а11(; таким образом, устойчивость метода эквивалентна тому, что число дрр мало или является медленно растущей функцией от п. На практике, дрр почти всегда не превосходит и. Поведение в среднем, по-видимому, соответствует функции пзуз или, возможно, даже п11г (242]. (См. рис.

2.1.) Это делает метод СЕРР наиболее предпочтительным выбором для многих задач. К сожалению, существуют немногочисленные примеры, где дрр может достигать величины 2" Предложение 2.1. Метод СЕРР гарантирует, что дрр < 2" 1. Эта оценка достижима. Доказательство. На первом шаге метода происходит перевычисление ада = ада — 1чиип где )111! < 1 и 1ига1 = 1а1а! < п1ах„, (аг,(; отсюда ~ауа( < 2 шахте ~ага~. Таким образом, каждый из и — 1 основных шагов метода может удвоить величину оставшихся матричных элементов, что приводит к общей оценке 2" Пример в вопросе 2.14 показывает, что эта оценка достижима.

П Это определение немного отличается от определения, обычно используемого в литературе, но, по существу, эквивалентно ему (см.(121, с.115)). 39 2.4. Анализ ошибок Собирая все приведенные выше оценки, получаем !)вАО, < ЗдррпзвЗАО (2.11) Здесь использованы неравенства ОД( < и и йЦ) < пдррОАЙ . Множитель Здррпз в оценке (2.11) заставляет ее почти всегда сильно переоценивать подлинное значение 96АЙ, даже если дрр = 1.

Например, при в = 10 т и и = 150, что соответствует матрице весьма умеренного порядка, имеем Зпзв > 1, что означает: потенциально все верные знаки в вычисленном решении могут быть потеряны. В примере 2.3 сравниваются графики функции Здррпзв и подлинной обратной ошибки, чтобы показать, насколько пессимистичен может быть первый из них. Величина ОвАО обычно имеет порядок 0(в) ОАЙ, поэтому можно сказать, что на практике метод СЕРР обратно устойчив, хотя мы и можем построить примеры, на которых он «проваливается».

В разделе 2.4.4 представлены практичные оценки ошибки в вычисленном решении системы Ах = Ь, много меньшие, чем оценка, получаемая при использовании неравенства йвАО, < ЗдррпэвйА$~, Можно показать, что метод СЕСР еще более устойчив, чем метод СЕРР. Оценка для коэффициента роста дор этого метода в наихудшем случае имеет вид (2б2, с. 213] "~ *'~ < и. 2. 3»!з . 4»!э,п»Л -П - и'/зы«а. в/4 гпахб (а,.! Эта верхняя граница также слишком велика для практических целей. Поведение коэффициента дор в среднем соответствует функции и'~~.