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

На правах рукописи

00505567/

Жлобич Павел Георгиевич

КВАЗИСЕПАРАБЕЛЬНЫЕ МАТРИЦЫ В ЛИНЕЙНОЙ АЛГЕБРЕ И ЕЕ ПРИЛОЖЕНИЯХ

01.01.07 — Вычислительная математика

АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата физико-математических наук

Москва — 2012

005055672

Работа выполнена в Федеральном государственном автономном образова тельном учреждении высшего профессионального образования "Москов ский физико-технический институт (государственный университет)"

Научный руководитель:

член-корреспондент РАН, профессор Тыртышников Евгений Евгеньевич Официальные оппоненты:

доктор физико-математических наук, профессор кафедры общей математики факультета ВМК МГУ Икрамов Хаким Дододжанович

кандидат физико-математических наук, руководитель отдела математического моделирования ОАО "Центральная Геофизическая Экспедиция" Книжнерман Леонид Аронович

Ведущая организация:

Федеральное государственное бюджетное учреждение науки Вычислительный центр им. А. А. Дородницына Россиийской академии наук

Защита состоится « 5 » октября 20 12 г. в 15:00 на заседании диссертационного совета Д 002.045.01 в Федеральном государственном бюджетном учреждении науки Институте вычислительной математики Российской академии наук, расположенном по адресу: 119333, г. Москва, ул. Губкина, д. 8.

С диссертацией можно ознакомиться в библиотеке Федерального государственного бюджетного учреждения науки Института вычислительной математики Российской академии наук.

Автореферат разослан сентября 20 12 г.

Ученый секретарь

диссертационного совета Д 002.045.01 доктор физико-математических наук

Г. А. Бочаров

Общая характеристика работы

Актуальность работы

Предлагаемая диссертация посвящена развитию теории матриц с малоранговой структурой особого вида, называемых квазисепарабелъными. Квазисепарабельная матрица, в некотором смысле, является дискретным аналогом функции Грина для одномерного оператора Штурма-Лиувилля. Так как функция Грина — это ядро соответствующего интегрального оператора, то очевидно, что квазисепарабельные матрицы являются подклассом более общих мозаично-скелетонных матриц Е. Е. Тыртышникова и иерархических матриц В. Хакбуша. Однако, в отличие от последних, квазисепарабельные матрицы обладают простой нерекурсивной линейной по размеру параметризацией, то есть, полное число ячеек памяти, необходимых для хранения квазисепарабельной матрицы, равняется 0[п). Численные методы с квазисепарабельными матрицами активно развивались в последнее десятилетие и, можно сказать, что к моменту написания диссертации были получены все основные алгоритмы линейной алгебры для этого класса матриц.

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

Цели работы

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

Перед автором были поставлены четыре задачи:

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

2. обобщить все ранние работы по обращению полиномиальных матриц Вандермонда на случай произвольной системы многочленов;

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

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

Методы исследования

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

а,

Р2Ч1 Рз^гЧ!

9^2 ¿2 РзЧ2

д1Ь2Н3 ЭгНз

а3

Э1Ь2. . .Ьп^Ьп дгЪз.. .Ьп^Ьп дзЪ4 — Ьп._1Нтг

.РпОп-1 • • • 1241 Рп<1п-1 • • • 4342 РпО-п—1... сцЯз

Здесь параметры {с1к, Цк, а^, рк, Эк, Ь^, Нк}, называемые генераторами, — это матрицы малых размеров. Как видно из этой параметризации, ква-зисепарабельная матрица является дискретным аналогом функции двух переменных, которые разделяются в разных подобластях области определения. Полное количество ячеек памяти, необходимых для хранения такой матрицы, растет линейно с ее размером.

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

Научная новизна работы

