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

АтАх = АтЬ. Это система из и линейных уравнений от и неизвестных, называемая системой нормальных уравнений. Почему вектор х = (АтА) 'АтЬ дает минимум функции 11Ах — ЬЦ? Это можно показать, замечая, что гессиан АтА положительно определен; следовательно, функция строго выпукла и всякая ее критическая точка есть точка глобального минимума. Или же, полагая х' = х+е, можно упростить функцию к виду суммы квадратов: (Ах' — Ь)г (Ах' — Ь) = (Ае + Ах — Ь) (Ае + Ах — Ь) = (Ае) (Ае) + (Ах — Ь)т(Ах — Ь) + 2(Ае)т(Ах — Ь) = 11АЦ, '+ 11Ах — Ь!), '+ 2.т(АтАх — АтЬ) = )(Ае)(г + ()Ах — Ь|)~~. Ясно, что сумма минимизируется выбором е = О.

Это не что иное, как теорема Пифагора, поскольку невязка т = Ах — Ь ортогональна к надпространству, натянутому на столбцы матрицы А, т. е, О = А т = А Ах — АтЬ. Это иллюстрируется рисунком; показанная на нем плоскость есть линейная оболочка столбцов матрицы А, поэтому векторы Ах, Ае и Ах' = А(х + е) принадлежат ей. Поскольку матрица Аг А — симметричная и положительно определенная, для решения нормальных уравнений можно применить разложение Холесского. Общая стоимость вычисления АтА и АтЬ и последующего разложения мат- 117 3.2. Л4атричные разложения рицы А А составляет и т + -'пз + 0[и~) флопов.

Так как т > п, то в общей стоимости преобладает цена пгт формирования матрицы АтА. 3.2.2. ЯВ.-разложение Теорема 3.1. [ЯН-разложение). Пусть А — матрица размера т х и, причем т > и. Предполооюим, что А имеет полный столбцевой ранг. Тогда существуют и едипственньг т хи-матрица Я с ортонормированными столбцами (т. е. Ц Я = Х„) и верхнетреугольная и х п-матрица Я с полоокительными диагональными элементами тп, гаакие, что А = ЯП. Доказательство . Будут даны два доказательства этой теоремы. Прежде всего, теорема представляет собой переформулировку процесса ортогонализации Грама — Шмидта [139]. Если этот процесс применить к столбцам а; матрицы А = [аз,аг,...,а„] в порядке возрастания их номеров, то будет получена система ортонормированных векторов оы...,9„, дающих другой базис того же надпространства.

Эти ортонормированные векторы являются столбцами матрицы с/. В процессе Грама — Шмидта вычисляются также коэффициенты г, = 9 а, в выражении каждого столбца как линейной комбинации векторов т ды..., йг: а; = 2 1 гзч с7 . Эти числа г;, суть элементы матрицы Л. Алгоритм 3.1. Классический алгоритм Грама — Шмидта (ССЯ) и модифицированый алгоритм Грома — Шмидта (МСЯ) для вычисления разложения А=ЯЛг /от г' = 1 ~о и /* вычислить г'-е столбцы матриц Я и П */ ог = аг /ог у' = 1 Со г — 1 /* вычесть из а; компоненту в направлении 41 */ г,; = д~аг СОЯ дг = о, — гйог епд /ог ги — [[чг [[2 1/ га = 0 /' а, линеупо зависит от аы..., а, 1 */ выход г~д 1/ д; = дг/га Предоставляем читателю в качестве упражнения проверку, что две формулы для г, в этом алгоритме математически эквивалентны [см. вопрос 3.1). Если А имеет полный столбцевой ранг, то ни один из элементов ги не будет нулем.

На рисунке процесс Грама — Шмидта проиллюстрирован для 2 х 2-матрицы А. Второе доказательство теоремы использует алгоритм 3.2, который будет представлен в рэзд. 3.4.1. П К сожалению, в арифметике с плавающей точкой алгоритм СОЯ численно неустойчив, если столбцы матрицы А почти линейно зависимы. Алгоритм МСЯ более устойчив и используется в алгоритмах, приводимых в этой книге далее. Однако, если А плохо обусловлена, то матрица ~ может сильно отличаться от ортогональной, т.

