Квазисепарабельные матрицы в линейной алгебре и ее приложениях тема диссертации и автореферата по ВАК РФ 01.01.07, кандидат физико-математических наук Жлобич, Павел Георгиевич

• Жлобич, Павел Георгиевич
• кандидат физико-математических наук
• 2012

Оглавление диссертации кандидат физико-математических наук Жлобич, Павел Георгиевич

Введение

1.1 Квазисепарабельные матрицы.

1.2 Историческая справка

1.3 Основные результаты данной работы

1.4 Содержание работы по

главам

Глава 1. Применение обобщенных матриц Фробениуса для вычисления корней многочленов

1.1 Введение.

1.2 ЫГ-разложение квазисепарабельных матриц.

1.3 дс1-алгоритмы для квазисепарабельных матриц

1.3.1 Стационарный qd-aлгopитм.

1.3.2 Прогрессивный дс1-алгоритм.

1.4 Дифференциальный дс1-алгоритм для квазисепарабельных матриц в хессенберговой форме.

1.4.1 с^ск алгоритм на языке многочленов.

1.5 Обобщенный процесс Грама-Шмидта и с^с^-алгоритм для квазисепарабельных матриц.

1.5.1 Включение сдвигов.

1.6 с^с1Б-алгоритм для обобщенной матрицы Фробениуса

1.6.1 Многочлен, заданный в базисе мономов.

1.6.2 Многочлен, заданный в базисе ортогональных многочленов.

1.6.3 Многочлен, заданный в базисе многочленов Сегё

1.7 Численные эксперименты.

1.8 Выводы.

Глава 2. Обращение полиномиальных матриц Вандермонда

2.1 Введение.

2.2 Принцип Теллегена.

2.3 Минимальная реализация динамической системы через граф потока сигнала.

2.3.1 Многочлены Горнера и обращение матрицы Вандермонда

2.3.2 Связь многочленов Горнера с собственными векторами матрицы Фробениуса.

2.3.3 Одночлены и прямой (наблюдаемый) граф

2.3.4 Многочлены Горнера и дуальный (управляемый) граф.

2.4 Неминимальная реализация динамической системы и обращение полиномиальной матрицы Вандермонда

2.5 Алгоритм обращения полиномиальной матрицы Вандермонда

2.6 Численные эксперименты.

2.7 Выводы.

Глава 3. Квазисепарабельные матрицы в задачах оптимального управления с уравнениями в частных производных

3.1 Введение.

3.2 Пример оптимизационной задачи с ограничениями в виде УЧП.

3.3 Многоуровневые квазисепарабельные матрицы.

3.4 Численные эксперименты.

3.5 Выводы.

Глава 4. Численная устойчивость быстрых алгоритмов решения систем с квазисепарабельными матрицами 93 ^ 4.1 Введение.

4.2 Основные идеи алгоритмов решения квазисепарабель-ных систем.

4.3 Алгоритм А

4.3.1 (^-разложение.

4.3.2 Обратная прогонка.

4.3.3 Анализ ошибок округления в алгоритме А

4.4 Алгоритм Б

4.4.1 Источники неустойчивости алгоритма Б

4.4.2 Численные эксперименты с алгоритмом Б

4.5 Выводы.

Введение диссертации (часть автореферата) на тему «Квазисепарабельные матрицы в линейной алгебре и ее приложениях»

Л. Квазисепарабельные матрицы

Данная работа посвящена развитию теории матриц с малоранговой структурой особого вида, называемых квазисепарабелъными, а также их применению для решения некоторых классических задач линейной алгебры. Прежде чем определить класс квазисепарабелъных матриц, сделаем несколько общих замечаний о матрицах с малоранговой структурой.

Матрицы с блоками малого ранга все активнее применяются в инженерной и научно-исследовательской практике. Например, при решении задач математической физики в интегральной форме, где возникают системы уравнений с плотными матрицами большого размера. Использование малоранговой структуры таких матриц позволяет не только «сэкономить» на хранении их коэффициентов, но и значительно ускорить некоторые операции линейной алгебры с ними (такие, как умножение на вектор). В качестве простейшего примера рассмотрим краевую задачу с оператором Штурма-Аиувилля:

С(и) (р(х)и,)'-ч(х)и = Пх), и^схп^сО + ^иЧа) = 0, (1) в2{и) := а2и(Ъ) + |32и'(Ъ) = 0.

В теории операторов Штурма-Лиувилля доказывается, что любое решение и(х) этой краевой задачи представимо в виде интеграла и(х) где д(х, £,) — так называемая функция Грина для (1). Причем функция Грина, для одномерного оператора Штурма-Лиувилля имеет вид п г /и-1(х)и2(£,), а ^ х ^ £,, где 111,112— частные решения однородного уравнения С{щ) = 0 с ¿^(-щ.) = 0. Представление (2) означает, что функция д(х, £,) — сепара-бельная в двух подобластях своей области определения. Следовательно, ее дискретный аналог — это матрица, в которой любая подматрица выше или ниже диагонали (с захватом диагонали) имеет ранг один.

Приближение больших областей в матрицах, возникающих при рассмотрении нелокальных операторов, матрицами малого ранга активно используется в современных численных методах [1, 2, 3]. Как и в случае задачи Штурма-Лиувилля, матрицы малого ранга связаны с разделением переменных в рассматриваемом операторе: г г f(*,y) = ^uk[x)vk{y) —> F = ^ukvJ. k=1 k=1

Существует значительное отличие одномерной задачи Штурма-Лиувилля от подобных задач в двумерных и трехмерных пространствах. Дело в том, что переменные в функциях Грина для большинства дифференциальных операторов в пространствах размерностей выше первой не разделяются аналитически. Примеры таких функций Грина:

• log |х — у |, 1 /\х — у | для уравнения Пуассона в 2D и 3D соответственно;

• ek'x~y'/lx — yl Для уравнения Гельмгольца в 3D.

Однако в этом случае можно построить разбиение возникающей матрицы на блоки такое, что каждый блок будет аппроксимироваться матрицей малого ранга с высокой точностью. В непрерывном же случае это означает, что первоначальная функция почти сепарабельна в разных подобластях своей области определения. Этод подход успешно реализован в мозаинно-скелетонном методе [1, 4, 5, 6] и методе иерархических матриц [7].

Квазисепарабельные матрицы — более узкий класс, чем мозаично-скелетонные или иерархические, и поэтому их область применения несколько уже. В то же время они обладают более простой и эффективной параметризацией. Хранение матрицы в квазисепарабель-ном формате требует 0(п) ячеек памяти, и алгебраическая сложность большинства алгоритмов с ними также 0[п). Напомним, что для мозаично-скелетонных и иерархических матриц количество памяти и сложность алгоритмов зависят как O(nloga(n) logb(£-1)) от размера п и точности е. Строгое определение квазисепарабельных матриц приведено ниже.

Определение 1. Матрица А называется квазисепарабелъной порядков р и q, если выполнено условие max rank A(k + 1 : n, 1 : k) ^ p, max rank Ail : k, k + 1 : n) ^ q при фиксированных р и q .

