ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2007, том 47, < 12, с. 1979-1987
УДК 519.612
О НЕКОТОРЫХ МЕТОДАХ ТИПА СОПРЯЖЕННЫХ НАПРАВЛЕНИЙ ДЛЯ РЕШЕНИЯ ПРЯМОУГОЛЬНЫХ СИСТЕМ
ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ1)
© 2007 г. Л. Ф. Юхно
(125047 Москва, Миусская пл., 4а, ИММ РАН) e-mail: yukhno@imamod.ru Поступила в редакцию 24.05.2007 г.
Исследуются вопросы применения модификаций некоторых известных методов типа сопряженных направлений для решения систем линейных алгебраических уравнений с прямоугольными матрицами. Показывается преимущество модифицированных методов по сравнению с исходными вариантами относительно накопления ошибок округления, особенно при плохой обусловленности матриц. Приводятся примеры эффективного использования этих модификаций для некоторых плохо обусловленных задач достаточно больших размеров. Библ. 3. Табл. 6.
Ключевые слова: система линейных алгебраических уравнений, метод сопряженных направлений, численная устойчивость.
1. ВВЕДЕНИЕ
При решении системы линейных алгебраических уравнений (СЛАУ) методами типа сопряженных направлений важным этапом является построение ортогональной в некоторой выбранной метрике системы векторов, по которой раскладывается искомое решение (см. [1]). Поскольку речь идет о конечномерном пространстве, то количество таких векторов при точных вычислениях конечно, т.е. такие методы являются конечными. При наличии ошибок округления, особенно в плохо обусловленных задачах, процесс построения такой "ортогональной" системы векторов может стать бесконечным, т.е. при этом методы можно рассматривать как чисто итерационные. В связи с этим важным свойством метода является его устойчивость в смысле невозрастания соответствующего функционала ошибки в процессе итераций.
Разработке эффективных модификаций для некоторых методов такого типа посвящены работы [2], [3]. Модифицированные методы при точных вычислениях полностью совпадают с исходными, т.е. также являются конечными, а при неточных - по той же причине становятся итерационными. Однако итерации строятся таким образом, что функционал ошибки на каждом шаге не возрастает, что гарантирует их устойчивость. Там же приводятся результаты расчетов для некоторых задач с квадратными матрицами.
Настоящая работа ставит целью применение некоторых из этих модифицированных методов и решение возникающих при этом вопросов к системам с прямоугольными матрицами. Рассматриваются как переопределенные, так и недоопределенные СЛАУ. Показывается преимущество применения методов непосредственно к исходной прямоугольной системе.
Приводятся результаты расчетов для некоторых модельных задач из теории аппроксимации функций. Возникающие в них СЛАУ характеризуются своей плохой обусловленностью. Расчеты демонстрируют эффективность модифицированных методов и их применимость для решения достаточно, на наш взгляд, сложных задач.
2. ПОСТАНОВКА ЗАДАЧИ
Будем рассматривать СЛАУ вида
Ax = b, (2.1)
где матрица A е С х m и вектор b е С заданы, вектор x е С искомый. Предполагается, если спе-
^ Работа выполнена при финансовой поддержке РФФИ (код проекта 06-01-00749).
1979
циально не оговорено, что А - матрица полного ранга, т.е. гапкА = тт(ш, п). Для пары векторов в Сп и в €ш будем использовать стандартное скалярное произведение
( и, V) = V* и.
При п > ш систему (2.1), возможно несовместную, заменим системой, полученной после первого преобразования Гаусса:
А * Ах = А* Ь. (2.2)
Так как det(A*A) > 0, то эта задача имеет единственное решение, которое будем обозначать через х*. Известно, что этот вектор дает минимум функционала (Ах - Ь, Ах - Ь).
При п < ш система (2.1) имеет ш - п-мерное линейное многообразие решений. Выделим из него решение (обозначим его через х**), ортогональное подпространству решений соответствующей однородной системы, т.е. решение, минимальное по норме. Используем для этого второе преобразование Гаусса (см. [1]). А именно, будем брать векторы, представимые в виде х = А*у. Относительно у из (2.1) получается задача
АА* у = Ь. (2.3)
Уравнение (2.3) имеет, и притом единственное, решение у*, которое и дает искомое решение х** = А*у* задачи (2.1).
Однако можно находить решения х* при п > ш и х** при п < ш непосредственно из системы
(2.1) без ее предварительного преобразования (см. [1]). Для этого достаточно применить к этой системе один из следующих методов типа сопряженных направлений: при п > ш - это метод ор-тогонализации столбцов и метод А*А-минимальных итераций; при п < ш - это метод ортогонали-зации строк, метод АА*-минимальных итераций и метод Крейга. Здесь перечислены те методы, которые рассмотрены в [2], [3] на предмет их оптимизации.
В настоящей работе остановимся на применении метода А*А-минимальных итераций для случая п > ш и метода АА*-минимальных итераций для случая п < ш. Иногда для краткости будем называть их А*А-метод и АА*-метод.
В каждом из этих случаев возможно модифицировать указанные методы так, что они будут являться численно более устойчивыми по отношению к влиянию ошибок округления, чем исходные.
В работе рассматриваются вопросы использования этих модификаций для решения задач с прямоугольными матрицами.
3. МЕТОД А*А-МИНИМАЛЬНЫХ ИТЕРАЦИЙ
Будем применять метод А *А-минимальных итераций к решению переопределенной системы (2.1), т.е. при условии п > ш. Рассмотрим формулы модифицированного варианта метода в этом случае (в [3] они приведены для ш = п).
Для некоторого произвольно взятого вектора 5 е С строится последовательность векторов
g2, ..., где gi е Сш, по следующим формулам:
£1 = 5, v 1= g2 = А *v 1- у 1 gl; (3.1) далее для г = 2, 3, ... вычисляются
g. + 1= А* V г - у^ - 6^_!, (3.2)
где
^ (А * VА *v1) 6 (VV.) _ п .. ..
1г = —(-)—, 6. = (-), V. = А£,, V. е С . (3.3)
(V, (V-1, V-1)
Такой метод, как было указано, рассматривается как чисто итерационный. При этом предполагается, что все векторы g1, g2, ... отличны от нуля.
Параллельно при г = 1, 2, ... вычисляются значения
хо = 0, х. = хг-1 + аgг, х. е С", (3.4)
где
а (b, v,) - (х, A* v,)
a, = -(-)-. (3.5)
( v „ v,)
Формула (3.5) составляет суть модификации (см. [2]).
Как уже было отмечено, при точной реализации метода, если векторы g1, ..., gm линейно независимы, вычислительный процесс заканчивается на m-м шаге, т.е. модифицированный метод, как и исходный, является конечным; в результате берется х* = xm. При наличии ошибок округления, особенно для задач с плохо обусловленными матрицами, метод становится итерационным. В [3] показано, что функционал ошибки, который для этого метода имеет вид
Fi(х) = (A(х _ х*), A(х _ х*)), (3.6)
при использовании модификации в процессе итераций гарантированно не возрастает на каждом шаге. Это было подтверждено приведенными там же расчетами (для примеров с квадратными матрицами).
Однако имеется следующее осложнение при решении задач больших размеров. Если матрица A содержит большие по модулю элементы, то ее норма велика, а тогда векторы g, и v, могут сильно расти по норме при возрастании номера i. Это может привести к выходу из допустимого диапазона чисел для разрядной сетки компьютера. Та же неприятность, как выход за границы диапазона и авост, может возникнуть и когда норма матрицы A достаточно мала.
Для устранения этого недостатка предлагается видоизменить формулы (3.1)-(3.5), применяя нормирование векторов vt на каждом итерационном шаге.
Пусть |z | обозначает норму вектора z, |z | = V(z, z). Приведем соответствующие формулы для модифицированного метода с нормированием:
gi = v j = -Agi, -i = |Ag1, g2 = A *v !_ y jgj; (3.7)
M-i
далее для i = 2, 3, ... вычисляются
g, + i = A* v, _ y,g, _ 8,g,_i, (3.8)
где
v, = M Agi, -, = |Ag, Y, = (A * vM A * ^ ), 8, = (3.9)
Параллельно при , = 1, 2, ... вычисляются
хо = 0, х, = х, _i + a,g,, (3.10)
где
(b, v,) _ (х,_i, A*v,) a, = -. (3.11)
г M,
Приведенные в разд. 5 результаты расчетов показали, что такой модифицированный метод с нормированием позволяет рассматривать матрицы достаточно больших размеров, избегая указанных неприятностей.
Сделаем замечание о выборе начального вектора 5.
При описании метода A*A-минимальныx итераций (см. [1]) предполагалось, что построенная совокупность m векторов s, A*As, (A*A)25, ..., (A*A)m - 15, подлежащая последующей A*A-ортого-нализации, является линейно независимой. Это требование может быть существенно ослаблено специальным выбором начального вектора s. Если взять в качестве вектора s вектор A*b, то достаточно предполагать, что rankA*A = m (т.е. что система (2.2) однозначно разрешима).
Докажем это утверждение. Поскольку A*A-метод эквивалентен A-методу после первого преобразования Гаусса, то для простоты изложения рассмотрим именно A-метод.
Итак, решается система Ax = b с матрицей A е Cm х m, где A = A* > 0. Рассмотрим для матрицы A аннулирующий многочлен наименьшей степени
ф( A ) = Ap + y i Ap + ... + Y pI = 0.
Согласно теореме Гамильтона-Кэли, для этой матрицы существует аннулирующий многочлен ш-й степени, поэтому р ^ ш. Так как А > 0, то detA Ф 0. Следовательно, ур Ф 0. Тогда можно выразить А-1 следующим образом:
А"1 = е11+02 А + ... + ерАр-1.
(Здесь ех = -Ур - х/Ур, е2 = -Ур - 2/Ур, ., ер = -1/Ур.) А отсюда следует
х = А'1 ь = е1 ь + е2Аь +... + е рАр-1ь .
Таким образом, для этого метода решение х лежит в подпространстве, порожденном векторами Ь, АЬ, ..., Ар - 1Ь. Следовательно, для метода А *А-минимальных итераций (см. (2.2)) решение х лежит в подпространстве, порожденном векторами А*Ь, А*А(А*Ь), ..., (А*А)р - 1(А*Ь), которые совпадают с указанной для этого метода совокупностью векторов при 5 = А*Ь. Поэтому если эти векторы линейно зависимы, то процесс вычислений просто досрочно прекратится.
Приведенные рассуждения предполагают точную реализацию метода. В силу его итеративного обобщения в виде предложенной модификации окончательный вывод верен и для приближенной реализации вычислений.
Таким образом, в качестве начального вектора 5 естественно брать вектор А*Ь.
Отметим также, что минимизация функционала (3.6) на каждом шаге при использовании А*А-ме-тода эквивалентна минимизации функционала нев
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.