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

3.5. Решение задачи ш1птх '0(Аг + бА)уг — ЬО11О методом усеченного БЪг11 при го1 = 10 О. Сингулярные числа матрицы АО указаны крестиками. Значения ~)бА!)О откладываются по горизонтальной осн. На верхнем графике показан ранг матрицы Ат+ бА, т.е. количество ее сингулярных чисел, превосходящих го!. На нижнем графиае представлено отношение 11уз — вг 5/!)з)(О, где зз — решение задачи при бА = О. 141 3.5. Задачи наименьших квадратов неполного ранга хорошо обусловлениал ситуация подтверждается малыми значениями ошибки слева от точки ]]бА][г — — со1 на нижнем графике рис.

3.4. С другой стороны, при ]]бА]]2 > со1 наименьшее ненулевое сингулярное число имеет порядок ОЯбА]]г). Вследствие малости этого числа, ошибка делает скачок к уровню ]]бА]]2/0[]]бА]]2) = 0(1), что видно из части нижнего графика рис. 3.4, находящейся правее точки ]]бА]]2 —— со1.

Ненулевые сингулярные числа матрицы Аг показаны на рис. 3.5 крестиками. Наименьшее из них равно 1.2.10 э, что чуть больше, чем со1. Поэтому при ][бА]]2 < со1 прогноз для ошибки есть ][бА]]2/10 э, что возрастает до 0[1) при ]]бА]]2 = со1. Это подтверждается нижним графиком рис. 3.5. О 3.5.2. Решение задач наименьших квадратов неполного ранга методом 14К-разложения с выбором главного столбца Более экономичной, но подчас менее точной альтернативой методу БЧП является 142-разложение с выбором главного снголбца. Пусть в точной арифметике ранг матрицы А равен т < п, и пусть ее первые т столбцов линейно независимы.

Тогда 14В;разложение А имеет вид ть11 ть12 1 А=ухе=Я ~ 0 0 0 0 где т хе-подматрица Ам невырожденна, а Л12 — подматрица размера т х (п — т). При машинных вычислениях мы могли бы надеяться получить матрицу вида В= 0 Лгг с блоком Нгг очень малой нормы [порядка е]]А]]г). В этом случае можно было бы просто положить Лгг = 0 и минимизировать ]]Ах — Ь]]2 следующим образом: пусть [Я, Я] — квадратная ортогональная матрица, тогда ]]Ах — Ь]]~ ~= -, (Ах — Ь) = ]]Дх — Ятб[~г + ][ЯтЬ]~г 3 ° . а=де,,е,],*=[,*! ° .„, „.я=["" ""].Фу„, ]]Ах — Ь]]2 = ]]В11х1+Лггхг — Я, Ь]]2+ ]]Я2Ь]]2+ ]ф Ь[]2 мииимизируется вектором х = ~ " ' ' 'г 2' при произвольном хг. Х2 Отметим, что выбор хг = 0 необязательно минимизирует ]]х]]2, но это — разумный выбор, особенно если матрица Н11 хорошо обусловлена, а произведение В11 хь12 МЗЛО.

142 Глава 3. Линейные задачи наименьших квадратов К сожалению, этот метод не надежен, так как В может быть близка к матрице неполного ранга, даже если ни один блок Взз в ней не мал. Например, для двухдиагональной п х п-матрицы имеем н м(А) 2 ", но А = Я В с Я = 1 и В = А, и ни один блок Вез в В не мал. Чтобы избежать подобньгх ситуаций, когда не удается распознать неполноту ранга, можно применить ЯВ-разложение с выбором главного столбца. Это означает, что вычисляется разложение АР = ЯВ, где Р— некоторая матрица- перестановка. Идея состоит в том, чтобы на шаге г (где г меняется от 1 до числа п столбцов матрицы А) выбрать столбец с наибольшей нормой в еще не приведенной части А (т.е.

в подматрице, стоящей на пересечении столбцов Ь 1+ 1,...,п и строк з, з + 1,...,т) и переставить его с з'-м столбцом. Далее выполняется обычное отражение, аннулнрующее элементы 1+ 1,..., ш столбца 1. Эта стратегия продиктована желанием получить как можно лучше обусловленный блок В~~ и как можно меньшие элементы в блоке Взз. Пример 3.9. Применим С4В;разложение с выбором главного столбца к матрице предыдущего примера (число .5 на главной диагонали и 1 на первой наддиагонали) при п = 11. Получим Вы и = 4.23 10 4, что является хорошим приближением к т ы(А) = 3.66.10 4. Заметим, что всегда В„„> н м(А).