Квазисепарабельные матрицы расширяют класс семисепарабель-ных матриц, определенных ниже.

Определение 2. Матрица А называется семисепарабельной порядков р и q, если существуют матрицы И] и такие, что tril(A1) = tril(Ri), rank(Ri) = р, triu(A-1) = triu(R2), rank(R2) = q, где 1г11 и 1;гш обозначают ниэюне- и верхнетреугольные части матриц (включая диагональ) соответственно.

Семисепарабельные матрицы — это матрицы, обратные к ленточным. В определении квазисепарабельных матриц условие на малоранговую часть матрицы ослабляется. семисепарабельная квазисепарабельная

Вместо требования, чтобы вся нижне-/верхнетреугольная часть матрицы была частью некоторой малоранговой матрицы, требуется, чтобы любая подматрица ниже и выше диагонали имела малый ранг. Это условие влечет за собой, каик будет показано далее, замкнутость класса квазисепарабельных матриц относительно операции обращения.

Теорема 3 (Густафсон, 1984). Пусть матрица А и обратная к ней А-1 разделены на блоки таких размеров, что определено их блочное произведение:

АЦ Ai2 , А"1 - В11 В12

А21 А22 В21 В22

Тогда размерности ядер соответствующих подматриц в А и А-1 равны: dim ker Ai i = dim ker В 22, dim ker A22 = dim ker В11, dim ker Au = dim ker В12, dim ker A21 = dim ker B21.

Следствие 4. Пусть A — квазисепарбельная матрица (определение 1). Тогда А-1 также квазисепарабельная, причем порядки квазисепарабельности для А и А-1 равны.

Замкнутость квазисепарабельных матриц относительно обращения позволяет построить алгоритмы их обращения линейной по размеру сложности. Многие другие операции линейной алгебры с квазисепа-рабельными матрицами также имеют линейную сложность. Полный перечень известных на данный момент алгоритмов приведен в таблице 1.

Таблица 1. Сложность алгоритмов с квазисепарабельными матрицами умножение на вектор 0{п) ьи 0(п) 0{п) обращение 0{п)

ЛВ ом

А±В 0{п) вычисление фробениусовой нормы 0(п) перевод поэлементного представления в квазисепарабельное (если последнее существует) 0(п2)

Линейная сложность в вышеприведенных алгоритмах достигается за счёт использования малопараметрических представлений квазисепарабельных матриц. Существует несколько альтернативных способов параметризовать квазисепарабельную матрицу. В данной работе мы будем использовать так называемое генераторное представление [8], приведенное ниже: di giH2 gib2H3 ••• gib2. .bnihn

P2qi d2 д2Нз • • • g2b3 •. bnihn.

Рза2Ц] РзЦ 2 d3 ■■■ g3b4.bniHR

Рп^Тг—1 • • • &2Я1 Vn Q-n—1 • • • &ЗЯ2 Pn ani . a4q3 ••• dn J

Параметры {d^, q^, a^pic, g^, b^, h*} — матрицы, чьи размеры приведены в таблице 2. Максимальный размер генераторов соответствует порядку квазисепарабельности: тах^ r\ = р, тах^ = q .

Благодаря своему простому виду и быстрой арифметике квази-сепарабельные матрицы находят свое применение во многих областях прикладной математики. В частности, в теории систем и теории управления, ортогональных многочленах и теории приближений,

Таблица 2. Размеры генераторов блочных квазисепарабельных матриц. dk qk dk Pk Як bk Hk строк 1 r* r[ 1 1 rg1 столбцов 1 1 rj.^ 1 интегральных уравнениях и статистике. Данная работа посвящена решению некоторых классических задач линейной алгебры, связанных с ортогональными многочленами и теорией оптимального управления. Прежде чем приступить к описанию основных результатов работы, расскажем об истории развития теории квазисепарабельных матриц, этапом которой, мы надеемся, является и данная работа. i.2. Историческая справка

Матрицы с квазисепарабельной структурой были впервые упомянуты в фундаментальной работе Ф. Р. Гантмахера и М. Г. Крейна «Осцилляционные матрицы и малые колебания механических систем» [9]. Продемонстрируем, как квазисепарабельные матрицы возникают в теории осцилляций. Рассмотрим натянутую струну длины I с концами, фиксированными в точках х = 0 и х = I. Обозначим через К(х, s) смещение струны в точке х под действием единичной силы в точке s. Для вывода формулы для К(х, s) воспользуемся рисунком 1. Приложим единичную силу F = 1 к точке 1 и обозначим через Т силу натяжения струны. Тогда

T(sin(a) + sin( (3) ) = 1.

При условии малости углов ос и (3 имеем , л Kfs,s) . K(s.s) sin(ct) - ^ ' \ sm(|3) = , ' s l — s

Следовательно,

Чтобы найти величину смещения К(х, s) в произвольной точке х, воспользуемся тригонометрическим тождеством

К(х, s) K(s,s) X s

Окончательно имеем

4(1-5) Л

К(х,в) — < з(1—х)

71

X ^ Б, Б.

5)

Сравнивая (2) и (5), легко видеть, что К(х,в) есть не что иное, как функция Грина для некоторой краевой задачи Штурма-Аиувилля.

Рис. 1. Вывод функции смещения для натянутой струны. У чЛ"

1 1 ^^ 1 |К(х, в) К(5,8)

1 1 1 х

Рассмотрим конечное число узлов на струне {х^}-^ . Тогда

К(Хг,хЛ ] ЛЬУ) \ ^

6) где гц = и V, I . Очевидно, матрица в (6) — квазисепарабель-ная. К тому же, как показано в работе Ф. Р. Гантмахера и М. Г. Крейна, она — осцилляционная.

Определение 5. Матрица А называется осцилляционной, если все ее миноры всех порядков неотрицательны.

