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

Эти три случая иллюстрируются рисунком. С помощью его правой части легко понять, почему число обусловленности бесконечно при 0 = х/2: в этом случае имеем х = 0 и почти любое, как угодно малое изменение в А или Ь приведет к ненулевому решению х, что является «бесконечно» больп(им относительным изменением. к ья = 0(к(А( ( к ьб =бес«сне«нос«» Оценке из теоремы 3.4 можно придать другую форму, в которой отсутству ет член с 0(е~) [258, 149): ||х — х||г екг/А) (г ||г||г '1 ||г — г|!г Здесь т — невязка возмущенной задачи: г = 1А + ЬА)х — (Ь + ЬЬ).

Мы увидим, что при правильной реализации оба метода — и 12й, и ЯУП,— численно устойчивы, т. е. дают приближенное решение х, которое 1точно) минимизирует некоторую функцию вида ||(А+ ЬА)х — 1Ь+ ЬЬ)||г, причем /||ЬА|| ||БЬ|!»( (( ||А|| ' ||Ь!|! Комбинируя этот результат с приведенными выше оценками возмущений, можем получить оценки для ошибки в решении задачи наименьп(их квадратов подобно тому, как были получены оценки для ошибки в решении системы линейных уравнений. Метод нормальных уравнений не столь точен.

Поскольку репгается система (АтА)х = АгЬ, точность определяется числом обусловленности кг(АтА) = кгг(А). Таким образом, ошибка всегда ограничена величиной типа кг|А)е, а не просто кг(А)е. Следует ожидать поэтому, что при решении нормальных уравнений может быть потеряно вдвое больше верных разрядов, чем в методах, основанных на Щ-разложении и ЯЪ'11. Кроме того, метод нормальных уравнений не всегда устойчив, т.е. вычисленное решение х, вообще говоря, не является точкой минимума функции вида |||А+ ЬА)х — 1Ь+ ЬЬ)||г с малыми ЬА и бЬ. Все же, если число обусловленности мало, метод, скорей всего, даст приближенное решение примерно того же качества, что и с)11-разложение и ЯН0. Поскольку задача наименьших квадратов быстрей всего решается именно посредством нормальных уравнений, этот метод является предпочтительным в случае хорошо обусловленной матрицы А.

В разделе 3.5 мы вернемся к вопросу о том, как следует решать очень плохо обусловленные задачи наименьших квадратов. 129 3.4. Ортогонааьные матрицы 3.4. Ортогональные матрицы Как было отмечено в разд. 3.2.2, процесс ортогонализации Грэма — Шмидта (алгоритм 3.1), примененный к почти линейно зависимым векторам, может дать матрицу Я, столбцы которой не ортогональны; таким образом, с его помощью нельзя вычислить ЯВ;разложение устойчиво. Поэтому наши алгоритмы будут основаны на использовании некоторых легко вычисляемых ортогональных матриц, нззываемьгх отражениями Хаусхолдера и враи1ениами Гивенса Подходящим выбором таких матриц можно получать нулевые компоненты в их векторных сомножителях.

Позднее мы покажем, что всякий метод, использующий эти матрицы для исключения элементов в матрице задачи, автоматически является устойчивым. Этот анализ ошибок будет приложйм к нашим алгоритмам ЯГс-разложения, а также ко многим алгоритмам вычисления собственных значений и ЯЧП в гл. 4 и 5. Несмотря на отмеченную возможность вычисления матрицы С1 с неортогональными столбцами, алгоритм МОЯ имеет важные применения в численной линейной алгебре. (Напротив, его менее устойчивая версия ССБ малополезна.) Среди них вычисление собственных векторов симметричных трехдиагональных матриц с использованием бисекции и обратной итерации (разд.

5.3.4), а также алгоритмы Арнольди и Ланцоша для приведения матрицы к некоторым «конденсированным» формам (рэзд. 6.6.1, 6.6.6 и 7.4). Алгоритмы Арнольди и Ланцоша используются в качестве основы алгоритмов для решения разреженных линейньгх систем и вычисления собственных значений разреженных матриц. Алгоритм МОЯ можно модифицировать таким образом, чтобы задача наименьших квадратов решалась устойчиво, однако столбцы матрицы Я по-прежнему могут сильно отклоняться от ортогональности (33).

3.4.1. Преобразования Хаусхолдера Преобразованием Хаусхолдера (или отражением) называется матрица вида Р = 1 — 2иит, где ЙЩ = 1. Легко видеть, что Р = Рт и РРт = (1 — 2ии«)(1— 2иит ) = 1 — 4иит + 4иитиит = 1, т. е. матрица Р симметрична и ортогональна. Она называется отражением, потому что вектор Рх является отражением вектора х относительно плоскости, проходящей через О перпендикулярно к и. Пусть дан вектор х. Тогда легко найти отражение Р = 1 — 2иит, аннулирующее в векторе х все компоненты, кроме первой: Рх = (с, О,...,0)т = с еы Это можно сделать следующим образом.

Имеем Рх = х — 2и(итх) = с еы поэтому и = -+-1(х — се|), т. е. и есть линейная комбинация векторов х и еы Так как ~)х~~т = !)Рх()з = )с), то и должен быть параллелен вектору й = х х 'ОхОтеы 130 Глава 3. Линейные задачи наименыиих квадратов откуда и = йДй~)з. Можно проверить, что при любом выборе знака получится вектор и, удовлетворяющий соотношению Рх = ег, если б ф О.

Мы будем поль- зоваться формулой =х + вяп(хг)ег, так как в этом случае не будет взаимного сокращения при вычислении первой компоненты в й. Итак, имеем хг + вбп(хг) !)х)~з хз й и и=— 1ЯЬ' й= Мы будем записывать это как и = Новее (х). (На практике можно хранить 5 вместо и, чтобы сэкономить на работе вычисления вектора и; в этом случае вместо Р = 1 — 2иит используется формула Р = 1 — (2ДйЯ)йб~.) Пример 3.5. Покажем, как с помощью отражений вычисляется цГс-разложение матрицы А размера 5 х 4. Из этого примера будет ясно, как следует действовать для произвольных т х п-матриц.

