Форсайт Дж., Малькольм М., Моулер К. - Машинные методы математических вычислений (1040536), страница 40
Текст из файла (страница 40)
(Во всяком случае оно предпочтительней быстрого алгоритма, дающего, возможно, неверный результат.) Поскольку ортогональные матрицы сохраняют норму, то ~'г,,'=) и (АРУ' — Ь)',~-.=) т.— (~~. Следовательно, Я'Р сводит общую задачу наименьших квадратов к задаче с диагональной матрицей. Легко видеть, что вектор г, дающий минимум для И, выражается соотношениями зт г =- —, если а.~О. 9, ' l г, произвольно, если ат —— О. Таким образом, Ь из уравнений диагональной формы решаются точно. Остальные уравнения приводят к тому, что вектор-невязка, возможно, отличен от нуля, причем его норма равна ()4' =~ч~г(,', где сумма берется по всем (, для которых а; =О либо (~п.
Обратное преобразование х= —.Рг решает тогда исходную задачу. Если задача имеет недостаточный ранг, то решение, минимизирующее йгй, неединственно. В такой ситуации часто бывает желательно добиться единственности, выбирая решение наименьшей длины. Это достигается, если положить г;=-О, когда О,=-О. В случае полного ранга решение единственно. В НАИМЕНЬШИЕ КВАДРАТЫ Иногда рассматривают модифицированные задачи наименьших квадратов, в которых минимизируется некоторая комбинация ~|г~! и 8х!1 Решения таких задач легко получаются из Я119 посредством методов, описанных в книге Лоусон, Хэнсон (1974). В качестве примера снова рассмотрим задачу с 4т 5 5 5 5 Мы уже видели, что система Ак=Ь несовместна и, следовательно, не может быть решена точно. Диагонализованные уравнения таковы: 35.127г, ж 10.716, 2.465г, ж — 0.872, Ог9 ж — 0.541, 0 ж — 0.193, 0 ж — 0.265. Символ ж означает, что эти уравнения нужно решить в смысле наименьших квадратов.
Первые два можно решить точно, полагая г, = 0.305, г, = — 0.354. Третье уравнение не определяет г,; возьмем г,=О, чтобы получить решение наименьшей длины. Последние три урав- нения показывают, что 8г|Р будет равна 0.54! '+0,1939-80.265'= =ОАОО. Заметим, что это можно вычислить без фактического на- хождения х; однако несколько более точный способ, особенно в случае малых невязок, состоит в вычислении х, а затем ЙЬ вЂ” Ах((. Решение х= 1'г есть х = 0.067 Тогда / 4.406 4.607 4.808 5.009 5.210 93. ПРИЛОЖЕНИЯ ззз что довольно близко к заданному вектору Ь. Опять-таки к вектору х можно прибавлять любые кратные последнего столбца матрицы Р', получая другие решения задачи наименьших квадг ратов, имеющие, однако, ббльшую длину. Псевдообратная матрица Матрица Х с размерами нхьч называется псевдообратной Мура — Пенроуза для матрицы А, если она удовлетворяет следующим четырем условиям: 1.
АХА =А. 2. ХАХ=-Х. 3. АХ симметрична. 4: ХА симметрична. Можно показать, что такая матрица Х всегда существует и единственна. Если А — квадратная иевырожденная матрица, то Х= =А-', очевидно, удовлетворяет этим условиям. Если А — прямоугольная и имеет полный ранг, то Х==(АРА) 'А'. Если А не имеет полного ранга, то такого простого представления нет, Для псевдообратной матрицы часто используется обозначение А'. Чтобы применить 5РВ для вычисления псевдообратной матрицы, начнем со случая 1х!, т. е.
с псевдообращения числа: ! —, если а~0, — о 1 О, если о=О. Заметим, что оФ =-0 либо 1 и что четыре условия Мура — Пенроуза выполняются очевидным образом. Теперь можно проверить, что псевдообратной к диагональной т х и-матрице »т, з нлимвньшив квлдехты 234 будет диагональная ахи-матрица 0 Отметим, что ЕХ~ =-Г„, т. е. это матрица с а единицами на диагонали и нулями в прочих позициях.
Наконец, используя трехчленный аналог известного тождества (АВ) '=-В-'А-', можно получить псевдообратную произвольной матрицы из ее сингулярного разложения по формуле Аг= рХ~цт Легко проверить, что при таком определении четыре уравнения Мура — Пенроуза выполняются. Если матрица А — квадратная и невырожденная, то все ее сингулярные числа отличны от нуля и, следовательно, Х~= =Х-' и А" =А-'. Псевдообратная матрица связана с задачей наименьших квадратов тем обстоятельством, что вектор х наименыпей длины, минимизирующий 11Ах--Ь11, можно выразить формулой к=А" Ь. Если основное соотношение А =- итт расписать поэлементно, то для 1=-1,..., т и 1=1, ..., и справедливо л а,, = У, а„и,Рты Р=! Конечно, суммирование можно ограничить теми значениями й, для которых охай.
Если определение Аз р~~цг также расписать, то для 1=-1, ..., п и 1=-1, ..., т ' хч,ь а;т = ах --~ с Здесь ограничение суммы на ненулевые о„становится необхо-; димым, Нужно подчеркнуть, что элементы псездообратной матрицы не являются непрерывными функциями элементов данной матрицы. В псевдообратной для меняющейся матрицы будут скачки всякий Рз пРилОжения 235 раз, когда матрица меняет ранг. Поэтому в вычислительной практике важно сравнивать малые сингулярные числа с некоторой границей, которая отражает ошибки в исходных данных и вычислениях. Эффективную певвдообратную матрицу можно определить посредством такой границы т следующим образом: ! —, если а.
т, о'=1 в ' 1 0 в противном случае. Это приводит к матрице А, имеющей элементы а >с сс Матрица Х=-Ас удовлетворяет лишь трем последним условиям Мура — Пенроуза, т. е. ХАХ =-. Х; ЛХ симметрична; ХА симметрична. Первое условие заменяется неравенством 11А ХЛ вЂ” А (~(т, где используется определенная выше норма. Если т~О, то матрица Х, удовлетворяющая этим чстырем соотношениям, может быть неединственной; однако Х= — А',с имеет то дополнительное свойство, что она минимизирует 11Х11.
Именно эта способность учитывать неточность входных данных и приближенность вычислений делает характеризацню псевдообратной матрицы посредством сингулярного разложения предпочтительней других подходов. Возвращаясь к нашему примеру, находим — 0.247 — 0.133 — 0.020 0.093 0.207 мй О 067 О 033 0 ООО 0 033 0 067 0.113 0.067 0.020 — 0.027 — 0.073 Заметим, что в этом примере АсА отнюдь не похожа на еди- ничную. Аппроксимирующие матрицы Если А — умеренно большая матрица, с:апример если т и и имеют порядок нескольких сотен, то хранение всех элементов А в машине может потребовать зн чительной доли всей доступной памяти.
Может оказаться желательным аппроксимировать матри- ». НАИМЕНЬШИЕ КВАДРАТЫ 23б цу А более «простой» матрицей, требующей меньше памяти. Разумно понимать под проси»ой матрицей матрицу сравнитечьно невысокого ранга, т. е. имеющую лишь несколько независимых столбцов. Кроме экономии памяти, такие матрицы могут давать важную информацию относительно конкретных задач. Уравнение А =(7Х'Р'г можно переписать в виде Л = о,Е,+о,Е»+... +а, Е„, где Е;=и;от» есть внешнее иро»шведение столбца матрицы (7 и соответствующего столбца матрицы Г. Каждая Е, — это простая матрица ранга 1, для хранения которой требуется лишь ли+и вместо ти ячеек памяти, Кроме того, произведение Е,х можно сформировать для любого вектора х всего лишь за т+и умножений, а не за иип, Чтобы оценить степень зкономии, заметим, что для т=и=300 будет т+п — 600, в то время как гни=90 000.
Сумма квадратов всех элементов матрицы Е; равна 1, так что ни одна из этих матриц не является численно более или менее важной, чем любая другая. Зато некоторые из коэффициентов о; могут быть значительно меньше других, и потому написанный ряд для А можно оборвать после нескольких членов. Сокращенная сумма определяет приближение к А, для которого, однако, хранение и оперирование обходятся дешевле. Разумеется, полезность этой процедуры зависит от конкретной матрицы Л и распределения ее сингулярных чисел. Недавно подобные приближения были с успехом применены в таких различных областях, как цифровая обработка образов и криптография; см. статьи Эндрюс, Паттерсон (!975) и Моулер, Моррисон (1977).
Вычисление определителей и нахождение вырожденных матриц Пусть А(г) — квадратная матрица, зависящая, быть может, сложным образом, от некоторого параметра г. Рассмотрим задачу нахождения того значения или тех значепий г, для которых Л (г) вырождена. Другими словами, найти г, для которых определи- '* тель А(г) равен нулю. Используя символ г(е1 для обозначения определителя, можем ~ написать де1(Л) = бе1((7) де1(Х) бе1(Р'т). Определитель ортогональной матрицы равен ~1, а определитель диагональной матрицы есть просто произведение ее диагональных элементов. Таким образом, г(е1(А)= — -~о, о, ... о„.
9.х ПРИЛ07КЕНИЯ 237 Вычисление определителя может оказаться непростым делом, потому что значения его меняются в очень широком интервапе, что часто приводит к переполнению либо машинному нулю. Для многих задач, формулируемых в терминах определите. лей, не требуется фактическое значение определителя, а лишь некоторое указание на то, когда он равен нулю. Теоретически возможно вычислить ~де1(А)~, беря произведение сингулярных чисел А. Однако в большинстве ситуаций достаточно использовать лишь наименьшее сингулярное число, избегнув тем самым переполнений, машинных нулей и прочих вычислительных затруднений.