Осцилляционные матрицы чрезвычайно важны при изучении колебаний механических систем. Например, если в узлах {Хг}-^ струны, изображенной на рисунке 1, расположены точечные массы {ггц.}^1, то все собственные частоты колебаний струны ш являются решениями уравнения с1е1;(1 — си2 diag(m1, т2,., тгцОК) = 0.

Ф. Р. Гантмахер и М. Г. Крейн в [9] также впервые доказали, что матрица, обратная к симметричной матрице Якоби (то есть трехдиа-гональной неприводимой), является семисепарабельной (определение

2) порядка 1. Важный шаг в обобщении этого результата был сделан У. Барреттом [10], который опустил условие неприводимости в трех-диагональной матрице и доказал, что в общем случае обратная к произвольной трехдиагональной матрице — квазисепарабельная (определение 1) порядка 1. Напомним, что квазисепарабельные матрицы включают семисепарабельные как свой подкласс. Дальнейшие обобщения на блочный случай и случай произвольных ленточных матриц были получены Ю. Икебе [11], В. Цао и У. Стюартом [12]. Следует также упомянуть достаточно общую теорему, доказанную Э. Асплундом в 1959 году [13].

Теорема 6 (Асплунд, 1959). Пусть матрица А = [а^] — р -хес-сенбергова, то есть а^ = 0, если ) > \ + р. Тогда ранг любой подматрицы матрицы А-1 выше р -ой поддиагонали или ниже р -ой наддиагонали не превышает р .

Следующий важный шаг в изучении связи малоранговой структуры матрицы и обратной к ней был сделан У. Густафсоном [14], который доказал теорему о размерностях ядер подматриц (теорема 3). Первоначальное доказательство было дано на языке модулей над кольцами. На матричном языке теорема была сформулирована и доказана М. Фидлером и Т. Маркхамом, [15]. Позже они использовали эту теорему, чтобы связать ранги подматриц в прямой и обратной матрице [16]. Этот результат, который мы приводим далее, является прямым обобщением теоремы 6.

Теорема 7 (Фидлер-Маркхам, 1987). Пусть матрица А обратъьма. Тогда г) все подматрицы В выше р -ой наддиагонали матрицы А имеют ранг не выше к: гапк(В) < к тогда и только тогда, когда и) все подматрицы С выше р-ой поддиагонали матрицы А-1 имеют ранг не выше р + к: гапк(С) < р + к.

При р = 0 из этой теоремы следует, что матрица, обратная к квази-сепарабельной, является квазисепарабельной того же порядка. И. Гох-берг, Т. Кайлат и И. Колтрахт [17] первыми вывели алгоритм обращения семисепарабельных матриц, а позже И. Гохберг и Ю. Эйдельман [8] расширили его на квазисепарабельный случай.

Многие матрицы ковариаций, встречающиеся в прикладной статистике, имеют квазисепарабельную структуру. Это объясняет интерес к этим матрицам со стороны ученых-статистиков [18, 19, 20, 21]. Рассмотрим два примера квазисепарабельных матриц из статистики.

Пусть Х1,., Хп — независимые, одинаково распределенные случайные величины, удовлетворяющие Р(Хг = ]') ='Р}> Р1 +р2 + - • .+Рк 1 . Интуитивное событие Хг = означает, что испытание с номером г привело к исходу ). Обозначим через У} количество испытаний из тг, приведших к исходу ). Тогда величины У-\,., У) удовлетворяют мультиномиальному распределению. Матрица ковариаций для У^ имеет вид

С = п

Р1П-Р1) -Р1Р2 -Р1Рз -Р2Р1 Р2(1-Р2) -Р2РЗ -РЗР1 -РЗР2 Рз(1-Рз)

РкР1

РкР2

РкРЗ

-РФк -р2рк -РзРк

Рк(1 -Рк).

Еще один пример возникает при рассмотрении простейшего авторегрессионного процесса АК(1): х1 = ср0 + + сц, сц — белый шум.

Если этот процесс стационарный в узком смысле, то матрица автоко-вариаций для равна

0 = ^ 1 Ф1 ф? ф?

Ф1 1 Фт ф? ф? Ф1 1 ф 1

Ф? ф? Ф1 1

Еще одна область, в которой широко используются квазисепара-бельные матрицы, — это теория ортогональных многочленов. Пусть {гк(х)}£=0 — система многочленов возрастающих степеней: с^т^х) = к. Предположим, что многочлены ортогональны на кривой у, то есть тч(х)г,(х)с!ц(х) =0, у

Тип кривой у влияет на вид рекуррентных соотношений, которым удовлетворяют многочлены. Например, многочлены, ортогональные на отрезке вещественной оси у = [а, Ь] (такие, как многочлены Чебы-шёва, Эрмита, Лежандра и пр.), удовлетворяют трехчленным рекуррентным соотношениям rk+i (х) = (х - ак) • Tic(x) - ßk • rk! (х), k = 0,1,2,., ri(x)=0, т0(х) = 1,

7) а многочлены, ортогональные на единичной окружности (также называемые многочленами Сегё), удовлетворяют двучленным рекуррентным соотношениям

Фк-1 М

X • ф1с1 (х) <Pi(x) " 1 -pk CPlc(x) -Pk 1 cp*k(z) = Zn(pk( 1/z). (8)

Используя вышеприведенные рекуррентные соотношения, в теории ортогональных многочленов доказывается, что матрицы осо (3] 1 ОС)

1 ••. •••

• 0Cn2 ßn-1 1 On-i u =

PoPi 1 О О О

Р0Р2 PlP2 1 О о

РоРз Р1РЗ -Р2РЗ 1 о

-р0р4 -Р1Р4 -Р2Р4 -РЗР4 1

-Р0Р5 — Р1Р5 -Р2Р5 -РЗР5 -Р4Р5

9) биективно связаны с многочленами в (7) и (8) соотношениями тк{х) = сМ(х1к - Тк,к), фк(х) = с1е1(х1к -ик)к), где Тк)к и Ик)к — главные к х к-подматрицы матриц Т и И соответственно.

Заметим, что обе матрицы Т и И из (9) — хессенберговы и ква-зисепарабельные порядка 1 сверху. Следовательно, многочлены, соответствующие квазисепарабельным матрицам, обобщают многочлены, ортогональные на отрезке вещественной оси и единичной окружности. Используя это свойство, можно обобщить многие алгоритмы, основанные на свойствах ортогональности многочленов, на случаи более общих кривых. Такие алгоритмы включают: численное интегрирование (квадратурные формулы Гаусса, Кленшо-Куртиса), интерполирование, Ьи-разложение матриц моментов (алгоритм Шура для матриц Тёплица и алгоритм Ланцоша для матриц Ганкеля). В частности, квазисепарабельная теория может быть использована для более эффективного и высокоточного способа нахождения корней многочленов (как собтвенных значений матриц в (9)). Приведем несколько примеров статей, посвященных связи квазисепарабельных матриц с ортогональными многочленами, и алгоритмам нахождения собственных значений для квазисепарабельных матриц: [22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36].

Невозможно упомянуть о всех работах по теории квазисепарабельных матриц и их приложениях на нескольких страницах. Эта область исследований насчитывает сотни публикаций и быстро развивается. Читатель, желающий углубить свои знания в теории квазисепарабельных матриц, найдет полезными две недавно вышедшие книги [37] и [38]. Хороший обзор литературы о квазисепарабельных матрицах можно также найти в [39]. ьЗ. Основные результаты данной работы

