Содержание к диссертации
Введение
Глава 1. Предварительные сведения 22
1.1. Обобщенный метод минимальных невязок 22
1.2. Обобщенный процесс Ланцоша 26
1.3. Метод MINRES-N 30
Глава 2. Метод MINRES-N2 32
2.1. Основной вариант 32
2.1.1. Детали реализации метода 32
2.1.2. Численные эксперименты 37
2.2. Линейные многочлены от унитарных матриц 41
2.2.1. Детали реализации метода 41
2.2.2. Численные эксперименты 45
Глава 3. Малоранговые возмущения эрмитовых систем 50
3.1. Нормальные возмущения 51
3.2. Анормальные возмущения ранга 1 59
3.3. Анормальные возмущения большего ранга 62
Глава 4. Восстановление матриц по их одноранговым возмущениям 70
4.1. Вспомогательные результаты 71
4.1.1. Решение задачи 4.3 71
4.1.2. Решение задачи 4.4 73
4.2. Восстановление эрмитовых матриц 77
4.2.1. Эрмитова матрица А 78
4.2.2. rank(A- А*) = 1 79
4.2.3. гапк(Л - А*) = 2 79
4.2.4. Нильпотентные матрицы R 80
4.3. Восстановление унитарных матриц 81
4.3.1. Унитарная матрица А 82
4.3.2. гапк(Л*Л-/п) = 1 85
4.3.3. тапк(А*А - /„) = 2 89
4.4. Восстановление комплексных симметричных матриц 94
4.4.1. Симметричная матрица А 95
4.4.2. гапк(Л- Ат) = 1 95
4.4.3. гапЦА- Ат) = 2 95
Литература 100
Приложение 1. Процедура MINRES-N2 (основной вариант) 103
- Обобщенный метод минимальных невязок
- Детали реализации метода
- Анормальные возмущения большего ранга
- Восстановление комплексных симметричных матриц
Введение к работе
Пусть А — заданная п х n-матрица. Фиксируем вектор х и рассмотрим степенную последовательность
х,Ах,А2х,...,Ат 1х . (0.1)
Линейная оболочка векторов (0.1) называется m-м крыловским подпространством, порожденным матрицей А и вектором ж, и обозначается К,т(А,х). Если An х известны из контекста, используется более простой символ Кт.
С ростом т размерность подпространства /Ст растет, пока, в типичном случае, не достигнет числа d(A) — степени минимального многочлена матрицы А. Если в разложении вектора х по корневому базису А отсутствуют какие-либо компоненты, то максимальная размерность крыловского подпространства может оказаться меньше, чем d(A).
Пусть нужно решить линейную систему
Ах = Ь (0.2)
с невырожденной матрицей А Є МП(С). Будем, для простоты, считать, что к системе (0.2) не применяется предобусловливание (или же что (0.2) есть система, полученная предобусловливанием некоторой первоначальной системы). Итерационный метод, в котором т-е приближение хт (т — 1,2,...) выбирается в подпространстве 1Ст(А, Ь), называется методом кры-ловских подпространств.
Наиболее известным среди этой группы методов является метод сопряженных градиентов. Он предполагает, что матрица А — симметричная (эрмитова) и положительно определенная, и описывается следующим образом:
Алгоритм 0.1. Алгоритм сопряженных градиентов: т = 0; х0 = 0; г0 = 6; р\ = Ь;
repeat
m = m+ 1 z = A • pm
Vm = K-l m-l)/&)
xm = xm—і -г vmpm
Гт = I m—I VmZ Pm+1 =Гт+ flm+lPm
пока число rm2 не станет достаточно малым
Здесь хт есть т-е приближенное решение, гт — о — Ахт — его невязка ирт- текущее направление спуска. Вектор хт обладает тем свойством, что вектор ошибки
е = хТ — х (0.3)
{хт — точное решение системы (0.2)), рассматриваемый как функция от х Є /Cm, имеет минимальную норму
\\е\\А = (Ае,еу/2 (0.4)
при х = хт. Норма (0.4) имеет смысл как раз потому, что А положительно определена.
Векторы хт, гт и вспомогательные векторы рт вычисляются посредством коротких (двучленных) рекурсий, что является одним из наиболее привлекательных свойств метода сопряженных градиентов. Это свойство удается сохранить для систем (0.2) с симметричными (эрмитовыми) незна-коопределенными матрицами А, изменяя способ выбора хт в подпространстве /Ст. Вектор хт можно искать, пользуясь условием Галеркина, т.е. так, чтобы соответствующая невязка гт была ортогональна JCm. Это приводит к методу SYMMLQ [1]. Другая возможность — найти хт, для которого евклидова длина невязки гтІ2 минимальна в /Ст. ЭТОТ принцип минимальных невязок определяет метод MINRES [1].
Предположим теперь, что матрица системы (0.2) несимметричная (неэрмитова) и по-прежнему невырожденна. Можно ли и для этого случая построить метод типа сопряженных градиентов, управляемый короткими рекурсиями? Этот вопрос, актуальный в начале 1980-х годов (см. [2]), был решен в СССР В. В. Воеводиным и Е. Е. Тыртышниковым [3, 4] и — независимо и несколько позже, — американскими математиками Фабером и Мантойфелем [5]. Нам удобно описать полученный этими авторами результат, пользуясь терминологией статьи [5].
Возвращаясь к методу сопряженных градиентов, заметим, что векторы спуска Рт (т = 1,2,...) ортогональны относительно билинейной формы задаваемой матрицей А (круглыми скобками обозначено стандартное скалярное произведение в Сп). Если в качестве основной принята обычная евклидова метрика в С", то вместо ортогональности в смысле (0.5) говорят об Л-сопряженности (или просто сопряженности) векторов рт, откуда и происходит название метода.
И в общем, несимметричном случае построение экономичного метода типа сопряженных градиентов эквивалентно указанию короткой рекурсии вида
т pm+l = Арт + J2 VimPi , (0.6)
i=m—s
связывающей векторы спуска рт (т = 1,2,...). Эти векторы должны давать базисы крыловских подпространств /Ст, ортогональные относительно той или иной формы
рв(х,у) = (Вх,у), В = В 0 . (0.7)
Ортогональность векторов рт обеспечивает минимальность соответствующих векторов ошибки ет = х\ — хт в Б-норме етв. Например, выбор В = А А отвечает принципу минимальных невязок. Однако, для простоты, мы в дальнейшем полагаем В = 1п.
Будем говорить, что матрица А принадлежит классу CG(s), если для любой системы (0.2) ортогональные базисы крыловских подпространств К.т могут быть построены посредством рекурсии (0.6). Результат Воеводина-Тыртышникова-Фабера-Мантойфеля состоит в описании класса CG(s).
Как известно, матрица А Є Мп(С) является нормальной тогда и только тогда, когда сопряженная матрица А представима многочленом от А:
А = р{А) . (0.8)
Степень многочлена р можно выбрать так, чтобы она не превосходила п—1. Как правило, минимальная степень р и равна п — 1, хотя для некоторых типов нормальных матриц р может быть многочленом малой степени. Так, p(z) = z, если А — эрмитова матрица.
Назовем А s-пормальной матрицей, если в соотношении (0.8) р есть многочлен степени s.
Теорема 0.1 [5]. Матрица А Є Мп(С), где п s + 2, тогда и только тогда принадлежит классу CG(s), когда А является s-нормальной матрицей либо степень d(A) ее минимального многочлена не превосходит s + 2.
К этому результату нужно добавить следующее описание s-нормальных матриц:
Теорема 0.2 [5]. Пусть А есть s-нормальная матрица. Тогда:
1) если s 1, то А имеет не более s2 различных собственных значений;
2) если s — 1, то А — матрица вида
А = аН + (31п, (0.9)
где Н — эрмитова матрица, ааи/J- комплексные числа.
Смысл теорем 0.1 и 0.2 в том, что методы типа сопряженных градиентов, управляемые рекурсией вида (0.6), возможны лишь для нормальных матриц с небольшим числом различных собственных значений либо для матриц, мало отличающихся от эрмитовых (см. (0.9)). Этот, по-существу отрицательный результат способствовал огромной популярности обобщенного метода минимальных невязок GMRES, предложенного для несамосопряженных систем Саадом и Шульцем [6].
В основе GMRES лежит метод Арнольди унитарного приведения матрицы к форме Хессенберга. Векторы Арнольди, конструируемые в этом методе, составляют ортонормированные базисы крыловских подпространств Кт(А,Ь).
GMRES указывает удобный способ реализации принципа минимальных невязок по отношению к этим подпространствам на основе информации, поставляемой процессом Арнольди. Поскольку матрица, выстраиваемая в этом процессе, хессенбергова, а не ленточная, длина рекурсии, управляющей векторами спуска (т.е. векторами Арнольди), не фиксирована, а линейно растет с ростом индекса т. Растет и количество арифметической работы, выполняемой на шаге процесса. Эти два обстоятельства составляют слабые стороны GMRES в сравнении с методами сопряженных градиентов.
Более подробный анализ GMRES будет дан в разделе 1.1. Этот анализ позволит нам упростить описание оригинального метода MINRES-N, развиваемого в данной работе.
Задача о построении экономичных методов типа сопряженных градиентов приобрела новый импульс в начале 1990-х годов с публикацией результатов У. Грэгга (см. [7]). Грэгг показал, что для унитарной матрицы U ортонормированные базисы крыловских подпространств можно построить с помощью двух коротких рекурсий:
Алгоритм 0.2
Vm+lPm+l = Upm +
Ifm+lPm Pm+1 = &m+lPm + Im+lPm+l ,
Jm+l = -(Upm,pm), (0.10)
Wl = (l-7m+l2)1/2 Основываясь на этих формулах, Ягеле и Райхель предложили эффективный алгоритм минимальных невязок для систем (0.2) с матрицами вида
A = aU + 0I, (0.11)
где U — унитарная матрица (см. [8]).
Исключая из формул (0.10) векторы рт, получим соотношение
1 ТТ Ут+1( ттт
Pm+1 = Upm t/Pm-1
0"m+l 7m&m+l
+ (7т±1 г + Ъ+Щрп . (0.12)
lm&m+l &m+i
Для матрицы вида (0.11) аналогичное соотношение выглядит так:
Pm+1 = Apm —Apm-i
aam+i Ot-lmOm+l
. (Im+l m P 7m+l7m\ . 7m+lO"m _ /n 1 o\
+ ( + )Pm H aPm-1 • (0.13)
lm&m+l &Om+l ( m+l Im m+lP
Эти соотношения схожи с (0.6) соответственно при S = 0 и s = 1. В то же время унитарная матрица U, все собственные значения которой различны, не может быть s-нормалыюй ни для какого малого значения s.
Приведенные факты все же не противоречат теореме 1. По сравнению с (0.6), формулы (0.12) и (0.13) содержат дополнительный вектор (Upm-.i или Apm-i), вычисленный на предыдущем шаге. Это приводит к расширенной постановке вопроса о существовании экономичных методов типа сопряженных градиентов: для каких матриц А ортогональные базисы крыловских подпространств можно строить посредством рекурсий вида
т т
Рт+1 = Е PimAPi - ] VimPi , (0.14)
i=m—t i=m—s
где sat — некоторые малые числа?
В статье [9] даны достаточные условия для положительного ответа на этот вопрос. Назовем нормальную матрицу А (, к)-нормальной, если найдутся многочлены pi{z) и qk(z) соответственно степени I и к, такие, что
A qk(A) = Pl(A) . (0.15)
Если, в частности, к = 0 и qo есть ненулевая константа, мы получаем -нормальные матрицы А. Выше было уже отмечено, что унитарные матрицы U и линейные многочлены (0.11) от них, вообще говоря, не являются -нормальными матрицами для малых L В то же время, равенство
U U = / (0.16)
означает, что всякая унитарная матрица U есть (ОД)-нормальная матрица. Подстановка в (0.16) соотношения (0.11), где коэффициент а предполагается ненулевым, дает
А {А - PI) = РА + {И2 - /52)/ . (0.17)
Таким образом, линейный многочлен от унитарной матрицы является (1,1)-нормальной матрицей.
Выделение класса (, /г)-нормальньіх матриц объясняется тем обстоятельством, что для матрицы А этого класса построение ортогональных векторов спуска с помощью рекурсии вида (0.14) возможно почти для любой правой части b (см. (0.2)). При этом можно считать, что
(s,t) = (,k) (0.18)
(см. [9]).
Следующая теорема из [10] дает описание класса (, &)-нормальных матриц:
Теорема 0.3. Пусть А есть (, /г)-нормальная матрица, причем и к — наименьшие числа, для которых справедливо соотношение (0.15). Тогда:
1) если к + 1, то d(A) 2\
2) если = к + 1, то d(A) 2 или =1,к = 0;
3) если к и свободный член / многочлена qk(z) отличен от нуля, то d(A) к2 + 1;
4) если к - 1 и / = 0, то d{A) к2;
5) если — к - 1 и / = 0, то d(A) к2 или = 0, к = 1;
6) если = к, то с?(Л) А;2 + 1 или = А; = 1.
Поскольку ( , А;)-нормальные матрицы представляют для нас интерес лишь при малых и к, то смысл теоремы 0.3 заключается в следующем: если исключить матрицы с небольшим числом различных собственных значений, то возможны лишь случаи (, к) = (1,0), (, к) = (0,1) и (, к) = (1,1). Первый из них соответствует линейным многочленам от эрмитовых матриц, а второй - унитарным матрицам. Можно показать, что (1,1)-нормальные матрицы, имеющие более двух различных собственных значений, — это линейные многочлены от унитарных матриц.
Итак, по сравнению с 5-нормальными матрицами, расширение класса (, &)-нормальных матриц произошло лишь за счет включения матриц вида (0.11). Оказывается, однако, что рекурсия типа (0.14) возможна для еще более широкого класса так называемых обобщенных (, к)-нормальных матриц. Будем говорить, что А Є МП(С) — обобщенная (, &)-ііормальная матрица, если найдутся многочлены pi{z) и qk{z) соответственно степени
и к, такие, что
гапк(Л %(Л) - ре(А)) = Х п . (0.19)
В частности, значение % = 0 соответствует (, /г)-нормальным матрицам.
Помимо (, &)-нормальных матриц, новый класс содержит их малоранговые возмущения. Это вытекает из следующей теоремы [9]:
Теорема 0.4. Пусть А — матрица вида
Л = Лі + Л2, (0.20)
где А\ есть (, /г)-нормальная матрица и гапкЛг = г. Тогда А является обобщенной (, &)-нормалыюй матрицей, для которой
Х (к + +1)г. (0.21)
Если, в частности, А\ — эрмитова или унитарная матрица, то % 2г. Для матрицы Лі, являющейся линейным многочленом от унитарной матрицы, х Зг.
В [9] показано, что для обобщенной (, /г)-нормальной матрицы Л построение ортогональных векторов спуска с помощью рекурсии вида (0.14) возможно почти для любой правой части 6 в системе (0.2). При этом параметры s и t формулы (0.14) должны удовлетворять неравенствам s- t-k x (0.22)
Мы дали почти исчерпывающее описание положения дел в области конструирования экономичных алгоритмов типа сопряженных градиентов. Для полноты картины следует упомянуть еще о двух приемах, используемых при решении несамосопряженных линейных систем:
1. Симметризация системы, т.е. переход от системы (0.2) к самосопряженной системе
А Ах = A b . (0.23)
К системе (0.23) можно применить стандартный метод сопряженных градиентов. Явное формирование матрицы А А не обязательно, поскольку произведения вида А Ар можно вычислять путем последовательного умножения вектора р на А и А .
2. Использование, наряду с подпространствами JCm(A,b), крыловских подпространств, порожденных сопряженной матрицей А , и построение пар биортогональных базисов в этих двух последовательностях подпространств. На этих принципах построены такие алгоритмы, как метод би-сопряженных градиентов или метод квазиминимальных невязок.
Недостатком первого приема является то обстоятельство, что число обусловленности матрицы А А есть квадрат числа обусловленности матрицы А. Поэтому переход к системе (0.23) не рекомендуется при решении плохо обусловленных систем. Недостатки второго приема связаны с часто наблюдаемой плохой обусловленностью биортогональных базисов и даже иногда их отсутствием (явление, называемое обрывом процесса ). Общий недостаток состоит в необходимости использования матрично-векторных произведений вида A q наряду с произведениями вида Ар. (Справедливости ради, следует отметить, что существуют варианты этих методов, такие, как алгоритмы CGS или Bi-CGStab, не требующие привлечения матрицы А ).
Результаты настоящей диссертации можно рассматривать как вклад в развитие теории экономичных итерационных методов для решения несамосопряженных систем линейных уравнений. Оригинальность используемого нами подхода состоит в том, что вместо стандартных крыловских подпространств принцип минимальных невязок применяется к обобщенным кры-ловским подпространствам. Эти последние были введены в другом контексте в статье X. Д. Икрамова и Л. Эльзнера ([11]; см. также [12]). Если рассмотреть обобщенную степенную последовательность
х, Ах, А х, А2х, А Ах, АА х, А 2х, А3х,... , (0.24)
то линейные оболочки ее конечны отрезков и называются обобщенными подпространствами Крылова, порожденными матрицей А и вектором х.
Мотивом для введения обобщенных крыловских подпространств в [11] было желание найти компактную форму, по возможности близкую к ленточной, для произвольной нормальной матрицы. Как известно, всякую эрмитову матрицу можно привести к трехдиагональному виду посредством конечного процесса, основанного на унитарных подобиях. Для такого приведения могут использоваться различные методы; к технике, используемой в настоящей работе, ближе всего находится метод Ланцоша, тоже относящийся к числу методов (стандартных) крыловских подпространств.
Если матрица А неэрмитова, то она может быть приведена к форме Хес-сенберга посредством метода Арнольди, представляющего собой обобщение процесса Ланцоша. (Этот метод уже упоминался выше в связи с алгоритмом GMRES.) Метод Арнольди не способен извлечь какие-либо выгоды из информации о нормальности матрицы А и дает для нее все ту же заполненную хессенбергову матрицу. Чтобы приблизить компактную форму к ленточной матрице, т.е. но возможности симметризовать ее, нужно, чтобы процесс приведения использовал не только Л, но и сопряженную матрицу А . Это рассуждение и приводит к обобщенным крыловским подпространствам. Более подробное их обсуждение дано в разделе 1.2.
В [11] показано, что компактная форма Я, в которую обобщенный процесс Ланцоша преобразует произвольную нормальную п х n-матрицу А, представляет собой ленточную матрицу с переменной шириной ленты (такие матрицы иногда называют профильными). Матрицу Н можно рассматривать и как блочно-трехдиагональную, причем порядок ее диагонального блока растет с ростом его индексов. Даже в наихудшем случае (т.е. тогда, когда рост порядка происходит скорее всего) общее число ненулевых элементов в Н не превосходит 3\/2п3 2. Эту величину нужно сопоставить с « у ненулевыми элементами хессенберговой матрицы.
Еще более важен следующий факт, также обнаруженный в [И]: пусть нормальная матрица А удовлетворяет соотношению
f(A,A ) = 0, (0.25)
где fix,у) — многочлен степени к. Тогда ширина ленты в матрице Н стабилизируется, начиная с некоторого индекса (зависящего от к). Таким образом, в этом случае компактная форма Н превращается в полноценную ленточную матрицу (или блочно-трехдиагональную матрицу, порядки диагональных блоков которой, начиная с некоторого момента, не возрастают).
Метод MINRES-N, составляющий основное содержание настоящей диссертации, предназначен для класса нормальных матриц, охватываемых (при различных к) соотношением (0.25). Этот класс включает в себя, в частности, эрмитовы матрицы (поскольку Н = Н есть уравнение (0.25) с к = 1), линейные многочлены от них, унитарные матрицы U и матрицы вида A = aU + (ЗІ (см. (0.16) и (0.17)), т.е. ио-существу, весь класс 5-нормальных матриц. При этом он гораздо шире этого последнего класса.
MINRES-N использует ленточность компактной формы Н для построения ортогональных базисов обобщенных крыловских подпространств посредством рекурсий фиксированной длины. К этим подпространствам применяется принцип минимальных невязок, причем условия для его применения гораздо более комфортные, чем в GMRES. Помимо уже упомянутой короткой рекурсии для генерирования векторов спуска (что является наиболее важным обстоятельством), задача наименьших квадратов, решаемая на каждом шаге, имеет ленточную, а не хессенбергову матрицу.
Вторая глава диссертации посвящена более подробному обсуждению метода MINRES-N2, который представляет собой специализацию MINRES-N для случая, когда в (0.25) к = 2. Иными словами, MINRES-N2 ориентирован на нормальные матрицы, удовлетворяющие уравнению вида
сцЛ2 + 2с12АА + с22А 2 + 2с10А + 2с20А + ст1п = 0 , (0.26)
где хотя бы один из коэффициентов квадратичной части сц, С\2 и с22 не равен нулю. При этом случаи
cn2 + c22V0 (0.27)
и
сп = С22 = 0 (0.28)
существенно различаются своей реализацией. Заметим, что первый из них включает в себя нормальные матрицы, спектры которых расположены на эллипсе, гиперболе, параболе или паре прямых, а второй — нормальные матрицы со спектрами, сосредоточенными на окружностях (т.е. линейные многочлены от унитарных матриц). Случай (0.27) разбирается в разделе 2.1, а случай (0.28) — в разделе 2.2.
Отметим, что изложение в этой главе основано на наших работах [13, 14].
В третьей главе диссертации показано, что метод MINRES-N, разработанный нами для специального класса нормальных матриц, применим, в действительности, к более широкому множеству нормальных и даже анормальных матриц. Чтобы описать это множество, напомним, что всякая квадратная комплексная матрица А может быть единственным образом представлена в виде
А = В + iC , (0.29)
где
В = В\ С = С . (0.30)
Эрмитовы матрицы
В=1-{АЛ-А ), С = (А-А ) (0.31)
называются соответственно вещественной и мнимой частями матрицы А. Нормальные матрицы А могут быть охарактеризованы тем свойством, что соответствующие им матрицы (0.31) перестановочны.
В разделе 3.1 рассматривается класс нормальных матриц А, выделяемый условием
к = rank С n , (0.32)
где п — порядок А. Спектр такой матрицы А принадлежит объединению вещественной оси и, самое большее, к прямых, параллельных мнимой оси, т.е. вырожденной алгебраической кривой порядка не выше к + 1. Следовательно, метод MINRES-N применим к матрицам этого типа. Мы показываем, что для матриц со свойством (0.32) MINRES-N имеет следующую привлекательную особенность: начиная с некоторого шага, рекурсия длины 3(к + 1), обычно выполняемая этим методом, может быть заменена трехчленной рекурсией, после чего метод фактически совпадает с эрмитовым процессом Ланцоша. Достигаемое этим путем ускорение сходимости иллюстрируется несколькими примерами.
В разделе 3.2 рассматриваются системы с матрицами (0.29), для которых
rank С = 1 . (0.33)
В отличие от раздела 3.1, не требуется, чтобы А была нормальной матрицей, т.е. ее вещественная и мнимая части В и С не обязаны коммутировать. Мы показываем, что к системам этого типа применим метод MINRES-N2.
Эффективность метода иллюстрируется несколькими примерами решения ленточных систем; на этих примерах производительность MINRES-N2 сопоставляется с производительностью библиотечной программы Matlab a, реализующей GMRES.
В разделе 3.3 условие (0.33) заменяется на
rank С = к 1 . (0.34)
От матрицы А по-прежнему не требуется нормальности. Доказано, что матрица, удовлетворяющая этому условию, может быть приведена посредством унитарного подобия к блочно-трехдиагональному виду, где порядок каждого диагонального блока не превосходит 1+А;. Это означает, что систему (0.2) с такой матрицей А можно решать методом MINRES-N(1 + к), т.е. вариантом метода MINRES-N, первоначально предназначенным для нормальных матриц, которые удовлетворяют уравнению (0.25) для некоторого многочлена степени 1 + к. Обсуждаются примеры, в которых производительность MINRES-N(1 + к) для систем, подчиняющихся условию (0.34), сравнивалась с производительностью GMRES.
Изложение в третьей главе диссертации основано на наших работах [15, 16]. Классы матриц, рассматриваемые в этой главе, играют по отношению к обобщенным крыловским подпространствам роль, сходную с ролью обобщенных (, А;)-нормальных матриц в стандартной теории. В частности, эти классы включают в себя малоранговые возмущения (, &)-нормальных матриц.
Для одного и того же типа (, &)-нормальных матриц длина рекурсии в MINRES несколько больше, чем длина соответствующей рекурсии вида (0.14). Так, для унитарной матрицы U правая часть (0.14) содержит только три слагаемых, поскольку
(s, t) = (€, к) = (0,1) . Максимальная длина рекурсии для той же матрицы в MINRES-N2 равна 6. Однако, в целом, различия в длинах рекурсий не столь существенны, чтобы повлечь значительные различия в общем времени решения системы. Гораздо важнее различия в скорости сходимости, вызываемые различием аппроксимациопных свойств стандартных и обобщенных крыловских подпространств. Но в этом отношении никаких однозначных выводов сделать нельзя: для одних типов матриц быстрее сходятся приближения, выбираемые из подпространств )Ст, для других — приближения, порождаемые обобщенными крыловскими подпространствами.
Что касается сравнения GMRES с MINRES-N, то здесь возможно более определенное заключение: даже в тех случаях, когда число итераций в GMRES значительно меньше, чем в MINRES-N, общее время решения системы в первом методе больше (часто гораздо больше), чем время решения той же системы вторым методом. В этом проявляется преимущество короткой фиксированной рекурсии по сравнению с линейно растущей длиной рекурсии в GMRES.
Отметим, что все численные эксперименты, обсуждаемые в главах 2 и 3, выполнялись на персональном компьютере Pentium IV с тактовой частотой 256 Мгц.
В разделе 4.2 четвертой главы рассматривается и решается следующая задача, мотивированная теоремой 0.4 и теоремами из главы 3: известно, что матрица А Є МП(С) является одноранговой коррекцией эрмитовой матрицы. Иначе говоря, А есть матрица вида
A = H + R, (0.35)
где
Н = Н , rank R = 1 . (0.36)
Однако сами матрицы Н и R неизвестны. Спрашивается, как найти эти матрицы?
В разделе 4.3 рассматривается и решается аналогичная задача, в которой эрмитова матрица Н заменена унитарной матрицей U. В основе изложения — наша статья [17].
В разделе 4.4 мы решаем аналогичную задачу по отношению к разложению
А = X + Y , (0.37)
где
X = Хт, rank Y = 1 . (0.38)
Это решение основывается на нашей публикации [18], а сама задача мотивируется желанием распространить метод типа сопряженных градиентов, построенный для комплексных симметричных матриц в [19], на одноранговые возмущения таких матриц. В заключительном разделе 4.5 обрисован путь, на котором можно получить желаемое обобщение. Вкратце, он состоит в замене унитарных подобий унитарными конгруэпциями или, в терминах обобщенного процесса Ланцоша, замене сопряженной матрицы А транспонированной матрицей Ат.
Обобщенный метод минимальных невязок
К системе (0.23) можно применить стандартный метод сопряженных градиентов. Явное формирование матрицы А А не обязательно, поскольку произведения вида А Ар можно вычислять путем последовательного умножения вектора р на А и А .
Использование, наряду с подпространствами JCm(A,b), крыловских подпространств, порожденных сопряженной матрицей А , и построение пар биортогональных базисов в этих двух последовательностях подпространств. На этих принципах построены такие алгоритмы, как метод би-сопряженных градиентов или метод квазиминимальных невязок.
Недостатком первого приема является то обстоятельство, что число обусловленности матрицы А А есть квадрат числа обусловленности матрицы А. Поэтому переход к системе (0.23) не рекомендуется при решении плохо обусловленных систем. Недостатки второго приема связаны с часто наблюдаемой плохой обусловленностью биортогональных базисов и даже иногда их отсутствием (явление, называемое обрывом процесса ). Общий недостаток состоит в необходимости использования матрично-векторных произведений вида A q наряду с произведениями вида Ар. (Справедливости ради, следует отметить, что существуют варианты этих методов, такие, как алгоритмы CGS или Bi-CGStab, не требующие привлечения матрицы А ).
Результаты настоящей диссертации можно рассматривать как вклад в развитие теории экономичных итерационных методов для решения несамосопряженных систем линейных уравнений. Оригинальность используемого нами подхода состоит в том, что вместо стандартных крыловских подпространств принцип минимальных невязок применяется к обобщенным кры-ловским подпространствам. Эти последние были введены в другом контексте в статье X. Д. Икрамова и Л. Эльзнера ([11]; см. также [12]). Если рассмотреть обобщенную степенную последовательность то линейные оболочки ее конечны отрезков и называются обобщенными подпространствами Крылова, порожденными матрицей А и вектором х.
Мотивом для введения обобщенных крыловских подпространств в [11] было желание найти компактную форму, по возможности близкую к ленточной, для произвольной нормальной матрицы. Как известно, всякую эрмитову матрицу можно привести к трехдиагональному виду посредством конечного процесса, основанного на унитарных подобиях. Для такого приведения могут использоваться различные методы; к технике, используемой в настоящей работе, ближе всего находится метод Ланцоша, тоже относящийся к числу методов (стандартных) крыловских подпространств.
Если матрица А неэрмитова, то она может быть приведена к форме Хес-сенберга посредством метода Арнольди, представляющего собой обобщение процесса Ланцоша. (Этот метод уже упоминался выше в связи с алгоритмом GMRES.) Метод Арнольди не способен извлечь какие-либо выгоды из информации о нормальности матрицы А и дает для нее все ту же заполненную хессенбергову матрицу. Чтобы приблизить компактную форму к ленточной матрице, т.е. но возможности симметризовать ее, нужно, чтобы процесс приведения использовал не только Л, но и сопряженную матрицу А . Это рассуждение и приводит к обобщенным крыловским подпространствам. Более подробное их обсуждение дано в разделе 1.2.
В [11] показано, что компактная форма Я, в которую обобщенный процесс Ланцоша преобразует произвольную нормальную п х n-матрицу А, представляет собой ленточную матрицу с переменной шириной ленты (такие матрицы иногда называют профильными). Матрицу Н можно рассматривать и как блочно-трехдиагональную, причем порядок ее диагонального блока растет с ростом его индексов. Даже в наихудшем случае (т.е. тогда, когда рост порядка происходит скорее всего) общее число ненулевых элементов в Н не превосходит 3\/2п3 2. Эту величину нужно сопоставить с « у ненулевыми элементами хессенберговой матрицы.
Еще более важен следующий факт, также обнаруженный в [И]: пусть нормальная матрица А удовлетворяет соотношению где fix,у) — многочлен степени к. Тогда ширина ленты в матрице Н стабилизируется, начиная с некоторого индекса (зависящего от к). Таким образом, в этом случае компактная форма Н превращается в полноценную ленточную матрицу (или блочно-трехдиагональную матрицу, порядки диагональных блоков которой, начиная с некоторого момента, не возрастают).
Метод MINRES-N, составляющий основное содержание настоящей диссертации, предназначен для класса нормальных матриц, охватываемых (при различных к) соотношением (0.25). Этот класс включает в себя, в частности, эрмитовы матрицы (поскольку Н = Н есть уравнение (0.25) с к = 1), линейные многочлены от них, унитарные матрицы U и матрицы вида A = aU + (ЗІ (см. (0.16) и (0.17)), т.е. ио-существу, весь класс 5-нормальных матриц. При этом он гораздо шире этого последнего класса.
MINRES-N использует ленточность компактной формы Н для построения ортогональных базисов обобщенных крыловских подпространств посредством рекурсий фиксированной длины. К этим подпространствам применяется принцип минимальных невязок, причем условия для его применения гораздо более комфортные, чем в GMRES. Помимо уже упомянутой короткой рекурсии для генерирования векторов спуска (что является наиболее важным обстоятельством), задача наименьших квадратов, решаемая на каждом шаге, имеет ленточную, а не хессенбергову матрицу.
Детали реализации метода
Получаемые в результате векторы Ланцоша дают ортонормированные базисы обобщенных крыловских подпространств т, подобно тому как век торы Арнольди образуют ортонормированные базисы стандартных кры-ловских подпространств К,т. Числовые коэффициенты процесса ортогона-лизации составляют компактную форму Н матрицы А.
Как мы видели в предыдущем разделе, эта форма представляет собой блочно-трехдиагональную матрицу. Для нормальной матрицы А общего вида порядок диагонального блока Нтт растет как т. Однако для матрицы А, удовлетворяющей уравнению (1.13) с многочленом / степени ;, порядки блоков стабилизируются, как только достигнут к. Вариант метода MINRES-N, рассчитанный на такие матрицы, называется MINRES-N&.
Отыскание наилучшего приближения хт в текущем подпространстве Ст. Оптимальность хт понимается в смысле принципа минимальных невязок, который, как и в случае GMRES, приводит к необходимости решать линейную задачу наименьших квадратов с числом неизвестных, равным размерности т подпространства Ст. В отличие от GMRES, матрица этой задачи имеет блочно-трехдиагональную (или, если угодно, ленточную), а не хессенбергову форму.
Более подробное обсуждение отдельных этапов метода будет дано следующей главе для наиболее важного его варианта — метода MINRES-N2. Из этого обсуждения будет ясна конструкция MINRES-N& для произвольного к 1.
Предположим, что нормальная матрица А не является эрмитовой и удовлетворяет уравнению вида (2.1) с ненулевым коэффициентом с22. Умножая (2.1) на вектор Ь, заключаем, что (А )2Ь есть линейная комбинация Л2Ь, АА Ь и векторов Ь, АЪ и А Ь, принадлежащих нулевому и первому слоям последовательности (1.15). Таким образом, чтобы получить базис подпространства 2, достаточно сохранить во втором слое векторы А2Ь и АА Ь, удалив из него (А )2Ъ. Заметим, что оставленные векторы получаются путем умножения матрицы А на векторы первого слоя АЬ и А Ь. Эти последние, в общем случае, линейно независимы. Итак, в рассматриваемой ситуации имеем ш0 = 1, ш\ = 2 и ы2 2; в типичном случае, ы2 — 2. Согласно теореме 1 из [11], Шк 2 для всех /г 2. Более того, нетрудно показать, что векторы Ланцоша каждого слоя к (к 2) можно получить, используя произведения векторов предыдущего слоя только с матрицей А. Произведения вида A q не нужны, начиная с четвертого шага алгоритма MINRES-N2.
Объединяя высказанные соображения, мы получаем алгоритм приведения допустимой нормальной матрицы Л к ее компактной форме Н. Запись алгоритма см. ниже. В ней 6 — заданный ненулевой вектор, в роли которого может выступать правая часть системы (0.2). Поскольку, начиная с к = 1, размерность подпространства & равна нечетному числу 2/с Н-1, алгоритм оформлен как процедура построения заданного числа 2т+1 первых столбцов матрицы Н и такого же числа ортонормированных векторов Ланцоша, которые здесь обозначаются через q\,...,#2m+i Мы переходим к обсуждению второй компоненты метода MINRES -отысканию наилучшего приближения в текущем подпространстве Се- Это обсуждение облегчается подробным разбором технологии GMRES, проведенным в разделе 1.1. Символы qi и Я будут теперь относиться к ортонормированным векторам Ланцоша и компактной форме Я, которые строятся приведенным выше алгоритмом. Если, как и прежде, Q — унитарная матрица, составленная но столбцам из векторов qi,... ,qn, то продолжает действовать равенство (1.1). Предположим, что выполнены га шагов алгоритма приведения к компактной форме (при этом первым шагом считается нормировка q\ — 6/Ц6Ц2). Число т будем считать нечетным: т = 21 + 1; тогда векторы 1,..., /т образуют базис обобщенного крыловского подпространства Се(А,Ь). Примем разбиение (1.4) для матрицы Н. Часть этого разбиения, а именно (2.4) соответствует первым га столбцам матрицы, изображенной на рис 2.1. Поиск приближения хт как вектора из подпространства Се, минимизирующего 2-норму невязки, с помощью разложения (1.2) сводится к поиску m-мерного вектора ут (см. (1.3)). Как и в разделе 2.1.1, вычисление ут равносильно решению задачи наименьших квадратов (1.5). Учитывая форму матрицы (2.4) (см. рис 2.1), это решение можно выполнить следующим образом: 1. Посредством вращения Ru аннулировать элемент (2,1) матрицы (2.4). 2. Для j = 2,..., га аннулировать поддиагональные элементы столбца j матрицы (2.4) посредством подходящего отражения Pj. 3. Применить те же преобразования, т.е. вращение Ru и отражения 2, Рт, к правой части задачи (1.5).
Анормальные возмущения большего ранга
Из (2.7) следует, что d Є R и хотя бы один из коэффициентов с и d ненулевой. Если, в частности, с ф 0, то выполнено условие (2.2). Таким образом, если спектр нормальной матрицы А лежит на любой кривой второго порядка, кроме окружности, то система с этой матрицей может быть решена методом MINRES-N2, реализованным так, как описано в данном разделе. Случай окружностей будет рассмотрен в следующем разделе.
Ниже мы описываем результаты численных экспериментов, в которых методы MINRES-N2 и GMRES сравнивались на системах линейных уравнений порядка 2000 с нормальными матрицами коэффициентов и спектром, лежащим на некоторой кривой второго порядка. Обусловленность систем, как правило, была очень хорошей. Метод GMRES был представлен библиотечной процедурой gmres системы Matlab. Для MINRES-N2 была составлена специальная Matlab-процедура, текст которой составляет Приложение 1.
Чтобы обеспечить принадлежность спектра желаемой кривой Г, мы выбирали А как диагональную матрицу, диагональные элементы которой суть точки, лежащие на Г. При исследовании сходимости метода такой выбор А не влечет какой-либо потери общности. Действительно, в некотором орто-нормировашюм базисе {е} пространства Сп нормальная матрица А приводится к диагональному виду D, а исходный итерационный процесс изомет-ричен такому же процессу, проводимому в базисе {е} с матрицей D. В то же время, нужно признать, что сравнение методов, проводимое лишь для диагональных систем, не дает вполне объективной картины. Для большой системы, в зависимости от степени заполненности матрицы А, время формирования матрично-векторных произведений вида Aq может составлять значительную или даже преобладающую часть общего времени решения системы.
Для выхода из итерационного процесса использовалось условие где є — заданное положительное число; обычно мы полагали є = Ю-8. Величина гш2 в алгоритме MINRES-N2 вычисляется по формуле (2.5). Пример 1. Спектр матрицы А системы (0.2) равномерно распределен на эллипсе Эта система неожиданно оказалась очень трудной для метода GMRES: при є = Ю-6 потребовалось 690 итераций, чтобы удовлетворить условие (2.9). В методе MINRES-N2 выход произошел после 34 итераций. В том, что число шагов в MINRES-N2 меньше, чем в алгоритме минимальных итераций GMRES, нет противоречия, поскольку методы используют крыловские подпространства различного типа. Время вычисления решения составило i = 0.28 сек для MINRES-N2 и t2 = 175 сек для GMRES. Известно, что ситуация, когда спектр системы целиком расположен в некоторой полуплоскости, граница которой проходит через нуль, благоприятна для GMRES. Это подтвердилось в нашем примере. Если спектр матрицы А равномерно распределен на правой половине эллипса (2.10), то при є = Ю-8 GMRES сходится за 64 итерации и требует ti — 22.66 сек. MINRES-N2 также ускоряется и сходится за 28 итераций {i\ — 0.25 сек). Пусть теперь спектр матрицы А равномерно распределен на верхней половине эллипса (2.10). GMRES по-прежнему сходится быстро (є = 10 8, 48 итераций, 2 = 18.8 сек), тогда как сходимость MINRES-N2 резко замедляется (150 итераций). Однако и в этом случае t\ = 2.7 сек t2, поскольку шаги MINRES-N2 требуют меньше времени, чем шаги GMRES. Предположим, наконец, что спектр матрицы А равномерно распределен на четверти эллипса (2.10), расположенной в первом квадранте. В этом случае очень быстро сходятся оба метода: MINRES-N2 требует 18 итераций и времени t\ — 0.16 сек, a GMRES — 23 итерации и г — 16.6 сек. Поведение обоих методов очень похоже на то, что мы имели в случае полного эллипса: для є = 10_6 GMRES потребовал 203 итерации и времени t2 = 38.2 сек, тогда как в MINRES-N2 были выполнены 66 итераций и затрачено время t\ = 0.92 сек. Если для тех же значений х сосредоточить спектр на верхних иолу-ветвях гиперболы (2.11), то оба метода улучшают свои показатели: при є = Ю-8 GMRES требует 59 итераций и времени t2 = 21.36 сек, a MINRES-N2 соответственно 38 итераций и t\ — 0.36 сек. Улучшение сходимости наблюдается и для спектра, сосредоточенного на одной ветви гиперболы, и тем более для случая, когда спектр принадлежит одной полуветви. Пример 3. Для значений ж, принадлежащих сегменту [0.1, 1.1], спектр матрицы А равномерно распределен на нижней и верхней дугах параболыВ этом примере GMRES значительно превосходит MINRES-N2 с точки зрения числа итераций: при є = Ю-8 первому методу достаточно 79 итераций, тогда как второму нужны 300 итераций. Тем не менее, и здесь MINRES-N2 имеет значительное преимущество в скорости: t\ — 8.84 сек, t i — 23.2 сек. Сосредоточим теперь спектр матрицы А на верхней дуге параболы (2.12), отвечающей отрезку [0.03, 1.2] на оси х. Обусловленность системы ухудшится, но не слишком сильно. Принадлежность спектра квадранту оказывается более важным фактором, поэтому оба метода улучшают свои характеристики. Однако качественно картина сходимости остается прежней: GMRES требует 45 итераций и времени ti — 18.17 сек, a MINRES-N2 при 110 итерациях затрачивает лишь t\ — 2.11 сек.
Восстановление комплексных симметричных матриц
Доказанная импликация вместе с оценкой (3.14) означает, что строгое возрастание размерностей в последовательности (3.15) возможно лишь для первых к значений индекса г. Для т, бблыыих, чем к, приращение размерности подпространства Ст по сравнению с
Cm-i возможно лишь за счет матрицы Бт. Это доказывает утверждение (3.12). Теорему 3.1 можно использовать для ускорения сходимости метода MINRES-N. Объясним это с помощью следующего примера.
Пусть А — диагональная матрица порядка 2000, все диагональные элементы которой, за исключением четырех, расположены в двух симметричных сегментах вещественной оси [-14,-10] и [10,14]. Остальные четыре элемента образуют две сопряженные пары « ±0.1г и « ±0.5г. Итак, спектр сосредоточен на координатных осях, т.е. принадлежит (вырожденной) кривой второго порядка
Поэтому к системе с матрицей А можно применить метод MINRES-N2. То обстоятельство, что диагональная система уравнений решается тривиально, не должно вызывать недоразумений. Действительно, всякая нормальная матрица А приводится к диагональному виду D в некотором ортонормированием базисе {е}, а метод MINRES-N2, проводимый с матрицей Д изометричен такому же процессу, применяемому к D в базисе {е}. Поэтому, с точки зрения исследования сходимости метода, ограничение диагональными матрицами не влечет за собой никакой потери общности. Правую часть Ь системы (0.2) мы выбирали в виде Таким образом, точным решением системы является вектор, составленный из единиц. На рис. 3.1 показаны графики евклидовых длин невязок для методов MINRES-N2 и GMRES, применяемых к описанной выше системе. В этих расчетах GMRES был представлен библиотечной функцией gmres системы Matlab, a MINRES-N2 — составленной нами Matlab — процедурой. Для обоих методов использовался один и тот же критерий остановки: где і— текущая невязка, а положительное число є задается пользователем. Для обсуждаемого примера є = Ю-8.
Результаты сравнения двух методов для данной системы можно считать благоприятными для MINRES-N2: при приблизительно одинаковом числе итераций (50 против 47 итераций для GMRES) выигрыш в общем времени счета происходит вследствие меньшего количества выполняемой арифметической работы. Однако обращает на себя внимание тот факт, что график г для MINRES-N2 состоит из ступенек ширины 4. Это означает, что после спуска на очередную ступеньку три последующих итерации не уменьшают норму невязки.
Объяснение состоит в следующем: MINRES-N2 исходит из предположения, что все числа шт равны двум. Между тем, согласно теореме 3.1, для рассматриваемой системы шт 1 при т 4. Пусть #1,(/2,--- — векторы ортонормированной последовательности, которую строит MINRES-N2, т.е. обобщенные векторы Ланцоша. При этом в качестве q\ берется нормированный вектор правой части 6, а каждая последующая пара векторов Ланцоша теоретически получается посредством ортопормализации, применяемой к векторам очередного слоя последовательности (1.15). В частности, векторы #1,... ,#9 образуют ортонормированный базис подпространства 4- После того как на 10-м шаге процесса построен вектор qw, MINRES-N2 пытается вычислить вектор 7ц, оставаясь в том же (пятом) слое последовательности (1.15). Однако о 5 = 1, т.е. в пятом слое нет вектора, независимого с qw и ортогонального к С\. В методе MINRES-N2 это обстоятельство проявляется в том, что вектор z, полученный из второго вектора пятого слоя ортого-нализацией к предыдущим векторам Ланцоша, имеет малые компоненты, фактически являющиеся накопленными ошибками округлений. Несмотря на это, MINRES-N2 нормирует z и принимает его в качестве q\\. Естественно, что такой вектор, случайный по отношению к подпространствам Ст, мало помогает в уменьшении длины невязки. Описанная ситуация повторяется на каждом последующем слое.
выполняемую перед нормированием, дающим очередной вектор Ланцоша. Число 5 задается пользователем; в описываемых экспериментах мы полагали S = 10 4. Если неравенство (3.19) выполнено, то текущий шаг завершается стандартным образом; в противном случае, уменьшаем значение и ширины текущего слоя на единицу, отбрасываем вектор z и переходим к вычислениям, соответствующим следующему шагу процесса ортогонализа-ции.
Рис. 3.2 показывает поведение норм невязок для GMRES и модифицированного алгоритма MINRES-N2. Эти графики соответствуют системе, построенной, как и выше: матрица А — диагональная со спектром, сосредоточенным на сегментах [-14,-10] и [10,14], кроме двух сопряженных пар « ±0.5г и « ±0.9г. Некоторое отличие этой матрицы от предыдущей связано с тем, что при заданном характере спектра (заданное число чисто мнимых сопряженных пар + равномерное распределение прочих собственных значений на заданных сегментах вещественной оси) конкретные значения диагональных элементов выбирались с помощью датчика псевдослучайных чисел.
Мы видим, что картина сходимости для GMRES почти полностью сохранилась; в то же время, модифицированный MINRES-N2 сходится значительно быстрее первоначального метода: критерий (3.18) теперь выполняется (при є = 10 8) после 30 итераций. Время счета в три раза меньше, чем для GMRES; ступеньки ширины четыре сменились ступеньками ширины два.
Почему ступеньки лишь уменьшили ширину, а не исчезли вовсе, автору в настоящее время не вполне ясно. Возможно, это связано с симметричным расположением спектра рассмотренных выше матриц относительно нуля. Во всяком случае, ступенек нет (после девятого шага ортогонализации) на рис. 3.3. Показанные здесь графики норм невязок соответствуют матрице А с двумя сопряженными парами собственных значений « ±0.2г и w ±0.4г и остальным спектром, равномерно распределенным на сегменте [1,4]. Для выполнения критерия (3.18) с тем же є = Ю-8 здесь достаточно 24 итераций. Заметим, что и GMRES сходится для этой системы быстрее (32 итерации), чем в предыдущих случаях.