е. величина [[ЯтЯ вЂ” 1[[ значительно больше числа в [31, 32, 118 Глава 3. Линейные задачи наименьших квадратов ЗЗ, 149]. Алгоритм 3.2 в разд. 3.4.1 является устойчивой альтернативой при вычислении разложения А = ДВ (см. вопрос 3.2). С помощью разложения А = ЯВ мы тремя несколько различающимися способами выведем формулу для вектора х, минимизирующего функцию ]]Ах — Ь]]2. Во-первых, всегда можно найти еще т — п ортонормированных векторов, составляющих матрицу ь„так, чтобы матрица ф, ф] была квадратной и ортогональной (например, можно взять дополнительные т — и линейно независимых векторов, составляющих матрицу Х, и применять алгоритм 3.1 к невырожденной п х п-матрице Я, Х]).

Тогда ]]Ах — Ь]]2 = ]](Я, 1е] (Ах — Ь)]Я согласно утверждению 4 леммы 1.7 -т %11х — Ь) О~---)ка ~х О Ь Д. ОтЬ ]]Вх — Я~Ь]Я + Я~Ь]Я ]ЩТЬц2 Систему Вх — Я~Ь = 0 можно решить относительно вектора х, потому что А и В имеют один и тот же ранг и; следовательно, матрица В невырожденна. Поэтому х = Л гЯтЬ, а минимальное значение функции ]]Ах — Ь]]2 равно ]]Оть]]2. Приведем чуть иное рассуждение, не использующее матрицу ф Запишем вектор Ах — Ь в виде .4х — Ь = Яйх — Ь = ~йх — ЯД +1 — ЯЯ )Ь = ~(Ях — ЯтЬ) — (1 — ~~Д~т) Ь. 119 3.2. Матричные разложения Заметим, что векторы Я(Вх — АтЬ) и (1 — ЯЯт)6 ортогональны, так как Я(Вх— Ю~Ь))~((1 — ЮЮ~)6) = ( — Я~Ь)~Я"(1 — ~Ы~))6 = ( — О~Ь)~ЙЬ = О Поэтому по теореме Пифагора, ()Ах — ЬДг = )(Я(Вх — Я ЬЯ~ + )((1 — ЯЯ )6~~э = ~! Вх — ать!!2 +! Н1 — Ют)6!!2, где использовано равенство )Яу(Я = '9уЯ (см. утверждение 4 леммы 1.7). Эта сумма квадратов будет минимальна, если ее первое слагаемое равно нулю, т.