В данной работе решены четыре классические задачи, относящиеся к разным разделам линейной алгебры и вычислительной математики в целом. Объединяющая тема работы — это использование теории квазисепарабельных матриц при решении поставленных задач.

Задача №1: вычисление корней многочленов с высокой точ

Пусть дан многочлен Р(х) степени п, разложенный по базису некоторых других многочленов {ric(x)}£=0 возрастающих степеней:

Р(х) = Гп(х) + rrin-i • Tni (х) + . -I- ттц • Гт (х) + то • г0(х). (10)

Тогда корни многочлена Р(х) совпадают с собственными значениями обобщенной матрицы Фробениуса. Например, в случае многочленов, ортогональных на вещественной оси (7), обобщенная матрица Фробениуса имеет вид

Стандартным способом решения задачи на собственные значения является применение С^И-алгоритма — одного из основных численно устойчивых алгоритмов линейной алгебры. Устойчивость (ЗЫ-алго-ритма означает, что вычисленные в арифметике с плавающей точкой ностью.

00 ßl

1 0И ß2

А= 1 ••• то —mi

И)

• • Oin-2 ßn-1 — тап2

1 ocni — mni собственные значения Л^ являются точными для малого возмущения первоначальной матрицы, то есть

ЛЪЛ2>.} € Л(А + ДА). (12)

Причем возмущение мало по норме: ||AA||f ^ C(u)||A||f для константы С (и) порядка машинного нуля и. Однако из (12) отнюдь не следует, что все Ai являются малыми относительными возмущениями точных собственных значений At первоначальной матрицы А, то есть, что имеет место

Ах - Atl ^ C(u) max{Ai, Ai}, i = 1,2,. (13)

Для плохо обусловленных матриц QR-алгоритм не может гарантировать высокую относительную точность малых собственных значений. Следовательно, малые корни многочлена Р(х) в (10), вычисленные QR-алгоритмом, примененным к обобщенной матрице Фробениу-са, могут иметь большую относительную ошибку. Причем эта ошибка никак не связана с высокой чувствительностью корней к возмущениям коэффициентов mic в разложении (10). Это свойство алгоритма, а не задачи.

Естественно попытаться построить алгоритм, который для матриц типа (11) находит приближения собственных значений с высокой относительной точностью (13). Такой алгоритм известен по крайней мере для одного класса матриц — симметричных трехдиагональ-ных положительно определенных. Алгоритм был получен К. Фернандо и Б. Парлеттом в 1994 году [40] и получил название дифференциальный qd-алгоритм (или dqds, для краткости). Относительная точность (13) алгоритма dqds была доказана только в симметричном трехдиагональном положительно определенном случае. Однако на практике она наблюдается даже в общем несимметричном случае (хотя не гарантирована).

В нашей работе мы обобщили dqds-алгоритм Фернандо и Парлет-та с трехдиагонального на хессенберговый квазисепарабельный случай. Хессенберговы квазисепарабельные матрицы связаны со многими важными системами многочленов: мономами; многочленами, ортогональными на отрезке вещественной оси (Чебышёва, Эрмита, Лагран-жа и др.); многочленами, ортогональными на единичной окружности (Сегё). Следовательно, наш алгоритм может быть применен для вычисления корней многочленов, представленных в базисах вышеописанных семейств многочленов. Хотя нам не удалось доказать свойство

13) для нового алгоритма, в проделанных нами численных экспериментах предложенный алгоритм продемонстрировал более высокую точность, чем стандартный алгоритм С^Ы.

Задача №2: обращение полиномиальных матриц Вандермонда.

Матрица Вандермонда

14) ух = " 1 1 XI • х2 • • х2

1 Хп ■ ' лтг . является, пожалуй, одной из наиболее известных структурированных матриц. Ни один учебник по линейной алгебре не обходится без примера или задачи (чему равен определитель?), посвященной матрице Вандермонда.

Известно, что для узлов {х]с}£=1, равномерно распределенных на отрезке, матрица Вандермонда плохо обусловлена. В 1994 году Е. Тыр-тышниковым было доказано, что ее число обусловленности растет экспоненциально с размером [41]. Поэтому инженеры пытаются избегать использования матрицы Вандермонда в вычислениях. Однако в 1966 Д. Траубом [42] был получен алгоритм обращения матрицы Вандермонда за 0{п2) арифметических операций с малой ошибкой в смысле прямого анализа ошибок округления. Казалось бы, это невозможно: «прямая» ошибка равняется обратной (не меньше машинного нуля) умноженной на число обусловленности (огромно!), однако весь «трюк» заключается в том, что алгоритм Трауба структурный (О (п2) против 0{п3)), и к нему неприменим стандартный результат теории ошибок округления.

Матрица Вандермонда (14) связана с задачей интерполяции в базисе мономов {хк}£~о . Заменим мономы произвольной системой многочленов возрастающих степеней (т^х)}]^ , с^г^х) = к. Тогда получим матрицу более общего вида, называемую полиномиальной матрицей Вандермонда:

Го[х2) Г!(х2)

Гп-1 г П-1 (х2)

1*0 (*п) Т1(Хп) ••• Гп|(Хп)

15)

Частными случаями этого класса матриц являются такие хорошо известные всем матрицы, как матрица дискретного преобразования Фурье и матрица дискретного преобразования косинусов.

Следуя успеху Трауба, многие исследователи в области численной линейной алгебры обратились к задаче обращения матрицы (15) для различных базисов многочленов. Л. Райхел и Г. Опфер [43] обобщили алгоритм Трауба на случай многочленов Чебышёва, а В. Ольшевский — на случай многочленов Сегё [44]. Прочие работы, посвященные этой проблеме, включают [45, 46, 47, 48].

В настоящей работе мы построили алгоритм обращения полиномиальной матрицы Вандермонда (15) для абсолютно произвольной системы многочленов {^(х)}^ , с^г^х) = к. Причем наш алгоритм быстрый (0{п2) арифметических операций) всегда, когда многочлены заданы короткими рекуррентными соотношениями, и тривиально сводится к алгоритмам, полученным ранее для частных случаев систем многочленов. То есть наш результат полностью обобщает все предыдущие результаты на эту тему.

