Дж. Деммель - Вычислительная линейная алгебра (1156793), страница 30
Текст из файла (страница 30)
Анализ ошибок для ортогональпых матриц Этот анализ устанавливает обратную устойчивость (ггь-разложения и многих алгоритмов вычисления собственных значений и сингулярных чисел, которые обсуждаются в дальнейшем. Лемма 3.1. Пусть Р— точное отражение (или вращение), а Р— его приближение в арифметике с плавающей точкой.
Тогда Я(РА) = Р(А+ Е) Щ/г = 0(г) ЦА)/г Я(АР) = (А+Р)Р, Щ!г = 0(в) ЦА)!г. Набросок доказательства . Применим обычное правило Я(аОЬ) = (аОЬ)(1+в) к формулам для вычисления Р и ее умножения на матрицу А. См. по этому поводу вопрос 3.16. П В словесном описании, лемма утверждает, что умножение на отдельную ортогональную матрицу является обратно устойчивым процессом. Стоимость ()В;разложения с использованием вращений вдвое превышает стоимость метода отражений. В дальнейшем вращения понадобятся нам и для других целей. Остановимся на некоторых деталях реализации. Как и в методе отражений, мы можем записать сомножнтели С4Ть-разложения на место матрицы А.
Используется тот же прием, т. е. информация, описывающая преобразования, помещается на места аннулируемых элементов. Поскольку вращение аннулирует лишь один элемент, для хранения информации о нем имеется только одно машинное слово. Оно используется следующим образом: пусть г = эш0 и с = сову. Если ~г! < ~с(, то в память записывается число г э8п(с), в противном случае, число -'-~-"~'-~. Обозначим хранимое значение через р. Числа г и с восстанавливаются из него так: если ~р~ < 1, то полагаем г = р и с = ~/1 — гг; в противном случае, с = — ' и г = ~/à — сг.
Если всегда запоминать значение г, а затем вычислять с = ~/Т вЂ” гг, то при г, близком к 1, значение с будет получено с большой погрешностью. Отметим еще, что вместо исходной пары с, г мы можем восстановить пару — с, — г; однако для наших целей обе пары эквивалентны. Имеется возможность выполнения последовательности вращений с затратой меньшего числа флопов, чем в описанном вьппе процессе. Она состоит в использовании так называемых бьгсгпрых вращений [7, 8, 33]. Поскольку при вычислении (4В-разложения быстрые вращения все же медленней отражений, мы не станем их рассматривать здесь. Глава 3.
Линейные задачи иаимеиыпих квадратов Теорема 3.5. Пусть к матрице Ао применлетсл последовательность ортогональныт преобразований. Тогда вычисленное произведение есть точное ортогональное преобразование матрицьг Ао + 6А, где !!бА!!г = 0(е)!!А!!г. Иначе говоря, процесс в целом обратно устойчив: 6(Р1Р1-г ° РгАоЯгЮг . Яз) = Рг. ° .Рг(Ао+ЕЮг .. ° Яз, причем !!Е!!г = 1. 0(е) !!А!!г. Здесь, как и в лемме 3.1, Р; и ь„ч суть ортогональные матрицы, вычисленные в арифметике с плавающей точкой, а Р, и (,>г суть точньге ортогональные матрицы.
Док зательство . Положим Р.— : Р ...Рг и Я = Яг...ф. Мы хотим пока- зать, что А = й(Р1А1 Яу) = Р, (А+Ест для некоторой матрицы Егч такой, что !!Ез!!г — — 10(е)!!А!!г. Будем рекурсивно пользоваться леммой 3.1. При 1 = 0 доказывать нечего. Предположим теперь, что утверждение теоремы верно для индекса 1 — 1. Тогда имеем В = б(РгА, = Р,(А,, +Е') по лемме 3.1 = Рз(Р1 г(А+Ез г)Я1, +Е') по предположению индукции = Р.(А+Е., + Р Е'Я,Я1 = Рз(А+ ЕиД, где !!Еи!!г = !!Ез д + Рг 1Е Я.
г!!г < !!Ез г!!г + !!Рг дЕ Яд г!!г = !!Ез — г !!г + !!Е !!г = 10(е)!!А!!г, поскольку !!Ег г!!г = (1 — 1)0(е)!!А!!г и !!Е'!!г = О(е)!)А!!г. Таким же образом анализируется правое умножение на матрицу Я,. П 3.4.4. Почему именно ортогональные матрицы3 Посмотрим, как вела бы себя ошибка, если бы матрицу Ас в теореме 3.5 вместо ортогональных матриц мы должны были бы умножать на последовательность неортогональных матриц.
Пусть Х вЂ” точная неортогонвльная матрица, а Х— ее приближение в арифметике с плавающей точкой. Тогда, применяя обычный анализ округлений матричного умножения, имеем й(ХА) = ХА+ Е = Х(А+ Х 'Е) = Х(А+ Р), где !!Е!!г < 0(е)!!Х!!г !!А!!г. Поэтому !!Р!!г < !!Х '!!г !!Е!!г < 0(е) кг(Х).Ц!АЦ!г. Итак, ошибка !!Е!!г приобретает множитель, равный числу обусловленности кг(Х) > 1. В более длинном произведении Хь... Хд АУю .. Уь ошибка умножилась бы на П,кг(Х;) кг(У,).
Этот множитель минимален (и принимает значение 1) тогда и только тогда, когда все матрицы Х; и Уг ортогонвльны (или являются скалярными кратными ортогонвльных матриц). 3.5. Задачи наименьших квадратов неполного ранга 135 3.5. Задачи наименьших квадратов неполного ранга До сих пор при минимизации функции ||Ах — Ь||г предполагалось, что А — матрица полного ранга. Что произойдет, если ранг матрицы не полон или же А «близка» к матрице неполного ранга? Такие задачи возникают во многих реальных ситуациях, например при выделении сигнала из зашумленных данных, решении некоторых типов интегральных уравнений, цифровом восстановлении изображений вычислении, обратного преобразования Лапласа, и т.
д. [141, 142|. Эти задачи очень плохо обусловлены и чтобы превратить их в хорошо обусловленные, приходится накладывать дополнительные ограничения на их решения. Превращение плохо обусловленной задачи в хорошо обусловленную путем наложения дополнительных условий на решение называется регуллризацией. Эта техника используется и в других областях численного анализа при возникновении плохо обусловленных задач.
Например, следующее предложение показывает, что для матрицы А, имеющей (в точной арифметике) неполный ранг, решение задачи наименьших квадратов не единственно. Предложение 3.1. Пусть А — матрица размера т х п, причем т > п, и пусть ранг т матрицы А меньше п. Тогда множество векторов х, минимизирующих ||Ах — Ь!|г, образует (п — г)-мерное линейное надпространство.
Доказательство. Пусть Аг = О. Если х минимизирует ||Ах — Ь!|г, то это же верно для вектора х + и. П Вследствие ошибок в элементах матрицы А или округлений при вычислениях, у А, как правило, не будет сингулярных чисел, в точности равных нулю, но будут одно или несколько очень малых сингулярных чисел. Следующее предложение показывает, что единственное решение такой задачи, скорей всего, очень велико и обязательно является очень чувствительным к погрешностям в правой части Ь (см. также теорему 3.4). Предложение 3.2. Пусть о' ш = о„а„(А) есть наименьшее сингулярное чис- ло матрицы А. Предположим, что о;„> О.
Тогда: 1. если х минимизирует ||Ах — Ь||г, то ||х||г > |ит Ь|/и я„, где и„— последний столбец матрицы У в разложении А = УКЪ'т; 2, замена Ь на Ь+ 6Ь приводит к замене х на х+ бх, где ||бх||а может достигать величине«||6Ь||а/а Иными словами, если А близка к матрице неполного ранга /т. е. число в ы мало), то решение х плохо обусловлено и, возможно, очень велико. Доказательство.
Для доказательства первого утверждения запишем х А«Ь = Ь'Х»П~Ь, откуда ||х||г = ||Е ~От ЬЦг > |(Е ~П~Ь)„| = |итЬ!/оаа . Чтобы доказать утверждение 2, выберем 6Ь, параллельное вектору и„. П Обсуждение регуляризации начнем указанием способа регуляризовать задачу наименьших квадратов, которая имеет неполный ранг в точной арифметике. Пусть А — матрица размера т х п и ранга г ( и. В (и — т)-мерном подпространстве решений задачи наименьших квадратов будем искать единственное решение с минимальной нормой (нормальное решение).
Это решение охарактеризовано в следующем предложении. 136 Глава 3, Линейные задачи наименьших квадратов Предложение 3.3. Для (точно) втярожденной матрицы А векторы х, минимизируютцие ||Ах — Ц|г, могут быть охарактеризованьг так: пусть А имеет ранг т < п и $ЪР А = ЬтЕ)тт. Запишем Я1тР в виде '= |О'~та| ~ О О ) |" "| = ~тЕ" | Ет О1 т т где Ет — невырожденная тхг-матрица, а каждая из матриц Ут и Ъ'т состпоит из г столбцов.
Обозначим наименьшее ненулевое сингулярное число матрицы А через о = опиь(Ет). Тогда: 1. все решения х описьгваются формулой х = )ттЕ ~Утз Ь+ Ъгг, где г — произвольный вектор; 2. решение х тогда и только тогда имеет минимальную В-норму, когда г = О. В этом случае х = 1ттЕ ~От Ь и ||х||г < ||Ь||г/о; 3. замена Ь на Ь+ бЬ может изменить нормальное решение х на величину, норма которой не превосходит ||бЬ||г/о. Иначе говоря, норма и число обусловленности единстивенного решения х с минимальной нормой (нормального решения) определяются наименьшим ненулевым сингулярным числом матрицы А.
Доказательство . Выберем У таким образом, чтобы т х т-матрица |У, Ьт] = [Уы Уг, Ц была ортогональной. Тогда ||Ах — Ь||г г= |||Ьт У)т(Ах — Ь)||гг Угу (У1Ет т х Ь) ттт ||Х 1ттх — тттЦг + ||тттЦг + ||тттЬ||э 1. Величина ||Ах — Ь(|г минимальна, когда ХдЪдтх = УттЬ, или х = 1ттХ тата Ь+ 1тгг, так как )тт~1тгг = О для любого вектора г. 2. Столбцы матриц 1тт и 1тг взаимно ортогональны.
Поэтому, по теореме ПифагоРа, ||хЯ = ||1ттЕ, 'Уттб||э + ||)таз|Я. Эта сУмма минимальна пРи г = О. 3. Замена Ь на бЬ изменяет х на величину, оцениваемую как Щ Х, 'УттбЬ||э < ||Х, '||г||бЬ!|г = ||бЬ!|г/о. 0 Смысл предложения 3.3 в том, что нормальное решение х единствешю и может быть хорошо обусловленным, если наименьшее ненулевое сингулярное число пе слишком мало. Это дает ключ к практическому алгоритму, обсуждаемому в следующем разделе. Пример 3.7.
Предположим, что проводится медицинское исследование воздействия некоторого препарата иа уровень сахара в крови. Для каждого пациента (им присвоены номера от т = 1 до т) регистрируются начзльиый (а, т) и 137 3.5. Задачи наименьших квадратов неполного ранга конечный (Ь,) уровень сахара в крови, количество введенного препарата (а, г) и прочая медицинская информация, в том числе вес в каждый из дней недельного курса лечения (а, г,..., асг).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
