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

Из вопроса 1.10 известна следующая оценка для оепибок округления при вычислении т: 1(Ах — Ь) — Я(Ах — Ь)! < (и+1)е()А) )х/+1Ь1). (2.15) Поэтому в оценке (2.14) можно заменить !т/ на 1т! + (и + 1)е(/А/ 1х/ + 1Ь1) (так сделано в ЬАРАСК-программе вйеввх), а в оценке (2.13) — ЦтЦ на ЦтЦ + (и + 1)е(ЦАЦ . ЦхЦ + ЦЬЦ). Множитель п + 1 обычно чересчур велик и при желании может быть опущен.

В-третьих, из-за округлений в гауссовом исключении, примененном к очень плохо обусловленной матрице, могут быть получены настолько неточные матрицы 1 и У, что оценка (2.14) окажется слишком малой. Пример 2.6. Приведем пример, принадлежащий Кахану (%. КаЬап), который иллюстрирует трудность получения подлинно гарантированных оценок ошибки. В этом примере матрица А будет шочно вырожденной. Следовательно, любая оценка для ~-~(;)Е))( не должна быть меньше единицы, указывая, что в вычисленном решении не может быть верных разрядов, поскольку истинное решение не существует. Вследствие округлений в ходе гауссова исключения будут получены невы- рожденные, но очень плохо обусловленные матрицы 5 и У.

В этом примере при вычислениях в МаПаЬ'е с 1ЕЕЕ-арифметикой двойной точности вычисленная невязка т из-за округлений оказывается тпочнмм нулем, поэтому обе оценки 69 2.4. Анализ ошибок (2.13) и (2.14) обращаются в нуль. Исправление оценки (2.13) путем добавления 4еИА(! )~хб + ЗЬ|О дает значение, большее 1, как и требовалось. К сожалению, вторая, «более точная» оценка (2.14) дает приблизительно 10 ~, ошибочно указывая, что семь разрядов вычисленного решения верны. Вот как строится этот пример.

Положим зг = 3/2зз, ~ = 2", 9 1553 10-з 1 6384 104 1 6384 104 6.1035 10 з 6.1035 10 з 0 6.1035 10 з — 3.4106 . 10 'з 6.1035 10 з и Ь = А . [1, 1+ е, 1]т. Матрица А может быть вычислена без каких либо округлений, но при вычислении Ь округления происходят, что означает: Ь лишь приближенно принадлежит линейной оболочке столбцов матрицы А, поэтому система Ах = Ь не имеет решений. Проводя гауссово исключение, получим 1 0 0 7, .66666 1 0 .66666 1.0000 1 2 9 1553, 10-з 1 6384, 104 1 6384 10 17 0 1.0923.

104 — 1 0923 104 0 0 1.8190.10 зз что дает для обратной матрицы значение 2 0480, 10з 5 4976, 10ы 5 4976, 10 А-~ 2 0480, 10з 5 4976, 10м 5 4976, 10 20480 10з 54976 10ы 54976 10м Отсюда следует, что в вычисленной матрице )А '( . )А! все элементы приблизительно равны числу 6.7109 10~, поэтому вычисленное кап(А) есть величина порядка 0(10"). Другими словами, наша оценка для ошибки указывает, что в вычисленном решении имеется 16 — 7 = 9 верных разрядов, тогда как, в действительности,ни один не верен. Если исключить сильный рост элементов, то можно показать, что феномен, иллюстрируемый данным примером, не может сделать слишком малой оценку (2.13) (где число Йгб надлежащим образом увеличено). Кахвн нашел еще семейство вырожденных и х п-матриц, для которых замена одного очень малого элемента (равного — 2 ") нулем понижает значение кся(А) до величины порядка 0(пз).

Аналогичным образом можно было бы построить примеры, где А не является точно вырожденной матрицей, так что оценки (2.13) и (2.14) верны в точной арифметике, но, вычисленные в машинной арифметике, становятся слишком малы. О 70 Глава 2. Решение линейных уравнений 2.5.

Улучшение точности приближенного решении Мы только что видели, что ошибка в вычисленном решении системы Ах = Ь может достигать уровня к(А)е. Что делать, если эта ошибка слишком велика? Одна из возможностей состоит в том, чтобы провести все вычисление повторно, используя арифметику более высокой разрядности, однако это может потребовать больших затрат и времени, и памяти. К счастью, при не слишком большом к(А) имеются гораздо более экономные методы получения более точного решения.