Задача №3: решение задач с седловой точкой. Оптимизационные задачи с ограничениями в виде уравнений в частных производных (УЧП) возникают повсеместно в процессе инженерного проектирования. Рассмотрим простейшую задачу такого типа: шт1||и-й||2 + (3||Я|2,

-У2и = !вП, (16) и = д на дС1.

Дискретизуем эти уравнения методом конечных элементов и поставим вопрос об отыскании стационарной точки лагранжиана полученной оптимизационной задачи. Тогда (16) сведется к системе линейных уравнений с блочной матрицей особого вида:

2(3 М. 0 -М" "0"

0 М. кт и = ъ

-М. К 0 Л £1 где М. — матрица масс, а К — матрица жесткости. Это — так называемая задача с седловой точкой. Ее решение представляет собой достаточно большую вычислительную сложность, так как матрица системы неопределенная, плохо обусловленная и часто большого размера. Огромное количество работ посвящено развитию методов для задач с седловой точкой. Хороший обзор результатов в этой области можно найти в [49].

Систему уравнений (17) чаще всего решают с помощью итерационных методов, а техники предобусловливания основаны на приближенном решении системы с матрицей блочного дополнения по Шуру для блока в позиции (3,3):

Б = — ^М + КМ1КТ^ . (18)

Эта матрица — отрицательно определенная, и для решения системы уравнений с ней может быть применен, например, метод сопряженных градиентов. Главная сложность заключается в том, что матрица Б плотная, и ее невозможно вычислить и/или сохранить в памяти. Один из подходов — это использовать только одно из слагаемых — ^ М. или —КМ.-1 Кт в качестве приближения матрицы Б. Очевидно, такое приближение в общем случае не очень хорошее, и число итераций в методе сопряженных градиентов может быть велико.

Нами предложен абсолютно новый подход к решению задач с сед-ловой точкой вида (17). Используя особую малоранговую структуру дополнения по Шуру Б из (18), нам удалось построить методы, которые позволяют найти приближение с высокой точностью плотной матрицы Б, ее Ъи-факторизацию и решить систему с ней за асимптотически 0[п) арифметических операций. Использование памяти также растет линейно с размером задачи. Наш подход в настоящее время развит только для тензорных сеток дискретизации, расширение его области применения — цель будущих исследований.

Задача №4: исследование численной устойчивости алгоритмов с квазисепарабельными матрицами.

Быстрые алгоритмы решения систем с квазисепарабельными матрицами были значительно развиты в последнее десятилетие. Множество новых алгоритмов было предложено несколькими исследовательскими группами из разных стран. Однако, несмотря на интенсивные исследования в этой области, ни для одного из этих алгоритмов не был проведен строгий анализ ошибок округления. Возможно, это объясняется тем, что алгоритмы с квазисепарабельными матрицами используют малопараметрические представления последних, и стандартные методы анализа ошибок округления (см. Н. Хайям «Точность и устойчивость численных алгоритмов» [50]) неприменимы к таким алгоритмам.

В настоящей работе нами получены первые результаты анализа численной устойчивости алгоритмов с квазисепарабельными матрицами. Мы рассмотрели два алгоритма решения квазисепарабельных систем уравнений с помощью С^И-разложения. Для краткости будем называть их алгоритмами А и Б. Алгоритм А в его наиболее общем виде был предложен в [51, алгоритмы 6.4 и 7.2], в то время, как алгоритм Б был выведен в [52, §5.3], а также описал в [37, глава 5]. Мы доказали численную устойчивость алгоритма А, а именно следующую теорему:

Теорема 8. Пусть А — п х тг квазисепарабельная матрица с диагональю, включенной в малоранговую структуру с одной стороны, и Ь е Мп. Предположим, система уравнений Ах = Ь решена в арифметике с плавающей точкой с помощью алгоритма А. Тогда вычисленное решение х является точным решением системы

А + ЛА)х = Ъ + ЛЬ, где

АА(:)))||2^ук1П2||А(:,))||2> ) = 1,.,п, ||ЛЬ||2 ^ Ук2п||Ь||2 и К1, К2 — малые (натуральные) константы.