Предлагаемый в работе алгоритм вычисления собственных значений хессенберговых квазисепарабельных матриц (алгоритм 4), является новым и принадлежит автору. Автору не известны другие алгоритмы решения той же задачи, обладающие сравнимой или более низкой алгебраической сложностью. Теорема 26, определяющая структуру матриц, обратных к полиномиальным матрицам Вандермонда, и задающая алгоритм их быстрого обращения, получена в соавторстве с В. Р. Ольшевским и Е. Е. Тыртышниковым. Идея использования графов потока сигнала и принципа Теллегена при ее доказательстве принадлежит В. Р. Ольшевскому. Класс многоуровневых квазисепарабельных матриц (определение 31) изобретен автором и не был известен ранее. Метод решения систем уравнений с седловой точкой, использующий многоуровневую квазисепа-рабельную структуру, является новым и принадлежит автору. Теорема 41, устанавливающая численную устойчивость алгоритма решения системы с квазисепарабельной матрицей, доказана автором совместно с Ф. М. Допи-ко и В. Р. Ольшевским.

Защищаемые положения

На защиту выносятся основные результаты работы:

1. Алгоритм LR-типа вычисления собственных значений хессенберго-вых матриц с малоранговой структурой в верхней части и вытекающий из него новый метод вычисления корней многочлена, разложенного по произвольному ортогональному базису, сложности 0[п) на корень.

2. Полное описание матриц, обратных к полиномиальным матрицам Вандермонда, а также алгоритм их обращения сложности 0{п2) в структурированном случае.

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

4. Теорема о численной устойчивости быстрого алгоритма решения системы уравнений с квазисепарабельной матрицей, предложенного в Eidelman Y. and 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.

Практическая значимость работы

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

Поиск корней многочленов является одной из важнейших задач алгебры, часто возникающих в инженерных вычислениях. С точки зрения линейной алгебры эта задача может быть переформулирована, как задача на собственные значения для обобщенной матрицы Фробениуса. Один из основных результатов диссертации — это алгоритм решения задачи на собственные значений таких матриц, использующий их структуру. Предложенный алгоритм обобщает алгоритм Фернандо и Парлетта для трехдиагональных матриц, который в симметричном положительно-определенном случае гарантирует высокую точность вычисленных собственных значений. У нас нет доказательства подобного утверждения для нового алгоритма, однако на практике его точность намного превосходит