При решении любого уравнения )'(х) = 0 можно попытаться использовать метод Ньютона, чтобы улучшить приближенное решение х; посредством формулы х>>> = х, — фф.. Применяя эту формулу к функции Дх) = Ах — Ь, получаем один шаг итерационного утпо гнения: т = Ах, — Ь решить систему Ад = т относительно вектора д хг~ь> = х, — д. Если бы мы могли точно вычислить вектор т = Ах, — Ь и решить систему Ад = т точно, то все закончилось бы в один шаг, чего и следует ожидать от метода Ньютона, применяемого к линейной задаче. Ошибки округлений препятствуют этой моментальной сходимости.

Алгоритм представляет теоретический и практический интерес именно в ситуации, когда матрица А настолько плохо обусловлена, что можно получить лишь очень неточное решение системы А4 = т (как и системы Ахо = Ь). Теорема 2.7. Предполохсим, чтпо вектор т вычисляется в арифметике двойной точности, и пусть к(А) е < с = з„4 — — < 1, где п — порядок матрицы А, а д — коэ4>рициент роста. Тогда, повторяя шаги итерационного уточпенил, получим вектор х„длл которого 'Охг — А >Ь|! !)А >ЬЙ Заметпим, чтпо в этой оценке ошибки не присутствует число обусловленности. Это означает, что решение системы можно втячислить с высокой точностью независимо от значения числа обусловленностпи при условии, что произведение к(А)е значитпельно меньше, чем 1.

(На т>рантике, с является слишком консерватпивной верхней границей и алгоритпм часто работаетп успешно, даже если к(А)е ) с.) Набросок доказательства. Чтобы доказательство оставалось прозрачным, будем учитывать лишь наиболее важные ошибки округлений. Для краткости, вместо 0 . '0 будем писать )! .

'й. Наша цель — установить неравенство Йх,+> — хЙ < ))хг — х)( = Цх, — хЙ. к(А)е с По предположению, ~ < 1, поэтому из неравенства следует, что ошибка ~~хт ь>— х~( монотонно убывает до нуля. (На практике, ее убывание в какой-то момент прекратится вследствие ошибок округления при присваиваниях х,в> — — х; — д„ которые в нашем анализе игнорируются.) 71 2.5.