Что касается алгоритма Б, то нами показано, что этот алгоритм численно неустойчив (несмотря на использование им (311-фактори-зации!). Мы объяснили причины его неустойчивости и подтвердили наши выводы результатами численных экспериментов.

1.4. Содержание работы по главам

Первая глава посвящена вычислению корней многочлена, разложенного по базису произвольных многочленов. В §1.2 мы объясним связь между этой задачей и задачей на собственные значения для хес-сенберговых квазисепарабельных матриц. Там же мы получим алгоритм Ьи-разложения квазисепарабельной матрицы линейной по размеру сложности. В §§1.3-1.4 мы докажем инвариантность квазисепарабельной структуры матрицы относительно итераций ЬЫ алгоритма и, основываясь на этом результате, обобщим с^с^-алгоритм Фернандо и Парлетта на случай хессенберговых квазисепарабельных матриц. Далее, в §1.5, мы дадим альтернативный вывод нового алгоритма через обобщенный процесс Грама-Шмидта. Этот вывод менее технический и проясняет значение некоторых промежуточных величин, возникающих при работе алгоритма. Главной задачей нашего исследования было отыскание корней многочленов с высокой относительной точностью. Поэтому в §1.6 мы уточняем новый алгоритм на случай обобщенных матриц Фробениуса для систем мономов и многочленов, ортогональных на отрезке вещественной оси. Глава заканчивается §1.7, в котором мы приводим результаты численных экспериментов с многочленами, такими, как многочлен Уилкинсона, заданными в базисах мономов и многочленов Чебышёва.

Во второй главе решается задача об обращении полиномиальной матрицы Вандермонда для произвольной системы многочленов. Мы выводим новый алгоритм, который обобщает многие алгоритмы для частных случаев, полученные другими исследователями ранее. Наш вывод основал на использовании междисциплинарной связи между тремя разделами науки: (1) электротехникой, (11) информатикой и (ш) линейной алгеброй. Из информатики мы берем так называемый принцип Теллегена, который мы рассматриваем в §2.2. Из электротехники нам потребуется понятие о графе потока сигнала (§2.3). В параграфе §2.3 мы свяжем два вышеперечисленных понятия с обращением полиномиальной матрицы Вандермонда и, использую эту связь, окончательно выведем новый алгоритм в §2.5. В заключительном §2.6 мы применяем новый алгоритм для обращения матрицы Вандермонда для системы многочленов, заданных I-членными рекуррентными соотношениями. Такой алгоритм, насколько нам известно, не был опубликован ранее. Результаты численных экспериментов демонстрируют высокую точность нового алгоритма.

В третьей главе мы применяем теорию квазисепарабельных матриц к решению задач с седловой точкой, возникающих в оптимизационных задачах с ограничениями в виде уравнений в частных производных. В §3.2 описывается одна такая модельная задача и объясняются трудности, возникающие при ее решении. В §3.3 мы вводим новый класс структурированных матриц — многоуровневые квазисепара-бельные матрицы. Новые матрицы естественным образом обобщают обыкновенные квазисепарабельные матрицы на случай задач многих переменных. Мы показываем, что многоуровневые квазисепарабельные матрицы могут быть применены при решении уравнений в частных производных и задач оптимального управления с ними. Результаты численных экспериментов такого рода приведены в §3.4. Они демонстрируют, что когда многоуровневые квазисепарабельные матрицы используются для построения предобусловливателя в итерационном методе решения системы с седловой точкой, общая сложность алгоритма и память растут линейно с ростом размера задачи.

Заключительная, четвертая глава посвящена изучению численной устойчивости двух: алгоритмов [51, алгоритмы 6.4 и 7.2] и [52, §5.3] решения систем с квазисепарабельными матрицами. Оба этих алгоритма основаны на быстром С^К-разложении и обладают линейной по размеру задачи сложностью. В §4.2 мы описываем общие идеи этих алгоритмов. Первый из алгоритмов детально рассматривается в §4.3, где мы доказываем, что этот алгоритм численно устойчив, а именно, что в арифметике с плавающей точкой он находит точное решение системы для малого (по норме) возмущения матрицы системы и вектора в правой части. Напротив, как мы показываем в §4.4, второй из расматриваемых нами алгоритмов неустойчив. В том же §4.4 мы объясняем причины его неустойчивости и подтверждаем наши выводы численными экспериментами. Главный вывод, который можно сделать из результатов этой главы, это то, что в случае квазисепа-рабельных матриц устойчивость алгоритма не есть прямое следствие его устойчивости в неструктурированном случае и требует отдельного исследования.

Заключение диссертации по теме «Вычислительная математика», Жлобич, Павел Георгиевич

4.5. Выводы

В этой главе мы получили первые результаты анализа ошибок округления быстрых алгоритмов для квазисепарабельных матриц. Мы проанализировали два алгоритма решения систем через структурное С^Ы-разложение для подкласса квазисепарабельных матриц порядка один, заданных уравнением (4.1), то есть матрицами, для

Рис. 4.2. Относительные невязки в 1000 экспериментах со случайными матрицами размера 100 х 100.

18 — 1« 14 12 10 8 6 4 2 0

Заключение

В заключение диссертации сформулируем ее основные результаты.

1. Предложен новый алгоритм нахождения корней многочленов, представленных в произвольных базисах. В случае, когда многочлены из базиса удовлетворяют коротким рекуррентным соотношениям, таким как трехчленные для ортогональных многочленов, сложность алгоритма составляет всего 0[п) арифметических операций на шаг. Результаты численных экспериментов подтверждают, что итерации нового алгоритма сходятся быстрее и приводят к более точным приближениям корней многочленов, чем стандартные методы.

2. Установлена связь между обращением полиномиальных матриц Вандермонда (линейная алгебра), обращением графа потока сигнала (электротехника) и принципом Теллегена (теоретическая информатика). Получен алгоритм обращения полиномиальной матрицы Вандермонда для произвольного базиса многочленов. Тем самым обобщены многие алгоритмы, предложенные ранее для частных случаев.

3. Определен и изучен новый класс малорангово-структурирован-ных матриц — многоуровневые квазисепарабельные матрицы. Новый класс матриц применен при решении задачи оптимального управления с ограничениями в виде уравнений в частных производных. Результаты численных экспериментов продемонстрировали линейный рост числа операции в предложенном методе с ростом размера задачи.

4. Изучена численная устойчивость двух алгоритмов решения системы уравнений с квазисепарабельной матрицей. Для одного из алгоритмов доказана теорема о численной устойчивости, в то время, как для другого построен контрпример, демонстрирующий численную неустойчивость, и объяснены ее причины.

Результаты, полученные в данной диссертации, демонстрируют, что использование квазисепарабельных матриц во многих областях вычислительной математики позволяет по новому взглянуть на существующие задачи и построить новые быстрые алгоритмы их решения.

Список литературы диссертационного исследования кандидат физико-математических наук Жлобич, Павел Георгиевич, 2012 год

1. Tyrtyshnikov Е. Mosaic-skeleton approximations // Calcolo. 1996. V. 33, №1-2. P. 47-57.

2. Hackbusch W., Nowak Z. On the fast multiplication in the boundary-elements method by panel clustering // Numerische Mathematik. 1989. V. 54, № 4. P. 463-491.

3. Rokhlin V. Rapid solution of integral equations of classical potential theory // Journal of Computational Physics. 1985. V. 60. P. 187-207.

4. Tyrtyshnikov E. Incomplete cross approximation in the mosaic-skeleton method // Computing. 2000. V. 4. P. 367-380.

5. Горейнов С.А., Замарашкин Н.Л., Тыртышников E.E. Псевдоске-летонные аппроксимации матриц // Доклады российской академии наук. 1995. Т. 343, № 2. С. 151-152.

6. Goreinov S., Tyrtyshnikov Е., Yeremin A. Matrix-free iteration solution strategies for large dense linear systems // Numerical Linear Algebra with Applications. 1996. V. 4, № 5. P. 1-22.

7. Hackbusch W. A sparse matrix arithmetic based on H -matrices. Part I: Introduction to H-matrices // Computing. 1999. V. 62. P. 89-108.

8. Eidelman Y., Gohberg I. On a new class of structured matrices I I Integral Equations and Operator Theory. 1999. V. 34. P. 293-324.

9. Гантмахер Ф.Р., Крейн, М.Г. Осцилляционные матрицы и малые колебания механических систем // М.-Л.: ГТТИ. 1941. 220 с.

10. Barrett W. A theorem on inverse of tridiagonal matrices // Linear Algebra and Its Applications. 1979. V. 27. P. 211-217.

11. Ikebe Y. On inverses of Hessenberg matrices // Linear Algebra and its Applications. 1979. V. 24. P. 93-97.

12. Cao W., Stewart W. A note on inverses of Hessenberg-like matrices // Linear algebra and its applications. 1986. V. 76. P. 233-240.

13. Asplund E. Inverses of matrices {ciij} which satisfy ciij = 0 for j > i + p j I Math. Scand. 1959. V. 7. P. 57-60.

14. Gustafson W. A note on matrix inversion // Linear Algebra and Its Applications. 1984. V. 57. P. 71-73.

15. Fiedler M., Markham T. L. Completing a matrix when certain entries of its inverse are specified // Linear Algebra and its Applications. 1986. V. 74. P. 225-237.

16. Fiedler M., Markham T. Rank-preserving diagonal completions of a matrix // Linear Algebra and its Applications. 1987. V. 85. P. 49-56.

17. Gohberg I, Kailath T., Koltracht I. Linear complexity algorithms for semiseparable matrices // Integral Equations and Operator Theory. 1985. V. 8, № 6. P. 780-804.

18. Graybill F. Matrices with applications in statistics, Ed. by 2nd. — Belmont, CA: The Wadsworth Statistics/Probability Series, 1983.

19. Barrett W., Feinsilver P. Gaussian families and a theorem on patterned matrices // Journal of Applied Probability. 1978. P. 514522.

20. Greenberg G., Sarhan A. Matrix inversion, its interest and application in analysis of data // Journal of the American Statistical Association. 1959. P. 755-766.

21. Roy S., Greenberg B., Sarhan A. Evaluation of determinants, characteristic equations and their roots for a class of patterned matrices // Journal of the Royal Statistical Society Series B (Methodological). 1960. P. 348-359.

22. Bini D., Daddi F., Gemignani L. On the shifted QR iteration applied to companion matrices // Electronic Transactions on Numerical Analysis. 2004. V. 18. P. 137-152.

23. Vandebril R. Semiseparable matrices and the symmetric eigenvalue problem: Ph.D. thesis / KU Leuven. — 2004.

24. Chandrasekaran S., Gu M. A divide-and-conquer algorithm for the eigendecomposition of symmetric block-diagonal plus semiseparable matrices // Numerische Mathematik. 2004. V. 96, № 4. P. 723-731.

25. Bini D., Gemignani L., Pan V. Fast and stable QR eigenvalue algorithms for generalized companion matrices and secular equations // Numerische Mathematik. 2005. V. 100, № 3. P. 373-408.

26. Eidelman Y., Gohberg I, Olshevsky V. Eigenstructure of order-one-quasiseparable matrices. Three-term and two-term recurrence relations // Linear Algebra and its Applications. 2005. V. 405. P. 140.

27. Eidelman Y., Gohberg I, Olshevsky V. The QR iteration method for Hermitian quasiseparable matrices of an arbitrary order // Linear Algebra and its Applications. 2005. V. 404. P. 305-324.

28. Mastronardi N., Van Camp E., Van Barel M. Divide and conquer algorithms for computing the eigendecomposition of symmetric diagonal-plus-semiseparable matrices // Numerical Algorithms. 2005. V. 39, № 4. P. 379-398.

29. Van Barel M., Van Camp E., Mastronardi N. Orthogonal similarity transformation into block-semiseparable matrices of semiseparability rank k // Numerical linear algebra with applications. 2005. V. 12, № 10. P. 981-1000.

30. Chandrasekaran S., Gu M., Xia J., Zhu J. A fast QR algorithm for companion matrices // Recent Advances in Matrix and Operator Theory. 2008. V. 179. P. 111-143.

31. Plestenjak B., Van Barel M., Van Camp E. A Cholesky LR algorithm for the positive definite symmetric diagonal-plus-semiseparable eigenproblem // Linear Algebra and its Applications. 2008. V. 428, №2-3. P. 586-599.

32. Bella T., Eidelman Y., Gohberg I, Olshevsky V. Classifications of three-term and two-term recurrence relations via subclasses of quasiseparable matrices // SIAM Journal of Matrix Analysis. 2009.

33. Bini D., Boito P., Eidelman Y. et al. A fast implicit QR eigenvalue algorithm for companion matrices // Linear Algebra and its Applications. 2010. V. 432, № 8. P. 2006-2031.

34. Vandebril R., Van Barel M., Mastronardi N. Matrix Computations and Semiseparable Matrices, Vol. I: Linear Systems. — The Johns Hopkins University Press, 2008.

35. Vandebril R., Van Barel M., Mastronardi N. Matrix Computations and Semiseparable Matrices, Vol. II: Eigenvalue and Singular Value Methods. — Johns Hopkins University Press, 2008.

36. Vandebril R., Van Barel M., Golub G., Mastronardi N. A bibliography on semiseparable matrices // Calcolo. 2005. V. 42, № 3. P. 249-270.

37. Fernando K., Parlett B. Accurate singular values and differential qd algorithms // Numerische Mathematik. 1994. V. 67, № 2. P. 191-229.

38. Tyrtyshnikov E. How bad are Hankel matrices? // Numerische Mathematik. 1994. V. 67, № 2. P. 261-269.

39. Traub J. Associated polynomials and uniform methods for the solution of linear problems // SIAM Review. 1966. V. 8, № 3. P. 277301.

40. Reichel L., Opfer G. Chebyshev-Vandermonde systems // Mathematics of Computation. 1991. V. 57. P. 703-721.

41. Calvetti D., Reichel L. Fast inversion of Vandermonde-like matrices involving orthogonal polynomials // BIT. 1993.

42. Gohberg I., Olshevsky V. Fast inversion of Chebyshev-Vandermonde matrices // Numerische Mathematik. 1994. V. 67, X2 1. P. 71-92.

43. Gohberg I., Olshevsky V. The fast generalized Parker-Traub algorithm for inversion of Vandermonde and related matrices // Journal of Complexity. 1997. V. 13, № 2. P. 208-234.

44. Kailath T., Olshevsky V. Displacement structure approach to polynomial Vandermonde and related matrices // Linear Algebra and Its Applications. 1997. V. 261. P. 49-90.

45. Benzi M., Golub G., Liesen J. Numerical solution of saddle point problems // Acta Numerica. 2005. V. 14. P. 1-137.

46. Higham N. Accuracy and Stability of Numerical Algorithms. — Philadelphia: SIAM, 2002.

47. Eidelman Y., Gohberg I. A modification of the Dewilde-van der Veen method for inversion of finite structured matrices // Linear Algebra and its Applications. 2002. V. 343. P. 419-450.

48. Van Camp E., Mastronardi N., Van Barel M. Two fast algorithms for solving diagonal-plus-semiseparable linear systems // Journal of Computational and Applied Mathematics. 2004. V. 164. P. 731-747.

49. Maroulas J., Barnett S. Polynomials with respect to a general basis // Journal of Mathematical Anedysis and Applications. 1979. V. 72. P. 177-194.

50. Vandebril R., Van Barel M., Mastronardi N. An implicit QR algorithm for symmetric semiseparable matrices // Numerical Linear Algebra with Applications. 2005. V. 12, № 7. P. 625-658.

51. Bevilacqua R., Bozzo E., Del Corso G. qd-type methods for quasiseparable matrices // SIAM Journal on Matrix Analysis and Applications. 2011. V. 32. P. 722-747.

52. Parlett B. The new qd algorithms // Acta Numerica. 1995. V. 4, № 1. P. 459-491.

53. Demmel J., Kahan W. Accurate singular values of bidiagonal matrices // SIAM Journal on Scientific and Statistical Computing. 1990. V. 11, № 5. P. 873-912.

54. Xia J., Chandrasekaran S., Gu M., Li X. Superfast multifrontal method for large structured linear systems of equations // SIAM Journal on Matrix Analysis and Applications. 2009. V. 31, № 3.1. P. 1382-1411.

55. Delvaux S., Van Baiel M. A Givens-weight representation for rank structured matrices // SIAM Journal on Matrix Analysis and Applications. 2007. V. 29, № 4. P. 1147-1170.

56. Eidelman Y., Gohberg I. On generators of quasiseparable finite block matrices // Calcolo. 2005. V. 42. P. 187-214.

57. Clenshaw C. A note on the summation of Chebyshev series // MTAC, v. 9. 1955. V. 17, №194. P. 118-120.

58. Parlett B., Reinsch C. Balancing a matrix for calculation of eigenvalues and eigenvectors // Numerische Mathematik. 1969. V. 13, № 4. P. 293-304.

59. Wilkinson J. The evaluation of the zeros of ill-conditioned polynomials. Part I // Numerische Mathematik. 1959. V. 1, № 1. P. 150-166.

60. Bella T., Eidelman Y., Gohberg I. et aI. A Fast Traub-like inversion algorithm for Hessenberg order one quasiseparable Vandermonde matrices // Journal of Complexity, to appear.

61. Bordewijk J. L. Inter-reciprocity applied to electrical networks // AppJ. Sci. Res. B. 1956. V. 6. P. 1-74.

62. Fiduccia C. M. On obtaining upper bounds on the complexity of matrix multiplication // Complexity of computer computations (Proc. Sympos., IBM Thomas J. Watson Res. Center, Yorktown Heights, N.Y., 1972). Plenum, New York: 1972. — P. 31-40.

63. Fiduccia C. M. On the algebraic complexity of matrix multiplication: Ph.D. thesis / Brown University. — 1973.

64. Biirgisser P., Clausen M., Shokrollahi M. A. Algebraic complexity theory. — Berlin: Springer-Verlag, 1997.

65. Parker F. Inverses of Vandermonde matrices // Amer. Math. Monthly. 1964. V. 71. P. 410-411.

66. Forney G. Concatenated codes. — Cambridge: The M.I.T. Press, 1966.

67. T.Kailath. Linear Systems. — Englewood Cliffs, NJ: Prentice Hall, 1980.1 ЧЧ —

68. A.V.Oppenheim, R.W.Schafer. Discrete-time Signal Processing. — Eglewood Cliffs, NJ: Prentice Hall, 1989.

69. Olshevsky V. Eigenvector computation for almost unitary Hessenberg matrices and inversion of Szego-Vandermonde matrices via discrete transmission lines // Linear Algebra and its Applications. 1998. V. 285. P. 37-67.

70. Геронимус Л. Полиномы, ортогональные на круге, и их приложения // Сообщения Харьковского математического общества. 1948. Т. 19. С. 35-120.

71. Barnett S. Congenial matrices // Linear Algebra and its Applications. 1981. V. 41. P. 277-298.

72. Bella Т., Olshevsky V., Zhlobich P. Signal flow graph approach to inversion of (H,m)-quasiseparable-Vandermonde matrices and new filter structures // Linear Algebra and its Applications. 2010. V. 432, № 8. P. 2032-2051.

73. Hackbusch W., Khoromskij В., Sauter S. On ^-matrices // Lectures on Applied Mathematics / Ed. by H. Bungartz, R. Hoppe, C. Zenger. — Berlin: Springer-Verlag, 2000. — P. 9-29.

74. Eidelman Y., Gohberg I. Linear complexity inversion algorithms for a class of structured matrices // Integral Equations and Operator Theory. 1999. V. 35. P. 28-52.

75. Martinsson P. A fast direct solver for a class of elliptic partial differential equations // Journal of Scientiñc Computing. 2009. V. 38, № 3. P. 316-330.

76. Агошков В.И. Методы оптимального управления и сопряженных уравнений в задачах математической физики. — М.: ИВМ РАН, 2003. 256 с.

77. Rees Т., Dollar Н., Wathen A. Optimal solvers for PDE-constrained optimization // SIAM Journal on Scientiñc Computing. 2010. V. 32, № 1. P. 271-298. doi: 10.1137/080727154. http://link.aip.org/link/SJ0CE3/v32/il/p271/sl&Agg=doi.

78. Elman H., Silvester D., Wathen A. Finite elements and fast iterative solvers: with applications in incompressible fluid dynamics. — Oxford University Press, USA, 2005.

79. Paige C., Saunders M. Solution of sparse indefinite systems of linear equations // SIAM Journal on Numerical Analysis. 1975. P. 617-629.

80. Gould N., Hribar M., Nocedal J. On the solution of equality constrained quadratic programming problems arising in optimization // SIAM Journal on Scientific Computing. 2002. V. 23, № 4. P. 1376-1395.

81. Chandrasekaran S., Dewilde P., Gu M., Somasunderam N. On the numerical rank of the off-diagonal blocks of Schur complements of discretized elliptic PDEs // SIAM Journal on Matrix Analysis and Applications. 2010. V. 31. P. 2261-2290.

82. Folland G. Introduction to partial differential equations. — Princeton, NJ: Princeton University Press, 1995.

83. Lions J., Mitter S. Optimal control of systems governed by partial differential equations. — Springer Berlin, 1971. — V. 396.

84. Chandrasekaran S., Gu M. Fast and stable algorithms for banded plus semiseparable systems of linear equations // SIAM Journal on Matrix Analysis and Applications. 2004. V. 25, № 2. P. 373-384.

85. Delvaux S.} Van Barel M. A QR-based solver for rank structured matrices // SIAM Journal on Matrix Analysis and Applications. 2008. V. 30, № 2. P. 464-490.