точность <311-алгоритма. К тому же, сложность предложенного алгоритма составляет всего 0{п) арифметических операций на корень, что на порядок меньше, чем 0{п2) для (^Я-алгоритма.

Другие важные с практической точки зрения результаты диссертации относятся к области оптимального управления в процессах, чья динамика описывается уравнениями в частных производных. Такие оптимизационные задачи возникают повсеместно в процессе инженерного проектирования, например, при построении летательных аппаратов и химических реакторов. Дискретизация задач управления с ограничениями в виде уравнений в частных производных приводит к системам с седловой точкой больших размеров. Численное решение таких систем — непростая задача из-за неопределенности и плохой обусловленности матриц систем. Нами предложен новый подход к решению систем с седловой точкой, использующий малоранговую структуру матриц и обратных к ним. Наш метод достаточно универсальный и может быть применен к задачам с разными дифференциальными операторами. Численные эксперименты с простейшим уравнением теплопроводности показали линейную сложность предложенного алгоритма относительно числа точек дискретизации. На данный момент область применения алгоритма ограничена двумерными задачами, дискретизованными на равномерных тензорных сетках. Однако мы надеемся, что идеи многоуровневой малоранговой аппроксимации матриц, предложенные в данной работе, в будущем будут развиты и применены к решению задач более общего вида.

Для практической реализации любого алгоритма на ЭВМ важно знать, является ли последний численно устойчивыми. Численная устойчивость означает, что, если все арифметические операции в процессе счета алгоритма произведены с малой ошибкой е, то результат вычислений является точным для малого О(г) возмущения первоначальных данных. Очевидно, что, если алгоритм не является численно устойчивым, то нет никакой гарантии, что ответ, найденный с его помощью им на ЭВМ, хоть как-то связан с входными данными. Нами изучена численная устойчивость двух известных алгоритмов решения систем уравнений с квазисепарабельными матрицами. Мы доказали численную устойчивость одного из алгоритмов и построили пример, в котором это свойство нарушается для другого алгоритма. Таким образом, мы прояснили какой из двух алгоритмов пригоден для практических вычислений на ЭВМ.

Апробация работы

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

кладной математической оптимизации и моделированию APMOD (Падер-борн, 2012), 16-ой (Пиза, 2010) и 17-ой (Брауншвейг, 2011) конференциях международного общества линейной алгебры ILAS, 24-ой международной конференции по численному анализу (Глазго, 2011), 3-ей международной конференции по матричным методам в математике и приложениях (Москва, 2011), конференции по линейной алгебре общества индустриальной и прикладной математики SIAM (Монтерей, 2009), а также на семинаре «Вычислительная математика и приложения» в ИВМ РАН, семинаре исследовательской группы по оптимизации ERGO Эдинбургского университета и семинаре «Ортогональность, теория аппроксимаций и приложения» в Мадридском университете имени Карлоса 3-го.

Публикации

По материалам диссертации опубликовано шесть работ, в том числе четыре (работы [1-4]) в журналах из Перечня ведущих рецензируемых научных журналов и изданий, рекомендованых ВАК, и две (работы [5,6]) в зарубежных рецензируемых сборниках.

Личный вклад автора

Результаты, описанные в главах 1 и 3 диссертации, получены автором самостоятельно. Результаты глав 2 и 4 получены в соавторстве, что отражено в соответствующей публикации [1]. Причем личный вклад автора в совместные работы основной — это теоремы 26 и 38 глав 2 и 4, соответственно, и результаты численных экспериментов.

Структура и объем работы

Диссертация состоит из введения, четырех глав основного текста и заключения. Объем диссертации — 135 страниц. Библиография включает в себя 96 наименований. Диссертация содержит 11 рисунков и 17 таблиц.

Краткое содержание работы

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

Матрицы с квазисепарабельной структурой были впервые упомянуты в фундаментальной работе Ф. Р. Гантмахера и М. Г. Крейна «Осцилля-ционные матрицы и малые колебания механических систем». Эти матрицы, в некотором смысле, являются дискретными аналогами функций Грина одномерных дифференциальных операторов. В качестве простейшего примера рассмотрим следующую краевую задачу с оператором Штурма— Лиувилля:

£(и):=(р(х)и')'-ч(х)и = Пх), Гй,(-и):=а,и(а) + (3,и'(а) =0, (1)

\Щи):=а2г1(Ъ) + р2и'(Ъ) =0,

В теории операторов Штурма-Лиувилля доказывается, что любое решение и(х) этой краевой задачи представимо в виде интеграла

U(x) =

g(x, £)f(£,)d£,

где д(х, £,) — так называемая функция Грина для (1). Причем функция Грина, для одномерного оператора Штурма-Лиувилля имеет вид

g(x,« = C-H*iUf!' Р)

[ui(£,)u2(x), £, < х < ъ,

где ui и u2 —ненулевые частные решения однородного уравнения £(гц) = 0 с Д(гц) = 0. Представление (2) означает, что функция д(х, £,) — сепа-рабельная в двух подобластях своей области определения. Дискретный аналог такой функции — это квазисепарабельная матрица, имеющая малоранговую структуру выше и ниже диагонали. Строгое определение квазисепарабельной матрицы приведено ниже.

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

max rankA(k + 1 : n, 1 : k) < р,

max rank Ail : k, k+ 1 : n) < q w

при фиксированных р и q .

Из этого определения вытекает, что квазисепарабельные матрицы задаются всего 0{п) независимыми параметрами. Сложность основных алгоритмов линейной алгебры с квазисепарабельными матрицами также 0{п), благодаря тому, что эти алгоритмы оперируют с параметрами, а не с элементами матрицы.

Квазисепарабельные матрицы находят свое применение во многих областях прикладной математики. В частности, в теории систем и теории управления, ортогональных многочленах и теории приближений, интегральных уравнениях и статистике. В главах 1 и 2 настоящей диссертации развиваются численные методы в теории ортогональных многочленов, основанные на использовании квазисепарабельных матриц. Глава 3 посвящена применению квазисепарабельных матриц к решению уравнений в частных производных и задач оптимального управления с ними. В главе 4 развивается теория квазисепарабельных матриц, а именно, анализ ошибок округления для алгоритмов с ними.

Первая глава посвящена вычислению корней многочлена, представленного в базисе произвольных ортогональных многочленов, таких, как многочлены Чебышёва, Лагранжа, Эрмита, Сегё и прочие. Пусть Р(х) и {rkM}£=0 — многочлен степени n и ортогональный базис соответственно, тогда

Р(х) = rn(x) + m-n-i ■ rn_, (х) + ... + тщ • ri (х) + ттю • г0(х).

(4)

Хорошо известно, что многочлены с единичным старшим коэффициентом, ортогональные на отрезке вещественной оси, удовлетворяют трехчленным рекуррентным соотношениям

Tk+i М = (х — ак) ■ rk(x) - ßk • гк_, (х), к = 0,1,2,....

(5)

В §1.2 мы показываем, что из (4) и (5) следует, что корни многочлена Р(х) совпадают с собственными значениями обобщенной матрицы Фробениуса

А =

00 ßi —то

1 а1 ß2 — mi i •• ••■ ;

*•• «n-2 ßn-1 — ТПП_2

1 «n-1 — ТПп-1

(6)

Одно из свойств этой матрицы — это то, что любая подматрица в ее верхнетреугольной части имеет ранг не выше 2, то есть А — квазисепарабель-ная. Казалось бы, при чем тут квазисепарабельность? Однако, в §§1.2-1.3

мы доказываем, что 1) факторы Ьи-разложения квазисепарабелыюй матрицы наследуют квазисепарабельную структуру, и что 2) эта структура сохраняется в итерациях ЬИ-алгоритма, то есть, при преобразованиях А в А', таких, что

