ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2008, том 48, < 12, с. 2092-2106
УДК 519.612
ГОМОТОПИЧЕСКИЙ МЕТОД ДЛЯ РЕШЕНИЯ СИММЕТРИЧНЫХ
ТЁПЛИЦЕВЫХ СИСТЕМ
© 2008 г. М. Ван Барел*, X. Д. Икрамов**, Ä. Ä. Чесноков*,**
(*Katholieke Universiteit Leuven, Department of Computer Science, Celestijnenlann 200 A, B-3001 Leuven, Belgium;
**119899 Москва, Ленинские горы, МГУ, ВМК)
e-mail: marc.vanbarel@cs.kuleuven.be; ikramov@cs.msu.su; ache@mccme.ru
Поступила в редакцию 29.12.2007 г. Переработанный вариант 22.05.2008 г.
Предлагается быстрый алгоритм решения симметричных тёплицевых систем. Метод постепенно и непрерывно преобразует единичную матрицу в обратную к заданной тёплицевой матрице. Требования к оперативной памяти составляют O(n), сложность метода O(log к( T)nlog n), где к(Т) равно числу обусловленности матрицы T. Приведены результаты численных экспериментов, подтверждающие эффективность метода. Библ. 29. Фиг. 3. Табл. 3.
Ключевые слова: тёплицевы матрицы, циркулянты, сверхбыстрый алгоритм, метод продолжения по параметру, итерационное уточнение, собственные значения.
1. ВВЕДЕНИЕ Рассмотрим линейную систему уравнений
Tx = b (1)
с невырожденной симметричной n х n-тёплицевой матрицей коэффициентов T и правой частью b е [n (или Cn). В настоящей статье представлен новый итерационный алгоритм решения такой системы со сложностью O( log к( T) n log n), где к(Т) обозначает число обусловленности T; требования к памяти - порядка O(n).
Матрица T е [n х n называется тёплицевой, если T = (ц) = (a;-_,):
T =
ao ax a2 .. an -
a-1 ao a1 •
a_2 a-1 ao • •• a2
a1
'-(n -1) a-2 a- -1 ao
(2)
Тёплицевы матрицы возникают в разных задачах прикладной математики, таких как аппроксимация дифференциальных уравнений (см. [1], [2]), аппроксимации Паде (см. [3]-[5]), локализация корней многочленов (см. [6]-[8]).
Тёплицева матрица определяется только 2п - 1 параметрами. Этот факт послужил основой для создания многих быстрых и супербыстрых алгоритмов для решения тёплицевых систем (1). Такие алгоритмы используют структуру матрицы. Два типа прямых быстрых методов со сложностью 0(п2) арифметических операций представляют собой методы типа Левинсона и типа Шура. Ссылки и другая информация относительно этих методов приведены в [9].
Как правило, такие методы требуют численной невырожденности ведущих главных подматриц в Т. Более того, если эти подматрицы являются невырожденными, но плохо обусловленными, то в случае реализации в арифметике с конечной точностью может возникать численная неустойчивость.
2092
Алгоритмы со сложностью менее 0(п2) называют супербыстрыми. Впервые идея супербыстрого метода для тёплицевых систем была предложена в [10]. Супербыстрые алгоритмы были разработаны в [11], [4], [12], [13]. Некоторые современные алгоритмы представлены в [14]-[17].
Один из лучших количественных способов для определения и использования структуры матрицы в вычислениях основан на концепции рангов смещения. Впервые это понятие появилось в [18]. Говорят, что матрица М имеет малый ранг смещения, если для некоторых матриц А и В ранг матрицы X = М - АМВ мал. При подходящем выборе матриц А и В тёплицева матрица Т может быть превращена в матрицу X ранга два. После этого X можно представить в виде произведения
X = ОН*, где О, Н е 1Кп х (или О, Н е Сп Х ). Матрицы О и Н называются генераторами смещения для Т. На основе О и Н появляется возможность вычислить генераторы для Т-1. Их размер также будет составлять п х 2.
При некоторых условиях на матрицы А и В можно вычислить произведение тёплицевой матрицы Т и обратной к ней Т-1 на вектор, используя только их генераторы смещения, за 0( п^ п) арифметических операций.
Таким образом, мы осуществляем все операции с тёплицевыми матрицами и обратными к ним с использованием концепции рангов смещения. Гомотопический метод, предложенный в настоящей статье, непрерывно преобразует генераторы смещения для обратной к некоторой (легко обратимой) тёплицевой матрице в генераторы смещения для обратной к заданной тёплицевой матрице Т. Для улучшения грубых приближений на каждом шаге гомотопического метода мы используем итерационное уточнение.
Итерационное уточнение сохраняет ранги смещения тех матриц, к которым оно применяется. Это дает возможность отказаться от применения алгоритмов сжатия, необходимых в других методах итерационного улучшения (см., например, [19]), что приводит к лучшей асимптотической оценке сложности (при сравнении, например, с недавним алгоритмом из [20]). Другим преимуществом настоящего метода является строгая верхняя оценка количества используемой памяти, так как для методик сжатия обычно требуется дополнительная оперативная память.
Настоящая статья имеет следующую структуру. В разд. 2 приведены различные процедуры итерационного улучшения. Общий гомотопический метод описан в разд. 3. Необходимые утверждения относительно рангов смещения представлены в разд. 4. Окончательная версия алгоритма приведена в разд. 5. Вопросы сходимости и сложности изучаются в разд. 6. Компьютерная реализация метода и результаты вычислительных экспериментов представлены в разд. 7.
2. ПРОЦЕДУРЫ ИТЕРАЦИОННОГО УЛУЧШЕНИЯ
Грубое начальное приближение X0 к обратной М-1 для невырожденной матрицы М может быть улучшено посредством различных процедур итерационного улучшения (в дальнейшем обозначаются как 11Р).
Такие процедуры включают в себя итерации по Ньютону (см. [21])
Xi + 1 = 2Xi - XiMXi, I = 0, 1,..., (3)
итерации по Ньютону с масштабированием (см. [22])
X¡ + 1 = а + 1(2X¡ - XiMXi), I = 0, 1,., (4)
формулу для уточнения на основе кубического полинома (см. [19])
Х{ +1 = aXi(MXi)2 + bXi(MXi) + cXi + М, г = 0, 1, ..., (5)
итерационное уточнение (см. [23]) для уравнения MX = I
Я, = I - MXi, Д = X 0 Я, X; +1 = X¡ + А, I = 0, 1, .... (6)
Для всех приведенных выше процедур обозначим через я, невязку I - MXi. Тогда их сходимость определяется 2-нормой начальной невязки ||Я01| = || I - MX01| (см. [24], [19], [23], [20]):
¡Я! ||Я0||Р , Р зависит от процедуры. (7)
Предположим, что известна матрица М0, обратная к ней М01 и некоторая матрица М. Назовем М(/ начальным приближением к М-1. Попробуем улучшить это приближение с использованием одной из процедур (3)-(6). Однако такая процедура может быть расходящейся, если ||Я01| 5г 1 или, другими словами, если начальное приближение является слишком грубым (расстояние от матрицы М0 до матрицы М слишком велико). Один из возможных способов решения этой проблемы представлен в следующем разделе.
3. ГОМОТОПИЧЕСКИЙ МЕТОД
В гомотопических методах выбираются простая задача, решение которой известно, и траектория между этой простой задачей и сложной задачей, которую на самом деле нужно решить. Решение простой задачи непрерывно преобразуется в решение сложной задачи при движении по выбранной траектории. Траектория может быть определена посредством введения в задачу дополнительного скалярного параметра.
Предположим, что нужно решить задачу F(x) = 0, x е в которой - некоторое линейное пространство. Выберем некоторую точку x0 е и определим функцию продолжения
H(x, а) = F(x) - (1- а)F(x0).
Будем затем решать задачу
H(x, а) = 0 (8)
для каждого значения а между 0 и 1. Когда а = 0, решением задачи (9) является в точности x = x0. Когда а = 1, имеем, что H(x, 1) = F(x) и, таким образом, решение задачи (8) совпадает с решением исходной задачи F(x) = 0. Существуют и другие способы определения траектории, более тесно связанные со структурой исходной задачи.
В настоящее время гомотопические методы имеют ряд приложений, таких как обобщенная проблема собственных значений (см. [25]), решение систем полиномиальных уравнений (см. [26], [27]), задачи многоцелевой оптимизации (см. [28]).
Теперь мы можем кратко сформулировать гомотопический метод в применении к задаче обращения матриц. Обозначим через А0 некоторую легко обратимую матрицу, через А1 - матрицу, обратную к которой нужно вычислить, и построим матричную функцию продолжения
М (а) = (1- а) А0 + а А1.
Предположим, что мы выбрали одну из процедур IIP (3)-(6). Тогда сущность гомотопического метода может быть описана следующим образом.
Алгоритм 1. Первый вариант гомотопического метода
вход : А0 - начальная матрица, А1 - заданная матрица
. -1
выход : А1
обозначения : М(а) = (1 - а)А0 + аА1 begin
а = 0; вычислить обратную Aj,1 = М(0)-1 while а < 1 do
увеличить а так, чтобы а =s 1 и выбранная процедура IIP сходилась улучшить приближение к М(а)-1 с использованием обратной матрицы, полученной на предыдущем шаге в качестве начального приближения для IIP end
return М(1)-1 = А-1 end
Таблица 1
Операторные матрицы Класс структурированных Ранг
Л В матриц М Ул, в(М)
¿0 тёплицевы и обратные к ним <2
Х 0 ганкелевы и обратные к ним <2
£(0 ¿0 вандермондовы <1
^0 £(0 обратные к вандермондовым <1
4. МАТРИЦЫ МАЛОГО СТРУКТУРНОГО РАНГА
Тёплицевы матрицы являются важным типом структурированных матриц. Вычисления со структурированными матрицами, проводимые надлежащим образом, имеют меньшие требования к памяти и меньшую вычислительную сложность по сравнению с аналогичными действиями для матриц общего вида.
4.1. Общая теория
Определим формально на множестве вещественных или комплексных т х п-матриц линей-
ные операторы Ь :
и типа Стейна (Ь = ДА, в)
типа Сильвестра (Ь = УА, в) Ь (М) = УА, в (М) = АМ - МВ
(9)
(10)
Ь (М) = ДА, в (М) = М - АМВ
для некоторой пары {А, В} операторных матриц.
Операторы типа Сильвестра и типа Стейна взаимосвязаны, и для них справедливо Предложение 1 (см. [24]). Пусть УА, в и ДА, в - операторы смещения, определенные в (9) и (10) соответственно. Тогда между ними существует следующая связь: УА, в = А Д^-1 в, если операторная матрица А не вырождена, и УА, в = - ДА -1 в, если операторная матрица в не вырождена.
Часто используются операторные матрицы Х/, 2* и В^). Здесь Х/ - простейший /-циркулянт (/ - некоторое число):
2/ =
0 1
/
1 0
и В(ч) = (V)я==о.
Зафиксируем матрицы А и в. Пусть Ь = УА, в или Ь = ДА, в. Ранг г матрицы смещения Ь(М) называется рангом смещения матрицы М. Если А и в выбраны
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.