е. если х = В-'Ятб Наконец, есть и третий вывод, который исходит из формулы для решения системы нормальных уравнений: х = (АтА)-'Агб (ВтцтоВ =В В- В"О. 6=В О,'Ь. Позднее мы покажем, что стоимость данного разложения и последующего вычисления решения х составляет 2пгт — г пз флопов, что примерно вдвое превышает стоимость решения нормальных уравнений, если т » и, и примерно равно этой стоимости в случае т = и. 3.2.3. Сингулярное разложение оЧР— это очень важное разложение, используемое, помимо решения задач наименьших квадратов, для многих других целей.

Теорема 3.2. (ЯЧР). Пусть А — произвольная тхп-матприца, причем т > и. Тогда справедливо представление А = ПХЪ'т, где матрица У имеетп размер т х и и удовлетворяет соотношению УтУ = 1, матрица 1т — квадратная порядка и и удовлетворяет соотношению 1тт1' = 1, а Е = бйай(аы...,ап), где аг » а„> О, Столбцы и„...,и„матарицы У назьгваются левыми сингулярными векторами (матрицы А). Столбцы иы..., ев матприцы У называются правыми сингулярными векторами.

Величины аг называются сингулярными числами. (При т ( и нужно рассматаривать БУР матрицы Ат.) Эта теорема допускает следующую геометрическую интерпретацию. Будем рассматривать т х и-матрицу А как линейный оператор, отображающий вектор х Е К" в вектор у = Ах Е К . Тогда можно выбрать ортогональную координатную систему для К" (где ортами осей являются столбцы матрицы И) и другую ортогональную координатную систему для К™ (со столбцами матрицы У в качестве ортов осей) таким образом, чтобы А приобрела диагональный вид (Е), т.е. вектор х = ~,", Ди, отображался бы в у = Ах = ~С,"., аедгиь Иначе говоря, всякая матрица станет диагональной, если выбрать подходящие ортогональные координатные системы в ее области определения и области значений.

Доказательство теоремы 3.2. Применим индукцию по тп и и: предположим, что ЯЧР существует для матриц размера (т — и) х (и — 1), и докажем его существование для т х и-матриц. Будем считать, что А ~ О; в противном случае, 120 Глава 3. Линейные задачи наимевьтвх квадратов -т .А [и Ч= -т -т Тогда т (Ав) (Аи) ]]Аз[Я ]]Ае]]г ]]Ав]]г и Пте = Ути и ]]Ав[]г = О. Мы утверждаем, что и итАг' =- О. Иначе мы имели бы о = ]]А]]г = []ПтАИ[]г > [][1,0,...,0]У~Ах']]г —— ]][о[и~А1т][[г > о, что невозможно.

(Здесь было использовано утверждение 7 леммы 1.7.) (с О ] [т О Итак, У АИ = ~ г-,тАЧ. ~ = ~ А . Применяя индуктивное предположение к матрице А, имеем А = П1Е1)т,, где Вы Е1 и Ъ~ — соответственно т матрицы размера (т — 1) х (и — 1), (и — 1) х (и — 1) и (и — 1) х (и — 1).

Отсюда 0 УгЕз1'~т 0 П~ 0 Е1 0 Ъ'~ О П, О Е, ' О или Это и есть требуемое разложение. П Б'х'Р имеет множество важных алгебраических и геометрических свойств. Наиболее важные из них мы сейчас сформулируем. Теорема 3.3. Пусть А = ПЕЪ'т есть Яхт хи-матприцы А, причем та > и. (Аналогичные факты справедливы для т ( и.) 1.

Пусть А — симметричная матрица с собственными значениями Л, и ортонормированными собственными векторами и;. Другими словами, А = ПЛПт есть спектральное разложение матрицы А; здесь Л сйай(Лы...,Л„), У = [иы...,и„] и ПБт = 1. Тогда $УР матрицы А имеет вид А = УЕ'х'т, где а, = [Л,] и и, = зйп(Лг)ио причем вйп(0) = 1. 2. Собственными значениями симметричной матрицы АтА являются числа ог. Правые сингулярные векторы иг суть соответствующие ортонормированные собстветсые векторы. можно положить Е = 0 и взять в качестве У и Ъ' произвольные ортогональные матрицы.

В качестве базиса индукции рассмотрим случай и = 1 (поскольку т > и). В этом случае положим А = УЕУт, где У = А/]]А]]г, Е = ][А]]г и 1т = 1. Для индуктивного перехода найдем вектор е, такой, что ]]и]]г = 1 и []А[]г — — ][Ае]]г > О. Такой вектор существует в силу определения ]]А]]г = тахг„г, 1][Аи]]г. Положим и =,,л„"„; этот вектор имеет единичную длину. Выберем У и Ъ' таким образом, чтобы У = [и, У] была ортогональной т х тматрицей, а Ъ' = [и, Ъ'] — ортогональной и х и-матрицей.

Теперь имеем 3.2. Матричные разложения 121 3. Собственными значениями симметричной матприцы ААт являются числа о~ и т — п нулей. Левые сингулярные вектортя и, суть ортонормированные собственнтяе векторъц отвечающие собственным значениям ог. Дополнительные тп — п ортогональных векторов можно взять в качестве собственнътх векторов для собственного значения О. О А 1 4. Пусть Н = 1 ~, где А — квадратная матрица, имеющая ЯЪтП А = УЕЪтт.

Положим Е = дйа8(оы..., оп), У = [иы ..,, и„] и Ъ' = [щ,..., и„]. Тогда 2п собственных значений матрицы Н вЂ” это числа хои а соотпветаствующие нормированные собственные векторы имеют вид —, ч'2 ~ з:ън 5. Если А — матприца полного ранга, то решением задачи пппДАх — Ь[]г является вектор х = Ъ'Е ~У~Ь.

6. [[А]]г = оы Если А — квадратпная и невырожденная матрица, то []А '[[г ' = а„и [)А[]г [[А '[)г = Яь. 7. Предположим, что ат » . от > от+т — — .- — — <т„= О. Тогда ранг мат рицтя А равен т. Ее ядро, т. е, надпространство векторов щ таких, что Аи = О, натянуто на столбцтл матприцы Ъ' с номерами от т + 1 до п, т. е.

является подпространстпвом арап(и„еы..., св). Образ матрицы А, т. е. надпространство векторов вида Аю для всех и, натянут на столбцы матрицы У с номерами 1,...,т, т. е. является подпространством арап(им..., ит) . 8. Пустпь Ев ~ — единичная сфера пространства К": Еи ~ = (х е К": [~х]~г = Ц. Пусть А Ео т есть образ сферы Я" т под действием матрицы А: А 5" т = (Ах: х б К" и [~х[[г — — 1). Тогда А ° Ео т является эллипсоидом с центром в нуле пространства К™ и с главными осями а;и,. 9. Положим 1' = [им из,..., и„] и У = [иыиг,...