А - сг1 = Ш, А' = Ш. (7)

Так как алгоритмы с квазисепарабельными матрицами обладают линейной по размеру сложностью, сложность одной ЬИ-итерации для матрицы из (6) также линейная.

Низкой алгебраической сложности недостаточно, чтобы считать алгоритм «хорошим», другая необходимая характеристика — численная устойчивость. К сожалению, ЬЯ-алгоритм (7) не является численно устойчивым. Однако, в 1994 году К. Фернандо и В. Парлетт получили аналог ЬК-алгоритма для трехдиагональных матриц, для которого известно, что в симметричном положительно определенном случае ошибка в вычисленных им собственных значениях Л{ мала относительно точных собственных значений Л^: ^

— Лг| ^ етах^ДО, 1=1,2,... (8)

Интересно, что даже такой широко известный численно устойчивый алгоритм, как <311, не обладает свойством (8).

В §1.4 мы обобщили алгоритм Фернандо и Парлетта на случай хессен-берговых квазисепарабельных матриц. Новый алгоритм использует все те же 0{п) арифметических операций на шаг. Так как нашей главной задачей было вычисление собственных значений матрицы А из (6) (корней многочлена (4)), в §1.6 мы уточнили новый алгоритм на этот случай. Хотя у нас нет доказательства, что в этом более общем случае предложенный алгоритм все ещё удовлетворяет (8), результаты численных экспериментов, представленные в §1.7, демонстрируют его более высокую точность в сравнении с алгоритмом СЗИ.

Во второй главе мы рассматриваем задачу об обращении полиномиальной матрицы Вандермонда

"г0(х,) г^х]) ■•• гп-Пх,)'

Уг= "'

.То(Хп) л(хп) ■■• Гп_! (хп)

Эта матрица обобщает хорошо известную обыкновенную матрицу Вандермонда, для которой гк(х) = хк. Задача об обращении матрицы Уг — классическая задача линейной алгебры, которая возникает, например, при интерполировании в базисе многочленов {^(х)}^ функции, заданной значениями в узлах {хк}£=1 •

(9)

Известно, что для узлов, равномерно распределенных на отрезке, классическая матрица Вандермонда плохо обусловлена. В 1994 году Е.Е. Тыр-тышниковым было доказано, что ее число обусловленности растет экспоненциально с размером. Однако в 1966 году Д. Траубом был получен алгоритм обращения матрицы Вандермонда за 0{п2) арифметических операций с малой прямой ошибкой. Казалось бы, это невозможно: прямая ошибка равняется обратной (не меньше машинного нуля) умноженной на число обусловленности (огромно!), однако весь «трюк» заключается в том, что алгоритм Трауба — структурный (0[п2) против 0(п3)), и к нему неприменим стандартный результат теории ошибок округления. Другими словами, структурное число обусловленности матрицы Вандермонда мало. Следуя успеху Трауба, многие исследователи в области численной линейной алгебры обратились к задаче обращения матрицы (9) для различных базисов многочленов. Л. Райхел и Г. Опфер обобщили алгоритм Трауба на случай многочленов Чебышёва, а В. Ольшевский — на случай многочленов Сегё.

Нами решена задача об обращении полиномиальной матрицы Вандермонда (9) для абсолютно произвольной системы многочленов (тк(х)}. Тем самым мы обобщили все результаты, полученные ранее для частных случаев. Наш подход основан на использовании междисциплинарной связи между тремя разделами науки: (1) электротехникой, (и) информатикой и (Ш) линейной алгеброй.

Из информатики мы берем так называемый принцип Теллегена, который мы рассматриваем в §2.2. В самом общем смысле принцип Теллегена гласит, что сложности умножения на квадратную матрицу без нулевых строк и столбцов и ее транспонированную равны. Доказательство этого факта основано на интерпретации умножения на матрицу через поток сигнала по взвешенному направленному графу.

В §2.3 мы связываем теорию графов с многочленами через граф потока сигнала, использующийся в электротехнике. Мы показываем, как обращение обыкновенной матрицы Вандермонда связано с обращением направления потока в таком графе. В следующем §2.4 мы обобщаем эту теорию на случай произвольных многочленов и полиномиальной матрицы Вандермонда и с помощью принципа Теллегена доказываем теорему 26, описывающую обратную к (9) матрицу. Из этой теоремы вытекает быстрый 0[п2) алгоритм обращения, который мы применяем в §2.5 к матрице, для которой аналогичный алгоритм не был известен ранее. Результаты численных экспериментов демонстрируют малую прямую ошибку в обратной матрице, найденной предложенным алгоритмом.

Третья глава посвящена задачам управления с ограничениями в виде уравнений в частных производных и методам их решения. В §3.2 рассмат-

ривается одна такая модельная задача, приведенная ниже

тт1||и-и||| + ЭМ2,

- у2и = f в П, и = д на ЗО.

Решение подобных задач — один из важных этапов современного инженерного проектирования.

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

(11)

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

2рМ 0 -М т "0"

0 М кт и - ъ

-М К 0 Л й

— (¿М + КМ -'к').

(12)

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

В §3.3 показывается, что плотная матрица (12) обладает малоранговой структурой особого вида, называемой многоуровневой квазисепарабелъ-ной (определение 31). В этом параграфе также объясняется как малопараметрическое представление матрицы Б из (12) может быть вычислено всего за 0[п) арифметических операций. Далее мы используем это структурное представление матрицы для построения ее предобусловливателя. Результаты численных экспериментов, представленные в §3.4, демонстрируют, что наш предобусловливатель обладает хорошими спектральными свойствами и, при использовании совместно с методом сопряженных градиентов, приводит к алгоритму решения (11) линейной по размеру задачи сложности.

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

Мы ставим вопрос о численной устойчивости алгоритмов решения систем уравнений с квазисепарабельными матрицами. Стандартный алгоритм решения системы уравнений — метод исключения Гаусса (LU-факто-ризация). Однако выбор ведущего элемента (перестановка строк и столбцов) в этом алгоритме разрушает квазисепарабельную структуру матрицы (а значит и линейную сложность), а без выбора ведущего элемента метод исключения Гаусса неустойчив. С другой стороны QR-факторизация — устойчивый алгоритм даже без перестановок строк и столбцов, и, в отличие от неструктурного случая, имеет линейную 0{п) сложность для квазисепарабельных матриц. Мы рассматриваем два алгоритма решения системы уравнений с квазисепарабельной матрицей через QR-фактори-зацию последней, представленные в научной литературе.

Алгоритм A: Eidelman Y. and 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.

Алгоритм Б: Van Camp E., Mastronardi N. and Van Barel M. Two fast algorithme for solving diagonal-plus-semiseparable linear systems // Journal of Computational and Applied Mathematics. 2004. V. 164. P. 731-747.

В §4.2 мы описываем общие идеи этих алгоритмов. Первый из алгоритмов детально рассматривается в §4.3. Главный результат этого параграфа — теорема 41 о численной устойчивости алгоритма А, приведенная ниже.

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

(А + ДА)х = b + ДЪ,

где

||АА(:,))||2<Тк,п2||А(:о)||2) j = 1,...,n, ||ДЪ||2 < Тк2п|№