Улучшение точности приближенного решения Оценим сначала ошибку в вычисленной невязке т. Имеем т = 6(Ах; — Ь) = Ах; — Ь+ 1, где, согласно результату из вопроса 1.10, )Д < цвз()А! ~х,~ + (ЬО + в)Ах, — Ь! тн е(Ах, — Ь|. Член порядка е~ возникает при вычислении т с двойной точностью, а член порядка е — при округлении вычисленного т до числа с обычной точностью.

Поскольку ез « в, то мы пренебрежем членом порядка вз в оценке для ~Л. Далее, имеем (А + бА)с1 = т, где, вследствие оценки (2.11), ()бАЙ < ув6А6 и 7 = Зпзд, хотя зта оценка обычно слишком груба. Как уже упоминалось, мы упрощаем анализ предположением, что присваивание хьы = х; — и' выполняется точно. Продолжая игнорировать члены порядка в~, получаем и' = (А+ бА) 'т = (1+ А 'бА) 'А 'т = (1+А 'бА) 'А '(Ах; — Ь+ 1) = (1+А 'бА) '(х; — х+А '1) (1 — А 'бА)(х; — х+ А '1) х; — х — А 'бА(х, — х) + А "1. Следовательно, х;~~ — х = х; — Н вЂ” х = А ~бА(х; — х) — А "1', а потому !)х;~.~ — хЙ < ()А ~бА(х; — х)Й + ((А ~Д < !)А ' )) !)бА)! .

)(х, — х!) + ((А ')! в '6 Ах, — Ь|) = '6А 'и' ЙбА(! йх; — х()+0А ~)!.в. ЙА(х, — х)!) < йА ')! .ув )(Ай '6х, — х6 +((А '6 йА6 е.))х; — х)! = ))А ')! ЙА() в (у+1) !)х, — хй. Таким образом,из неравенства ~ = )/А '!) . /(А!! в(у+ 1) = к(А)в/с < 1 вытекает сходимость итераций. П Итерационное уточнение или иные варианты метода Ньютона могут использоваться для улучшения приближенных решений и во многих других задачах линейной алгебры. 2.5.1. Итерационное уточнение в арифметике обычной точности Этот раздел может быть опущен при первом чтении книги. Иногда для проведения итерационного уточнения отсутствует предусматриваемая им возможность выполнять вычисления с двойной точностью. Если, например, уже входные данные суть числа с двойной точностью, то нам потребовалось бы вычислять невязку с учетверенной точностью, что не всегда возможно. На некоторых машинах (скажем, РС с процессором 1псе1 Репгппп) доступна расширенная двойная точность, прибавляющая к обычной двойной точности 11 дополнительных двоичных разрядов (см.

рвзд. 1.5). Это меньше, чем нужно для учетверенной точности (которая потребовала бы, по меньшей мере, 2 ЗЗ = 106 двоичных разрядов), но все же ощутимо повышает основную точность. Глава 2. Решение линейных уравнений Однако если ни одна из этих возможностей не доступна, итерационное уточнение все-таки можно использовать, вычисляя невязку г с обычной точностью (т. е. той точностью, какую имеют входные данные). В этом случае теорема 2.7 теряет силу. С другой стороны, теорема, приводимая ниже, показывает, что при некоторых технических предположениях имеет смысл выполнить один шаг итерационного уточнения в обычной точности, поскольку это уменьшает по- компонентную относительную обратную ошибку, определенную в теореме 2.3, до уровня 0(е).

Если соответствующее относительное число обусловленности ксл(А) = ]~]А ~~ ]А]~] (см. равд. 2.2.1) значительно меньше обычного числа обусловленности к(А) = ]]А ']) ° ])А]), то полученное решение будет точней исходного. Теорема 2.8. Предположим, что вектор г вычисляется с обычной точно- стью, и пусть ])А ')] .)]А]) ., ' '.е(1. ш1пг()А! ]х]); Тогда один шаг итерационного уточнезгия дает вектор хд, такой, что (А + бА)хз = Ь+ бЬ, где ]бай] = 0(е)~а,.] и ]бЬ;] = 0(е)]Ь,].

Ин ми словами, по- компонентная относительная обратная ошибка для х1 имеет мип мальпый возможный порядок малости. Это означает, например, что при разреженных А и 6 возмущения БА и 66 имеют ту же структуру разреженности, что и, соответственно, А и Ь.

Доказательство теоремы можно найти в [149]; относительно других деталей см. ]14, 225, 226, 227]. Итерационное уточнение с обычной точностью и оценка ошибки (2.14) реализованы в ЬАРАСК'е, например, программой вйевчх. Пример 2.7. Рассмотрим те же матрицы, что в примере 2.5, и проделаем один шаг итерационного уточнения в той же точности, с какой выполнялись остальные вычисления (е 10 'в). Для этих матриц обычное число обусловленности к(А) 10'4, тогда как кои(А) = 1, поэтому нужно ожидать большого улучшения точности. В самом деле, покомпонентная относительная ошибка метода СЕРР и соответствующая ошибка из формулы (2.14) становятся меньше 10 'г. Ма$1аЬ-программу для этого примера можно найти по адресу НОМЕРАСЕ/Мас1аЬ/р1чоаш.

О 2.5.2. Уравновешиванне Имеется еще один общий прием для повышения точности при решении линейных систем, называемый уравновешиванием. Под ним понимают выбор подходящей диагональной матрицы Р с последующим решением системы РАх = .РЬ вместо Ах = Ь. Матрипу выбирают с таким расчетом, чтобы для РА число обусловленности было меньше, чем для А (если это возможно). Так, в примере 2.7, беря в качестве ди число. обратное 2-норме 1-й строки матрицы А, получим матрипу РА, очень близкую к единичной; тем самым число обусловленности снижается с 10ы до 1. Можно показать, что подобный выбор матрицы Р всегда уменьшает число обусловленности произведения РА до уровня, не более чем в ,/п раз превышающего минимум по всем возможным выборам диагональной 2.6.