Действительно, замена элемента В„„на 0 приводит к уменьшению ранга; между тем, о ьз1А) — это ноРма наименьшего возмУщениЯ, способного понизить значение ранга. О В общем случае можно лишь утверждать, что — ~'"'2) ( 2", но обычно Воо является приемлемым приближением к н ы(А). Наихудший случай, однако, столь же плох, как и наибольший рост главного элемента в методе СЕРР.

В самое последнее время предметом интенсивных исследований стали более сложные схемы, так называемые (1В-алгоритмы с выявлением ранга (гапйтеоеайпо ОВ а1догчйтз). Эти алгоритмы распознают значение ранга более надежно и подчас быстрее, чем ЯВ-разложение с выбором главного столбца.

Алгоритмы с выявлением ранга построены в ~28, 30, 48, 50, 109, 126, 128, 150, 196, 236]. Мы вернемся к их обсуждению в следующем разделе. ЯВ-разложение с выбором главного столбца реализовано ЬАРАСК-подпрограммой яйечр1. 1 АРАСК включает в себя еще несколько аналогичных разложений: ЯВ (вйег91), Щ (вйе19ь) и ЯЬ 1ябе91Х).

Будущие версии ЬАРАСК'а будут содержать улучшенные варианты ЯВ-разложения. 143 3.6. Сравнениепроизводительности методов 3.6. Сравнение производительности методов для решения задач наименьших квадратов Какой метод быстрее всего решает плотные задачи наименьших квадратов? Из обсуждения в равд. 3.2 следует, что наиболее быстрым является метод нормальных уравнений, за ним следует метод С4Рь и только потом ЯЧП. Если А— хорошо обусловленная матрица, то метод нормальных уравнений по точности сравним с прочими методами.

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

Мы дадим краткое изложение недавней работы [206], где было проведено сравнение производительности нескольких алгоритмов друг с другом, а также с наиболее быстрым устойчивым алгоритмом для задач полного ранга, а именно С4й-разложением без выбора главного столбца, реализованным на основе отражений, как описано в равд. 3.4.1, и с иерархией памяти, оптимизированной в соответствии с вопросом 3.17. Сравнение проводилось на компьютере 1ВМ ВЯ6000/590 в режиме вычислений с двойной точностью.

Среди сравниваемых методов были Яй-алгоритмы с выявлением ранга„ упоминавшиеся в разд. 3.5.2, а также различные реализации метода ЯЧП (см. разд. 5.4). В тестировании использовались матрицы различных порядков и с разнообразным распределением сингулярных чисел. Мы обсудим результаты, полученные для двух типов матриц: тнп 1: случайные матрицы с элементами, равномерно распределенными на отрезке [ — 1, 1]; тип 2: матрицы, сингулярные числа которых геометрически распределены на отрезке [1,е] [иначе говоря, г-е сингулярное число равно у', где 7" = е).

Матрицы типа 1, в общем случае, хорошо обусловлены, а матрицы типа 2 близки к матрицам неполного ранга. В тестировании участвовали квадратные матрицы как малого, так и большого порядка [соответственно, и = т = 20 и и = т = 1600). Использовались именно квадратные матрицы, потому что для т х п матрицы А с т, значительно ббльшим, чем п, имеется более дешевая альтернатива: выполнить (~В;разложение в качестве подготовительного шага, а затем применить к полученной матрице В алгоритм с выявлением ранга либо метод ЯЧП.

[Так делается в ЕАРАСК'е.) Если т » и, то стоимость начального ЯВ-разложения доминирует над ценой последующих операций с и х п-матрицей В, а потому все указанные алгоритмы имеют приблизительно одну и ту же стоимость. Наиболее быстрым вариантом С4К-разложения с выявлением ранга был метод из [30, 196]. Для матриц типа 1 отношение времени работы этого метода к времени ЯВ-алгоритма без выбора главного столбца изменялось от 3.2 при и = т = 20 до всего лишь 1.1 при и = т = 1600. Для матриц типа 2 это отношение изменялось от 2.3 [и = т = 20) до 1.2 (п = т = 1600).

Подобное же отношение для программы йбе<рг из текущей версии ЬАРАСК'а колебалось менялу 2 и 2.5 для матриц обоих типов. 144 Глава 3. Линейные задачи наименьших квадратов Среди вариантов БЧВ наиболее быстрым был метод из [58], хотя при п = т = 1600 примерно столь же быстро работал алгоритм, основанный на стратегии «разделяй и властвуй» (см. разд. 5.3.3). (Этот алгоритм, кроме того, использовал гораздо меныпую память.) Отношение времени работы БЧВ к времени ЯН-алгоритма без выбора главного столбца изменялось для матриц типа 1 от 7.8 (при п = т = 20) до 3.3 (при п = т = 1600). Для матриц типа 2 это отно|пение изменялось от 3.5 (п = т = 20) до 3.0 (п = т = 1600).

Для алгоритма 66е1ав из текущей версии ЬАРАСК'а подобное же отношение колебалось от 4 (матрицы типа 2 при п = т = 20) до 97 (матрицы типа 1 при п = т = 1600). Это огромное замедление следует, по-видимому, отнести к эффектам иерархической организации памяти. Мы видим, таким образом, что при решении задач неполного ранга имеется баланс между надежностью и скоростью: ЧВ;разложение без выбора главного столбца — это самый быстрый, но и наименее надежный метод, алгоритм БЧ — самый медленный, но зато самый надежный, а методы цН с выявлением ранга находятся посредине между тем и другим. Если т » п, то стоимость всех алгоритмов примерно одинакова. Выбор конкретного метода зависит от сравнительной важности скорости и надежности для пользователя.