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

В целом для каждого пациента измеряются л ( т медицинских величин. Целью при этом является прогноз значения Ь, ло величинам а; ы..., а, „, что рассматривается как задача наименьших квадратов ппп лАх — Ьлг. Планируется использовать решение х для предсказания конечного уровня сахара Ьд в крови будущего пациента у посредством вычисления скалярного произведения ~ г, аггхг.

Поскольку обычно вес человека незначительно меняется изо дня в день, столбцы матрицы А с 3-го по 9-й, т.е. те, что содержат веса, скорей всего, будут очень похожи. Для целей нашего рассуждения примем, что столбцы 3 и 4 одинаковы (так действительно может быть, если значения весов округлены до ближайшего целого числа фунтов). Это означает, что ранг матрицы А не полон и вектор ха = [О, О, 1, — 1, О,..., 01~ принадлежит ее ядру. Поэтому, если х есть (нормальное) решение задачи наименьших квадратов пппгОАх — Ьйг, то х+ дхо также является решением (не обязательно нормальным) для любого числа 3, например, 3 = 0 и Д = 10г.

Есть ли какие-нибудь основания для того, чтобы предпочесть одно значение Д другому? Ясно, что выбор 3 = 10г неудачен: если будущий пациент у приобретет фунт веса с первого дня на второй, то в предсказание конечного уровня сахара в крови, т.е. в величину ~ ", а,гхг, эта разница в один фунт войдет умноженной на 10г. Выбор 3 = О, что соответствует нормальному решению х, гораздо более разумен. О Дополнительные аргументы, обосновывающие использование нормальных решений в задачах неполного ранга, можно найти в (141, 142). Если А — квадратная невырожденная матрица, то единственным решением задачи Ах = Ь является, конечно, вектор х = А 'Ь.

Если в А строк больше, чем столбцов, а ранг, возможно, не полон, то единственное нормальное решение задачи наименьших квадратов может быть записано в сходном виде х = А+Ь. Участвующая здесь псевдввбратиал матрица Мура — Пгнрврза А+ определяется следующим образом: Определение 3.2. (Псевдообратная Мура — Пенроуза А+ для матрицы А, возможно, неполного ранга.) Положим А = ПЕ$'т — ПгЕгЪ'т, как в равенстве (3.1).

Тогда Аэ УгЕ, 'Пгт. Это можно записать и в виде А+ = ЬгЕ+Пг, гдг Е+ Е~ 0 Е,' 0 Таким образом, решение задачи наименьших квадратов всегда дается формулой х = А+Ь. Если А — матрица неполного ранга, то такой вектор х имеет минюгальную норму среди всех решений задачи. 3.5.1. Решение задач наименьших квадратов неполного ранга посредством Б'г'гг Наша цель состоит в том, чтобы вычислить нормальное решение х в условиях приближенных вычислений. В предыдущем разделе мы видели, что нормальное решение единственно и его число обусловленности зависит от наименьше- 138 Глава 3. Линейные задачи наименьших квадратов го ненулевого сингулярного числа. Поэтому вычисление нормального решения предполагает знание наименьшего ненулевого сингулярного числа, а потому и ранга матрицы А.

Основная трудность при этом состоит в том, что ранг является разрывной функцией от матрицы. Рассмотрим, например, вырожденную 2 х 2-матрицу А = бйа8 (1, 0). Ее наименыпее ненулевое сингулярное число есть о = 1. Согласно предложению З.З, нормальным решением задачи ш1п,'ОАх — Ь!)т с 5 = (1,1)т является вектор х = [1, О)т, а соответствующее число обусловленности равно 1/о = 1. Если посредством произвольно малого возмущения перейти к матрице А = е(1а8 (1, е), то о уменьшается до е, а вектор х = (1,1/е) становится огромным, квк и число обусловленности 1/е.

Как правило, такие малые возмущения величины 0(е)йА()з н происходят при округлениях. Как мы только что видели, они способны увеличить число обусловленности с 1/о до 1/е. Алгоритмически эта разрывность в поведении учитывается следующим образом: в общем случае, всякое вычисленное сингулярное число У, удовлетворяет оценке )о; — о;! < 0(е)((А5т. Это является следствием обратной устойчивости, а именно вычисленное сингулярное разложение есть точное БЧП для слабо возмущенной матрицы А = ОЕУт = А+ 5А, причем йоА!) = 0(е) .

~)Ай. (Более детальное обсуждение вопроса дано в гл. 5.) Стало быть, всякое д;, удовлетворяющее неравенству а; < 0(е) 8А~(з, может рассматриваться как нуль, потому что д; не отличимо от нуля в пределах точности вычислений. Применительно к нашему 2 х 2-примеру, это означает, что мы должны были бы заменить е в матрице А нулем до решения задачи наименьших квадратов. Это увеличило бы наименьшее ненулевое сингулярное число с е до 1; соответственно, число обусловленности уменьшилось бы с 1/е до 1/о = 1. Более общо, будем считать, что пользователем задана мера Со1 неопределенности в элементах матрицы А.

Из наличия округлений следует, что Фо1 > е '5 А 5, но значение 1о1, в зависимости от происхождения А, может быть и много большим, чем эта граница. Положим теперь ее; = дм если д; > со1, и д, = О, в противном случае. Пусть Е = сйа8(д,). Произведение 1) Ю~ назовем усеченным сингулярным разложением матрицы А, поскольку сингулярные числа меньшие, чем 1о1, заменены нулями. Решим задачу наименьших квадратов, пользуясь усеченным ЯЧО вместо исходного. Обоснованием этого служит оценка ()УЕЪ УЕЪ ит — яУ(Е Е)Ъ из < ьо1, г.

е. изменения в А, вызванные заменой о; на Ю,, меньше, чем изначально присутствующая неопределенность в данных. Основным доводом в пользу Е является то, что среди всех матриц, находящихся от Е на расстоянии, не превышающем Фо1, Е максимизирует наименьшее ненулевое сингулярное число о. Иначе говоря, такой выбор минимизирует величину нормального решения х и его число обусловленности. На рисунке проиллюстрированы геометрические взаимосвязи между исходной матрицей А, матрицей А = УЕк'т и матрицей А = УЕк'т. Каждая из матриц изображается точкой евклидова пространства К™". Матрицы неполного ранга образуют в этом пространстве поверхность, что также иллюстрируется рисунком.

3.5. Задачи наименьших квадратов неполного ранга 139 Пример 3.8. Проиллюстрируем описанную выше процедуру с помощью двух матриц неполного ранга, имеющих размер 20 х 10: А1 ранга г, = 5 и Ат ранга гэ — — 7. Пусть А, = У,Е,Ъ',т (1 = 1, 2) — сингулярные разложения этих матриц, где общим размером матриц Ьн Е, и )г, является ранг г, матрицы Аб это тот же способ записи, что и в предложении 3.3. Ненулевые сингулярные числа матрицы А; (или матрицы Е,) указаны крестиками на рис. 3.4 (для А|) и рис.

3.5 (для Ат). Заметим, что А1 имеет пять больших ненулевых сингулярных чисел (все они чуть больше единицы, поэтому соответствующие им крестики сливаются в один, находящийся на правой границе рисунка). В то же время, семь ненулевых сингулярных чисел матрицы Аз на рис. 3.5 распределены на отрезке, левым концом которого является число 1.2. 10 э ьо1. Выберем г;-мерный вектор х',. и положим х; = Ч.х',ч Ь; = А;х, = У,Е;х',; таким образом, х, есть точное нормальное решение задачи ппп,9А;х — Ь;9з.

Рассмотрим теперь последовательность возмущенных задач с матрицами А, + 6А, где возмущение 6А выбирается случайно, но так, чтобы значение 96А(! варьировалось в широких пределах. Полагая ьо! = 10 е, применим к задачам ш1пй(А, + 6А)у; — ЬДз процедуру усеченного ЯЧВ. Сплошные линии на рис. 3.4 и 3.5 указывают значение вычисленного ранга матрицы А, + 6А (т.е. количества найденных сингулярных чисел, б/'ольших, чем ьо1 = 10 э) как функцию от й6Айт (см. верхние графики) и ошибку 9у; — х;5з/9х,)~э (нижние графики). Маь1аЬ-программа, создающая эти графики, находится в НОМЕРАСЕ/Ма$1аЪ/В.ап1с Пебс1епс.ш. Ситуация на рис. 3.4 проще, поэтому мы начнем обсуждение с нее. Матрица А1+6А имеет пять сингулярных чисел вблизи или чуть правее 1 и еще пять сингулярных чисел, не превышающих 96А~)э.

Если Й6А)~а < ьо1, то вычисленный ранг будет для А1 + 6А тем же, что и для Аы т. е, будет равен 5. Ошибка тоже медленно увеличивается от уровня порядка машинного эпсилон (- 10 'е) до приблизительно 10 'о вблизи 96Айз = Фо1. Для ббльших значений )~6А(~т обе величины делают скачок соответственно к 10 и 1. Это согласуется с нашим анализом в предложении З.З, который показывает, что число обусловленности обратно наименьшему ненулевому сингулярному числу, т. е.

наименьшему сингулярному числу, превосходящему ьо1. При 96А9з < ьо1 это наименьшее ненулевое сингулярное число близко к (или даже чуть превышает) 1. Поэтому предложение 3.3 предсказывает ошибку порядка 96А!)з/0(1) = й6А~(м Эта 140 Глава 3. Линейные задачи наименьших квадратов Ранг возмущенной матрицы; исходные сингулярные числа указаны символом х; 1о1=1е-09 ~О О ~О" ОО" ~О" Норма возмущения 1У1.х1 1/1х1 ! ОО ю" ~а" 10-ч Оа" ГО Норма возмущения ОО' Рис.

3.4. Решение задачи шшгн '0(А~ + бА)уг — 6~((О методом усеченного отгР прн го1 = 10 О. Сингулярные числа матрицы Аг указаны крестиками. Значения (фА))О откладываются по горизонтальной оси. На верхнем графике показан ранг матрицы Аг + бА, т. е. количество ее сингулярных чисел, превосходящих 1о1. На нижнем графике представлено отношение суг — хг))т/схг~!О, где хг — решение задачи при бА = О. Ранг возмущенной матрицы; исходные сингулярные числа указаны символом х; 1о1=! е-09 ю О ю" ГО" ю Норма возмущения Ьг-хгйхЫ геч ю" ГО" юч ОО" ю Норма возмущения Рис.