86. Olshevsky V., Strang G., Zhlobich P. Green's matrices // Linear Algebra and its Applications. 2010. V. 432. P. 218-241.

87. Vandebril R., Van Barel M., Mastronardi N. A note on the representation and definition of semiseparable matrices // Numerical Linear Algebra with Applications. 2005. V. 12, № 8. P. 839—858.

88. Gondzio J., Zhlobich P. Multilevel quasiseparable matrices in PDE-constrained optimization. — submitted to Numerische Mathematik.

89. Dopico F., Olshevsky V., Zhlobich P. Stability of QR-based fast system solvers for a subclass of quasiseparable rank one matrices. — submitted to Mathematics of Computation.

90. Zhlobich P. Differential qd algorithm with shifts for rank-structured matrices. — SIAM Journal on Matrix Analysis and Applications (SIM AX).

91. Olshevsky V., Tyrtyshnikov E., Zhlobich P. Tellegen's principle, non-minimal realization of systems and inversion of polynomial Vandermonde matrices I I Russian Journal of Numerical Analysis and Mathematical Modelling. 2012. V. 27, № 2. P. 131-154.

92. Bella T., Olshevsky V., Zhlobich P. A quasiseparable approach to five-diagonal CMV and Fiedler matrices // Linear Algebra and its Applications. 2011. V. 434, № 4. P. 957-976.