В приводимьгх ниже выкладках Р; суть ортогональные 5 х 5-матрицы, х означает произвольный ненулевой элемент, а символ о используется для нулевых элементов. хххх оххх оххх оххх А, =РА= 1. Выбрать матрицу Рг так, чтобы оххх х ххх оххх оохх о охх 2. Выбрать матрипу Рз = ~0 Р,~( так, чтобы Г1 01 Аз = РзАг —— о охх хххх оххх о охх о о ох Р 01 3. Выбрать матрипу Рз — — ~ 1 ~ так, чтобы Рз Аз = РзАз = о о ох хххх оххх 4. Выбрать матрицу Р4 = так, чтобы А4 — = РаАз = оохх о о ох оооо Здесь отражение Р,' выбиралось таким образом, чтобы аннулировать поддиагональные элементы столбца 1. Это не изменяет нулей, уже введенных в предыдущие столбцы. Обозначим полученную верхнетреугольную 5 х 4-матрицу через В = А4.

Тогда А = Рг Рз Рз Рс В = ЯВ, где Я образована первыми четырьмя столбцами матрицы РгтРзтРзтРс — — РгРзРзР4 (так как все матрицы Р; симметричны), а В состоит из первых четырех строк матрицы В. О Приведем общий алгоритм ЯК-разложения, основанный на использовании отражений. 131 ЗА. Ортогональные матрицы Алгоритм 3.2. ЯЯ-разложение посредством отражений: 7ог 4 = 1 Со шш(т — 1,п) и; = Нопяе (А(с': т,1)) Р! = 1 — 2и;ир А(г: т, 1: и) = Р,"А(1: т, г': п) епд 7ог Обсудим некоторые детали реализации метода.

Матрицы Р, нигде не требуется формировать в явном виде, так как умножение на них производится в соответствии с более экономичной формулой (Х вЂ” 2и,и~)А(г'; т,д: п) = А(1: т,г: п) — 2и,(и~А(1: т,1: п)). Для хранения матрицы Р; достаточно запомнить лишь вектор и, нли же й, н число ~(йД. Эта информация может храниться в столбце 1 матрицы А; в действительности, в этом столбце ничего не нужно менять! Таким образом, ЯК-разложение может быть записано на место матрицы А, причем Я хранится в факторизованной форме Р,... Р„ы а Р, хранится в виде вектора й; в подлиагональных позициях столбца 1 матрицы А. (Поскольку диагональные позиции заняты элементами ги, нужен дополнительный массив длины п для хранения первых элементов векторов й,.) Вспомним, что для решения задачи наименьших квадратов ш1пЙАх — ЬЙя с помощью ЯК-разложения А = ЧК нужно вычислить вектор ЯтЬ.

Это делается следующим образом: ЯгЬ = ЄР~...РяЬ, поэтому Ь нужно последовательно умножать на Р„Ря,... Рьс 1ог 1 = 1 Со п 7= — 2 и,'Ь(1: Ь(4: т) = Ь11: т) + 7и; енса Сог Стоимость этого алгоритма — п скалярных произведений 7 = — 2 и~Ь и п операций типа яахру Ь+ 7и.. Стоимость вычисления разложения А = ЯК этим способом составляет 2пэт — яэпз флопов; цена последующего реянения задачи наименьших квадратов при известном С4К-разложении равна уже только 0(тп) флопов. В ЬАРАСК'е решение задачи наименьших квадратов методом ЯК-разложения осуществляется программой яяе1я. ЯК-разложение можно реорганизовать так, чтобы в нем использовались матричные умножения н другие операции уровня 3 ВЬАБ, подобно тому, как это было сделано для гауссова исключения я разд.

2.6; см, по этому поводу вопрос 3.17. Если в МаС1аЬ'е т х п-матрица А имеет больше строк, чем столбцов, а Ь имеет размер т х 1, то оператор ЬФ решает задачу наименыпих квадратов. Если требуется сяК-разложение само по себе, то оно вычисляется оператором [Я,К]=с1г(А). 3.4.2. Вращения Гнвенса ~ сояд — сбпй ) Вращение Гивенса (или просто вращение) К(о) ив я ~ . о о ~ поворачи- ~ яйпй сояэ' вает всякий вектор з й КЯ на угол о против часовой стрелки (см. рисунок). 132 Глава 3.

Линейные задачи наименьших квадратов (Е)и Нам потребуется также определить вращение на угол В в координатной плоскости (ю,у) пространства К": соя 0 — яйп 0 В(г',у,В) = соя 0 я1п В При заданных г и у и векторе х можно аннулировать компоненту х, выбирая соя В и яш В таким образом, чтобы т. е. соя 0 = — *- — и я! и 0 = — — *-,к —. / 2е т / 2~ та' Алгоритм ЯК-разложения, основанный на вращениях, аналогичен алгоритму, использующему отражения, однако при обработке столбца г элементы в нем аннулируются по одному за шаг (в направлении, скажем, снизу вверх).

Пример 3.6. Проиллюстрируем два промежуточных шага в вычислении ЯК- разложения матрицы 5 х 4 посредством вращений. Чтобы перейти от х х х х х х х х о х х х о о х х о о о х о о о х о о х х производим умножения х х х х о х х х х х х х о х х х с — я я с о х х х о о х х о о х х о о х х о о х х о о х х о о х х о о х х о о о х 133 ЗА. Ортогонэльныв матрицы х х х х х х х х о х х х 1 с — г ! г' с' 1 о х х х о о х х о о х х о о х х о о о х о о о х о о о х 3.4.3.