и К], Кг — малые (натуральные) константы, у^ — общепринятое обозначение для накопившихся ошибок округления. Обозначим, через и машинный ноль, тогда у*. := ки/(1 — ки).

В §4.4 мы проводим анализ ошибок округления для алгоритма Б. Нами построен пример размера 4x4, который демонстрирует, что алгоритм Б численно неустойчив. Мы объясняем причины неустойчивости и подтверждаем наши выводы результатами численных экспериментов.

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

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

1. Выведен новый алгоритм вычисления собственных значений хессен-берговых матриц с малоранговой структурой в верхней части, обобщающий с^ск-алгоритм Фернандо и Парлетта для трехдиагональ-ных матриц. На его основе построен новый метод вычисления корней многочлена, разложенного по произвольному ортогональному базису, сложности 0{п) на корень.

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

3. Определен и изучен новый класс матриц с малоранговой структурой — многоуровневые квазисепарабельные матрицы.

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

5. Доказана теорема о численной устойчивости алгоритма (ЗИ-факто-ризации квазисепарабельной матрицы и основанного на нем метода решения систем уравнений.

Публикации по теме диссертации

1. Olshevsky V., Tyrtyshnikov Е., Zhlobich P. Tellegen's principle, non-minimal realizations of systems and inversion of polynomial Van-dermonde matrices // Russian Journal of Numerical Analysis and Mathematical Modelling. 2012. V. 27, № 2. P. 131-154.

2. Bella Т., 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.

3. 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.

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

5. Bella Т., Eidelman Y., Gohberg I. et al. Classifications of recurrence relations via subclasses of (H,m)-quasiseparable matrices // «Numerical Linear Algebra in Signals, Systems and Control» / Ed. by P. Van Dooren, S. Bhattacharyya, R. Chan et al. — Springer, 2011. — V. 80 of «Lecture Notes in Electrical Engineering». — P. 23-53.

6. Bella Т., Eidelman Y., Gohberg I. et al. A Traub-like algorithm for Hessen-berg-quasiseparable-Vandermonde matrices of arbitrary order // «Numerical Methods for Structured Matrices and Applications» / Ed. by V. Mehrmann, D. Bini, V. Olshevsky et al. — Birkhauser Basel, 2010. — V. 199 of «Operator Theory: Advances and Applications». — P. 127-154.

Отпечатано в типографии МГУ 119991, ГСП-1, г. Москва, Ленинские горы, д. 1, стр. 15 Заказ № 1077. Тираж 70 экз.

Введение

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. Изучена численная устойчивость двух алгоритмов решения системы уравнений с квазисепарабельной матрицей. Для одного из алгоритмов доказана теорема о численной устойчивости, в то время, как для другого построен контрпример, демонстрирующий численную неустойчивость, и объяснены ее причины.

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

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.