«матриц ...»

Российская академия наук

Институт вычислительной математики

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

Оселедец Иван Валерьевич

УДК 519.6

Нелинейные аппроксимации матриц

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

ДИССЕРТАЦИЯ

На соискание учёной степени

кандидата физико-математических наук

Научный руководитель чл.-корр. РАН, проф. Тыртышников Е. Е.

Москва 2007 1 Содержание Введение 2 i.1 Нелинейные аппроксимации матриц: зачем и как.... 3 i.2 Основные результаты работы................ i.3 Содержание работы по главам............... Глава 1. Метод чёрных точек и наилучшие циркулянтные предобуславливатели 1.1 Введение............................ 1.2 Задачи C+R и D+R аппроксимации............ 1.3 Чёрные точки, малые ранги и скелетоны......... 1.4 Адаптивная версия метода чёрных точек......... 1.5 Тёплицев случай........................ 1.5.1 Быстрое вычисление образа Фурье для тёплицевой матрицы...................... 1.6 Существование C+R аппроксимации для некоторых классов тёплицевых матриц.................... 1.7 Численные эксперименты.................. 1.8 Метод чёрных точек для произвольного шаблона.... 1.9 Неизвестный шаблон..................... 1.10 Выводы............................. Глава 2. Нестандартные вейвлет-преобразования 2.1 Введение............................ 2.2 Основные понятия и определения.............. 2.3 Вейвлет-пространство. Масштабирующие и лифтинговые коэффициенты....................... 2.4 Основная система....................... 2.5 Решение основной системы.................. 2.6 Нахождение масштабирующих коэффициентов...... 2.7 Алгоритм вычисления вейвлет-преобразования...... 2.8 Численные эксперименты.................. 2.8.1 Пример 1........................ 2.8.2 Пример 2........................ 2.9 Выводы............................. Глава 3. Тензорные аппроксимации матриц со структурированными факторами 3.1 Введение............................ 3.2 Масштабированные циркулянтные предобуславливатели 3.3 Приближённое обращение структурированных матриц. 3.4 Методы построения приближённой обратной матрицы............................ 3.4.1 Метод Ньютона с аппроксимациями........ 3.4.2 Модифицированный метод Ньютона........ 3.5 Численные результаты.................... 3.5.1 Масштабированный циркулянтный преобуславливатель........................ 3.5.2 Предобуславливатели на основе метода Ньютона 3.6 Выводы............................. Глава 4. Супер-быстрое обращение двухуровневых тёплицевых матриц 4.1 Введение............................ 4.2 TDS формат.......................... 4.3 Арифметика TDS формата.................. 4.3.1 Основные арифметические операции........ 4.4 Основные арифметические операции в тензорном формате 4.4.1 TDS-рекомпрессия.................. 4.4.2 Оператор обрезания................. 4.5 Метод Ньютона и выбор начального приближения.... 4.6 Численные результаты.................... 4.7 Структура обратных к двухуровневым матрицам специального вида.......................... 4.7.1 Так почему же 5?................... 4.7.2 Обобщение на случай большего числа слагаемых 4.8 Выводы............................. Заключение Литература Введение i.1. Нелинейные аппроксимации матриц: зачем и как К решению линейных систем уравнений основной задаче линейной алгебры и матричного анализа сводится подавляющее большинство практических вычислительных задач. Однако, несмотря на наличие универсальных методов, многие приложения приводят к большим системам, которые не могу быть решены даже на современных суперкомпьютерах.

В данной диссертации развиваются эффективные методы вычислений с плотными матрицами, в которых сами матрицы и результаты матричных операций аппроксимируются матрицами специальной структуры, определённой относительно малым числом параметров.

Зависимость от параметров носит нелинейный характер, поэтому речь идёт о методах нелинейной матричной аппроксимации.

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

Часто структура матрицы видна сразу или следует из физических свойств задачи. Например, в задачах с оператором, инвариантным относительно сдвига, получающиеся матрицы имеют тёплицеву (или блочно-тёплицеву в многомерных задачах) структуру, т.е. элемент матрицы зависит лишь от разности индексов: aij = bij. Для тёплицевых матриц существуют быстрые алгоритмы, основанные на БПФ. Тёплицевы матрицы классический пример матриц с линейной структурой. Можно привести другие примеры: ганкелевы матрицы ( aij = bi+j ), ленточные матрицы, разреженные матрицы.

Ещё один важнейший класс матриц матрицы малого ранга, т.е. матрицы вида U Rnr, V Rmr, где ранг r m, n. Это пример матрицы с нелинейной структурой: её элементы зависят от параметров (элементов матриц U и V ) нелинейно.

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

Тёплицевы матрицы соответствуют одномерным интегральным уравнениям, где использование сеток большой размерности не является необходимым. На практике значительно более интересным представляется решение многомерных уравнений. Для ядер, инвариантных относительно сдвига и дискретизации на равномерной сетке, получаются многоуровневые тёплицевы матрицы. Такие матрицы тоже можно умножать на вектор за квазилинейное время, однако до сих пор универсальных прямых методов решения таких систем за то же время неизвестно. Существующие формулы (формулы Гохберга-Хайнига) содержат не O(N) параметров, а O(N3/2 ), и, видимо, удобных формул с меньшим числом параметров не существует. Что же делать?

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

i.2. Основные результаты работы При решении любой задачи всегда хочется получить сначала некоторый общий подход. В данной диссертации предложены два таких общих подхода для двух классических задач матричного анализа обращение матриц и построение предобуславливателей. Напомним, о чём идёт речь. Пусть нам нужно решить линейную систему с помощью какого-нибудь стандартного итерационного метода. Однако часто бывает так, что требуется большое количество итераций для сходимости, поэтому решают эквивалентную систему вида где матрица P легко обратима, (т.е. P1 f можно вычислить достаточно быстро). Матрица P называется предобуславливателем. Как проверить, что матрица P хорошая? Обычно хотят добиться того, чтобы матрица AP1 была хорошо обусловлена. Однако задача оптимизации числа обусловленности является довольно сложной, и поэтому её заменяют гораздо более простой задачей аппроксимации вида где P принадлежит некоторому классу быстрообратимых матриц, а некоторая (обычно фробениусова) матричная норма. Такой подход, например, применяется для построения циркулянтных предобуславливателей для тёплицевых матриц1. Однако, как известно из теории итерационных методов, хорошая обусловленность достаточна, но отнюдь не необходима для быстрой сходимости итерационных методов. Большую роль играет наличие кластеров собственных значений предобусловленной системы около 1. Это означает, что подавляющее большинство собственных значений, за исключением может быть конечного числа, находится в окрестности 1. А как проверить наличие кластеров? Оказывается, что все теоремы о существовании кластеров явно или неявно опирается на представление матрицы A в виде где P предобуславливатель, R матрица малого ранга и E матрица малой нормы, т.е. выполняется приближённое равенство где R поправка малого ранга. Наше предложение (и общий подход!) состоит в том, чтобы использовать (2), а не (1) в качестве отправной точки для построения предобуславливателя P. Если мы зафиксируем класс предобуславливателей (например, циркулянтные матрицы), то мы получим некоторую задачу нелинейной матричной аппроксимации, в которой необходимо находить и P, и R. Обычно находят P, а R оценивают; однако, как увидим позднее, гораздо более эффективно находить P и R одновременно. Выбирая различные классы матриц P, мы получаем целую россыпь новых задач нелинейной матричной аппроксимации, для которых можно пытаться придумать эффективные алгоритмы решения. В данной диссертации рассмотрены следующие Такие предобуславливатели называются предобуславливателями Т. Чэна (T. Chan) классы матриц, в которых ищутся предобуславливатели: циркулянтные матрицы, разреженные матрицы с известным шаблоном, матрицы вида где T матрица какого-либо быстрого преобразования (Фурье, синус-преобразование, косинус-преобразование, вейвлет-преобразование), аS разреженная матрица. Идея построения основана на методе чёрных точек, который является обобщением крестового метода [?, 40, 13] для приближения матрицами малого ранга. Фактически, с помощью метода чёрных точек можно для заданной матрицы A построить за число операций порядка O(N) приближение (если оно, конечно, существует) вида где S разреженная матрица с известным шаблоном, R матрица малого ранга. Область применимости метода не ограничивается только предобуславливанием. Нетрудно увидеть в (4) задачу приближения матрицы A матрицей малого ранга за исключением малого числа элементов. Такая задача возникает при заполнении пропусков в больших массивах данных, например данных наблюдений или данных, связанных с исследованием ДНК. Также в диссертации предлагается обобщение метода на случай, когда положение испорченных элементов неизвестно и требует определения.

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

Нетрудно показать, что Xk A1 с квадратичной скоростью (для всех начальных приближений X0, достаточно близких к A1. Для обычных матриц такой алгоритм, конечно, непрактичен сложность его составляет O(N3 ) операций на одну итерацию.

Однако для структурированных матриц, допускающих быстрое умножение, метод оказывается очень эффективен. При этом возникают дополнительные сложности, связанные с восстановлением структуры матрицы после каждой итерации. В итоге, итерационный процесс можно записать в виде где R некоторый (нелинейный) оператор проектирования. В диссертации показано, что при достаточно общих предположениях на оператор R, метод сохраняет квадратичную скорость сходимости. Оказывается, что для структурированных матриц можно построить модификацию метода Ньютона, которая оказывается гораздо более быстрой; в точной арифметике два итерационных процесса эквивалентны, а в приближённой модифицированный метод требует умножения матриц гораздо более простой структуры, что приводит к уменьшению вычислительных затрат.

Используя различные форматы данных, теперь можно получать различные методы. Для многоуровневых матриц удобным форматом оказываются матрицы малого тензорного ранга:

где U V блочная матрица вида Выгода от такого представления очевидна вместо N2 = n4 параметров требуется всего лишь rn2 = N параметров, а r обычно порядка 10-20. Как же вычислять такое представление? Оказывается, что после простой перестановки индексов задача сводится к задаче аппроксимации переставленной матрицы матрицей малого ранга. Для решения этой задачи особенно эффективен метод неполной крестовой аппроксимации, который и применяется в дальнейшем. После этого, можно поставить вопрос о дальнейшем сжатии матриц-факторов Ui Vi. Это необходимо, так как умножение на вектор в тензорном формате всё ещё требует более чем линейное по N число операций O(N3/2 ). Это уже не очень приятно при N = 10000. Возможно несколько подходов.

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

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

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

Однако оценок рангов смещения для факторов получить не удалось.

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

При исследовании вопроса о структуре матриц, обратных к дважды тёплицевым был (сначало экспериментально) обнаружен интересный факт. Пусть матрица A имеет вид где R – матрица ранга 1, а D диагональная матрица. Тогда тензорный ранг A не превышает 5. Это утверждение удаётся обобщить на матрицы вида Доказательство является конструктивным и даёт способ представления обратной матрицы. То есть, обнаружен класс матриц малого тензорного ранга, замкнутый относительно обращения.

Как это связано с дважды тёплицевыми матрицами? Обращение такой матрицы мы начинаем с аппроксимации матрицей малого тензорного ранга с факторами вида C + R. Из вышеозначенных результатов вытекает, что соответствующая обратная матрица имеет относительно малый тензорный ранг. Таким образом, получен теоретически обоснованный алгоритм для обращения дважды тёплицевых матриц линейной сложности. Если бы у нас получилась оценка на ранг смещения для факторов, то мы получили бы полностью обоснованный алгоритм сублинейной сложности (по порядку зависимости от размера матрицы). Заметим, что такая сублинейная зависимость наблюдалась нами на модельных задачах.

i.3. Содержание работы по главам Первая глава посвящена классической задаче теории структрированных матриц построение циркулянтных предобуславливателей для тёплицевых (и не только!) матриц. И тёплицевы, и циркулянтные матрицы матрицы линейной структуры и в многочисленных предыдующих работах использовались лишь методы на основе наилучших (в некоторой норме) линейных приближений циркулянтами заданной тёплицевой матрицы. Мы же сформулировали задачу, как задачу нелинейной аппроксимации заданной матрицы суммой циркулянта и матрицы малого ранга. Раздел 1.1 содержит краткое описание исходной задачи и состояния дел на данный момент. Там же сразу формулируется основная задача, которую мы будем решать задача C+R аппроксимации. В разделе 1.2 даются различные формулировки задачи C + R аппроксимации и показывается, как эта задача связана с восстановлением неизвестных элементов в малоранговой матрице (матрице с чёрными точками ). В разделе 1.3 формулируется алгоритм чёрных точек для решения задачи D + R аппроксимации, к которой сводится задача C + R аппроксимации. Этот алгоритм не является итерационным и доказана теорема о том, что алгоритм восстанавливает подавляющее большинство пропусков за конечное число итераций, на практике порядка 10-20. В разделе 1.4 формулируется практическая, адаптивная версия метода чёрных точек, позволяющая строить C+R и D+R аппроксимации без какой-либо дополнительной информации на вход надо лишь подать исходную матрицу и нужный параметр точности аппроксимации. Для матриц общего вида, описываемых n2 параметрами, получается алгоритм сложности O(n2 ).

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

В разделе 1.6 приведены теоретические результаты, касающиеся существования C + R аппроксимаций для класса тёплицевых матриц.

Оказывается, такие аппроксимации существуют для практически всех матриц, встречающихся в литературе по супер-линейному предобуславливанию. В разделе 1.7 приведены численные эксперименты по построению C+R аппроксимаций для различных тёплицевых матриц.

В разделе 1.8 идея метода чёрных точек получает своё естественное продолжение. Приведён вариант метода, позволяющий восстанавливать неизвестные элементы в малоранговой матрице с произвольным расположением этих самых неизвестных элементов. Однако, в отличие от диагонального шаблона, надо внимательно следить за тем, какие элементы удалось восстановить, а какие нет. Для решения этой проблемы предложены два способа. В разделе 1.9 сделано самое общее возможное обобщение метода чёрных точек на случай, когда положение неизвестных элементов неизвестно. Для этого предложено максимизировать разреженность матрицы A R с помощью минимизации специального функционала.

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

В разделе 2.2 даются основные понятия и определения, необходимые для дальнейшего изложения. В разделе 2.3 вводится самый важный в главе объект вейвлет-пространство и описывается так называемая лифтинговая схема построения вейвлет-преобразований с требуемыми свойствами. В разделе 2.4 формулируются основные требования на вейвлет-преобразование: наличие заданного количества нулевых моментов и выписывается система линейных уравнений специального вида на коэффициенты, определяющие искомое преобразование. В разделе 2.5 формулируется основной результат главы и выписывается явная формула для решения основной системы. Раздел 2.6 посвящён нахождению масштабирующих коэффициентов показано, что их можно находить с помощью уже описанного алгоритма по аналогичным формулам. В разделе 2.8 описан конкретный, пошаговый способ реализации вейвлет-преобразования, требующий O(n) операций.

Также описаны алгоритмы вычисления обратного и обратного транспонированного преобразования они активно используются в численных расчётах для восстановления исходных данных по преобразованным. В разделе 2.8 приведены численные эксперименты, сравнивающие новые преобразования с преобразованиями Добеши. Показано, что выигрыш по степени сжатия составляет в различных примерах от 30% до 50% процентов.

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

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

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

Глава 4 посвящена теоретическому и практическому изучению вопроса обращения двухуровневых тёплицевых матриц. В рамках развиваемого подхода построен построен метод Ньютона с аппроксимациями для обращения двухуровневых тёплицевых матриц с использованием введённого TDS-формата (Tensor-displacement-structure). Сам новый формат вводится в разделе 4.2. В разделе 4.3 описываются основные арифметические операции над матрицами в TDS формате сложение, умножение, и показывается, что их можно выполнить за O( n) log n) операций. Важнейший элемент одного шага метода Ньютона оператор обрезания описан в том же разделе. Показано, что задача сводится к вычислению фробениусова скалярного произведения двух TDS-матриц, и это вычисление может быть проведено очень быстро. В разделе 4.5 приводится напоминание о работе метода Ньютона и описывается способ выбора начального приближения. В разделе 4.6 приведены численные эксперименты. И, наконец, в разделе 4.7 впервые получены теоретические результаты о структуре обратных матриц к матрицам малого тензорного ранга специального вида.

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

Глава 1.

Метод чёрных точек и наилучшие циркулянтные предобуславливатели 1.1. Введение Начнём мы с рассмотрения одной классической задачи построения циркулянтных преобуславливателей и покажем, как она сводится к задаче нелинейной матричной аппроксимации.

Идея использования циркулянтов в качестве предобуславливателей к тёплицевым матрицам была впервые предложена Гильбертом Стрэнгом в 1986 году [31]. Его идея состояла в том, чтобы строить циркулянт, используя половину элементов из первого столбца и первой строчки тёплицевой матрицы. Другой популярный подход так называемый оптимальный предобуславливатель Т. Чэна (T. Chan) ближайший во фробениусовой норме циркулянт к заданной (тёплицевой) матрице [6]. Эти предобуславливатели легко построить, но в некоторых случаях они не работают (число итераций может сильно расти с увеличением размера матрицы n ). Для плохих случаев было построено несколько методов и алгоритмов (см. обзор [8]). Однако, наиболее эффективные подходы явно используют информацию о так называемом символе (производящей функции) тёплицевой матрице (чуть ниже мы дадим определение, что такое символ). По этой причине они могут быть названы функциональными, а не матричными (см. [25]) подходами. Более того, метод, подходящий для симметричных положительно определённых матриц, может не работать для незнакоопределённых или несимметричных матриц существующие подходы построения циркулянтных преодобуславливателей не являются универсальными.

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

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

Начнём с самого начала. Если линейная система решается с помощью какого-либо итерационного метода (такого, как CG или GMRES) и наблюдается медленная сходимость (что случается часто), тогда известное лекарство состоит в переходе к предобусловленной системе где P называется предобуславливателем. Анализ качества предобуславливателя обычно начинается с погружения выбранной системы в последовательность систем (матриц коэффициентов, правых частей, предобуславливателей), параметризованных размером матрицы n. Далее, для того, чтобы предобуславливатель был хорошим, добиваются выполнения следующих свойств (a) Для AP1 существует равномерная оценка на число обусловленности по n ;

(b) Собственные значения AP1 имеют кластер в 1. Это, значит, что подавляющее большинство собственных значений матрицы AP1 находятся в –окрестности 1.

По крайней мере, для эрмитовых положительно определённых матриц и при некоторых дополнительных предположениях в общем случае, свойство (a) означает линейную сходимость, а свойство (b) дат так называемую суперлинейную сходимость (см. [42]). Поэтому свойство (b) особенно интересно. Когда же существует кластер и как доказывается его существование? Существование кластера тесно связано с разложением вида [39] где rank R = r n и ||E||. Отметим, что матрицы в правой части (1.1) зависят от n и.

Но если представление (1.1) так привлекательно, то почему бы не исходить прямо из него? Мы предлагаем строить предобуславливатели P основываясь прямо на (1.1). Если мы будем выбирать P из некоторого матричного класса, то (1.1) становиться задачей аппроксимации со следующей (пока нестрогой) формулировкой:

C+R аппроксимация Аппроксимировать матрицу A, суммой вида где C циркулянт, (в (1.1) соответствует P ) и R матрица малого ранга.

Рассмотрим, например, тёплицевы матрицы A = [aij ] размера n = 128, 256, 512, порождённые символом f = x4 (это означает, что ak являются коэффициентами Фурье для f ). Пусть P = C в (1.1) является либо предобуславливателем Стрэнга, либо предобуславливателем Т.

Чэна. Тогда, установив точность = 102, найдём R с помощью отбрасывания сингулярных чисел матрицы по порогу так, что ||E||2. В этом случае получаются следующие ранги для R :

Таблица 1.1. Зависимость rank R от n ( = 102 ) Для этого фиксированного, rank R практически не зависит от размера матрицы n. Однако исследуем, как этот ранг зависит от для фиксированного n.

Таблица 1.2. Зависимость rank R от ( n = 256 ).

Как мы видим, оба предобуславливателя сработали неудовлетворительно: кластера собственных значений нет собственные значения не групируются около единицы. Ситуация с предобуславливателем Стрэнга выглядит получше, однако видно, что ранги растут как, 1. Матрица A плохо обусловлена, поэтому мы должны аппроксимировать её с высокой точностью и ни один из двух рассмотренных предобуславливателей не даёт собственный кластер. Однако данная матрица может быть очень точно аппроксимирована суммой циркулянта и матрицы достаточно малого ранга. Но соответствующий циркулянт не имеет ничего общего ни с преобуславливателем Стрэнга, ни с предобуславливателем Т. Чэна. Более того, будет доказано, что для довольно общего класса тёплицевых матриц (включая все примеры в статьях по суперлинейным предобуславливателям) существуют аппроксимации суммой циркулянта и матрицы малого ранга[63] с оценкой вида Поэтому мы можем быть уверены, что хороший циркулянт существует и заинтересованы в том, как его вычислить.

1.2. Задачи C+R и D+R аппроксимации Так как каждый циркулянт диагонализуется с помощью дискретного преобразования Фурье (Discrete Fourier Transform DFT) где F это матрица дискретного преобразования Фурье, а D диагональная матрица из собственных значений, то задача C + R аппроксимации может быть переформулирована следующим образом:

Поэтому, общая задача C + R аппроксимации легко сводится к задаче D + R аппроксимации, где D диагональная матрица. В силу унитарности матрицы Фурье F задачи D + R и C + R аппроксимации эквивалентны.

Теперь давайте строго определим, что же такое приближённо и малый ранг. Если мы зафиксируем ранг малоранговой части разложения, то получится следующая задача аппроксимации:

D+R I: Для заданной матрицы A и целого числа r > 0 найти матрицу B вида B = D + R, где rank R r,и диагональную матрицу D, которые минимизируют ||A B||F.

Мы можем исключить из формулировки либо R либо D. Если мы исключим R, то получится оптимизационная задача в терминах сингулярных чисел.

D+R II: Для заданной матрицы A и целого числа r > 0 найти диагональную матрицу D, которая минимизирует Важно отметить, что это негладкая, невыпуклая задача оптимизации, которая, по-видимому, имеет много локальных минимумов (нам не известно о каких либо эффективных методах её решения).

Если же мы исключим D, то получится следующая формулировка:

D+R III: Для заданной матрицы A = [Aij ] и целого числа r > найти матрицу R = [Rij ] ранга не выше r, которая минимизирует Задача D + R аппроксимации была впервые рассмотрена в [3], где был предложен итерационный метод её решения. Это был вариант метода переменных направлений, названный ADR (Alternating Diagonal Rank) с двухшаговой итерацией следующего вида:

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

1. D = arg min ||A D R||F;

Видно, что на каждом шаге невязка ||A D R||F невозрастает. К сожалению, по-видимому, это единственное достоинство метода ADR.

Часто для его сходимости требуется большое количество итераций.

Иногда он застревает в локальном минимуме. Он имеет большую вычислительную сложность делает его (в вышеприведённом виде) неприемлемым для практического использования. Существует способ видоизменить этот метод (совсем нетривиально) так, чтобы он сходился к глобальному минимуму, однако вычислительная сложность всё равно остаётся очень высокой. Тем не менее, если имеется хорошая аппроксимация полученная с помощью какого-нибудь другого метода (например, с помощью алгоритма, предлагаемого в данной диссертации), можно попытаться построить некоторую быструю версию ADR для улучшения заданной аппроксимации к решению.

Давайте посмотрим более внимательно на формулировку D + R III. Вспомним, что нам интересен случай, кого матрица A хорошо аппроксимируется суммой диагональной матрицы и матрицы ранга r.

Поэтому начнём с предположения, что матрица A является точной суммой диагональной матрицы и матрицы ранга r. Как можно восстановить D и R, зная только их сумму? Ответ мы узнаем совсем скоро.

1.3. Чёрные точки, малые ранги и скелетоны Задача формулируется следующим образом. Пусть матрица A является точной суммой диагональной матрицы и матрицы ранга r.

Зная, что A = D + R, как восстановить D и R по A ?

Очевидно, что в матрице R мы знаем все внедиагональные элементы. Поэтому, всё что нужно определить это диагональные элементы матрицы R. Перед описанием общего метода, рассмотрим следующий простой пример матрицы ранга 2 и размера 6 6 :

(Это действительно матрица ранга 2, так как aij = i + j ).

Предположим теперь, что мы не знаем диагональных элементов Диагональные элементы обозначены чёрными точками. Вопрос такой: как дополнить внедиагональную часть, заменив чёрные точки числами так, чтобы получившаяся матрица имела ранг 2? Можно применить следующую простую идею. Возьмём подматрицу, образованную столбцами 4,5,6 и строками 2,3,4:

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

Как нетрудно видеть, выбор этих строк и столбцов неединственнен. В данном примере, различные способы выбора опорных строк и столбцов приведут к одному и тому же результату. На практике же, выбор правильных строк и столбцов, используемых для восстановления очень важен, и неправильный выбор может привести (и часто приводит) к неустойчивости.

Опишем теперь приведённую процедуру в общем случае и докажем, что она действительно восстанавливает чёрные точки.

Рассмотрим произвольную матрицу B ранга r, возьмём r линейно независимых строк и r линейно независимых столбцов из B и образуем матрицы L Rnr (из столбцов) и U Rrn (из строк). Пусть B обозначает подматрицу размера r r, находящуюся на пересечении этих выбранных строк и столбцов. Тогда подматрица B невырождена и матрица B может быть представлена в виде который называется скелетным разложением.

Основное утверждение состоит в том, что матрица ранга r однозначно определяется по своим линейно независимым r столбцам и строкам. Построим теперь скелетное разложение для матрицы с чёрными точками на диагонали. Для нашего примера, строки 3,4 и столбцы 1,2 дают нам невырожденную подматрицу на пересечении, и поэтому мы можем написать Для того, чтобы определить, как работать с чёрными точками, введём следующую арифметику чёрных точек :

где x обычное число. Поэтому, умножение матриц в (1.3) даёт Это означает, что мы нашли (подчёркнутые) диагональные элементы (5,5) и (6,6). Так как все внедиагональные элементы в A заданы, теперь нам известны два полных столбца с номерами 5, 6 и две полных строчки с номерами 5, 6 из исходной A и можно снова использовать скелетное разложение для восстановления оставшихся неизвестных элементов полной матрицы.

В общем случае мы делаем то же самое.

Метод чёрных точек. Пусть задана матрица A, допускающая представление в виде A = D + R с диагональной матрицей D и матрицей R ранга r. Тогда для нахождения по крайней мере n 2r столбцов и строк неизвестной матрицы R нужно действовать следующим образом:

1. Взять в A невырожденную подматрицу A размера r r, которая не содержит диагональных элементов A. Пусть строчки и столбцы этой подматрицы имеют в A индексы i1,..., ir и j1,..., jr соответственно, и пусть матрицы L и U размерами n r и r n составлены из этих строк и столбцов.

2. Образовать матрицу Q = LA1U (всё ещё с чёрными точками) и обнаружить, что элементы больше не являются чёрными точками. Следовательно, к этому моменту нам известны по крайней мере n 2r диагональных элементов матрицы R.

Алгоритм основан на следующей простой теореме.

Теорема 1 Элементы определённой выше матрицы Q удовлетворяют (1.4).

Доказательство. Используя определение L и U, получаем, что Если i = j1,..., jr и j = i1,..., ir, тогда ни один из элементов вида Ri,jk, Rj,il не находится на главной диагонали. Поэтому, все эти элементы известны и соотвествущие элементы Q должны совпадать с соответствующими элементами R. В наших приложениях r n, поэтому первые два шага метода чёрных точек позволяют восстановить большинство диагональных элементов. Для нахождения оставшихся элементов R нужен ещё один шаг:

(3) Если n достаточно велико ( n 2r r, или, что эквивалентно, n 3r ), тогда нам известно по крайней мере r столбцов и r строчек матрицы R. Предположим, что эти r столбцов и строчек линейно независимы. Тогда, используя их для построения скелетного разложения, восстановим оставшиеся чёрные точки в R.

Отметим, что это шаг основан на предположении, что первые два шага дали нам r линейно независимых столбцов и строчек с известными элементами. Для этого достаточно потребовать, что в A есть две невырожденные подматрицы размера r r, не содержащие общих столбцов и строчек и не содержащие также диагональных элементов матрицы A.

1.4. Адаптивная версия метода чёрных точек Осталось несколько необсуждавшихся проблем. Во-первых, матрица A может не являться точной суммой диагональной матрицы и матрицы малого ранга. Вместо этого, пусть где E такое, что ||E|| может рассматриваться как некий шум. Более того, значение r ранга R (зависящее от ) может быть неизвестно заранее. Поэтому мы получаем задачу, где ранг r нужно определить, исходя из заданного порога точности.

В случае наличия шума, выбор столбцов и строчек (или, что одно и то же, выбор опорной подматрицы), по которым строится скелетное разложение, играет основную роль. Вопрос в том, какая подматрица является наилучшей. Если бы чёрных точек не было, ответом является так называемый принцип максимального объёма [13]: если подматрица A имеет максимальный объём (т.е. максимальный по модулю определитель) среди всех подматриц размера r r, тогда оценка на ошибку скелетного разложения выглядит следующим образом:

где r+1 (A) (r + 1) сингулярное число матрицы A. Мы предполагаем, что при наличии чёрных точек такой же хорошей подматрицей будет подматрица максимального объёма среди всех подматриц без чёрных точек.

Так как нахождение подматрицы максимального объёма составляет нетривиальную задачу, мы можем попытаться заменить её на подматрицу с достаточно большим объёмом. Такая подматрица может быть получена с помощью различных вариантов метода неполной крестовой аппроксимации (см. [14, 38, 1]).

Приведём описание крестового метода. Он работает следующим образом (здесь R это ненулевая матрица, которую нужно аппроксимировать малоранговой матрицей):

1. Инициализация: k = 0, Rk = R.

2. Выбор ведущего элемента: (i0, j0) = arg maxi,j |Rk |.

3. Вычислить крест-скелетон:

где uk это строка с номером i0 матрицы Rk и vk столбец с номером j0 в матрице Rk.

4. Вычислить новую невязку Rk+1 = Rk Ck.

5. Если норма невязки ||Rk+1|| достаточно мала, тогда остановиться и вернуть i=0 C как аппроксимацию ранга r к матрице R.

Иначе, увеличить k на 1 и перейти к шагу 2.

На каждом шаге мы вычитаем из матрицы одну матрицу ранга 1, называемую скелетоном. Она определяется по выбранным столбцу и строчке так, чтобы получившийся скелетон совпадал с матрицей Rk по одному столбцу и одной строчке. Это некоторый дискретный аналог интерполяции.

Так как мы интерполируем на каждом шаге одну строчку и один столбец, ведущие элементы элементы (т.е. элементы Rk0 j0 ) нужно выi бирать так, чтобы исключить большие элементы в матрице невязки. Когда алгоритм завершается, он на самом деле возвращает вариант скелетного разложения с предположительно хорошей подматрицей (выбор ведущего элемента обычно приводит к подматрице достаточно большого обьёма).

Суммарная стоимость приведённого варианта составляет O(n2 r) из-за того, что на каждом шаге ищется максимальный элемент по всей матрице (полное пивотирование). Изначально, однако, метод неполной крестовой аппроксимации был создан в надежде на то, что он может аппроксимировать матрицу малого -ранга, используя только лишь малое число её элементов [14]. Если матрица имеет точный ранг r, тогда гауссово исключение с выбором ведущего элемента даёт нулевой ведущий элемент точно после r шагов. На самом деле, точно такой же подход может быть применён и при наличии шума. Однако, можно реализовать различные стратегии методов выбора ведущих элементов (например, выбор по строке или столбцу). Особую вычислительную выгоду дают методы, использующие лишь небольшое число элементов матрицы. С частичным выбором ведущего элемента мы можем, конечно, получить больший коэффициент перед r+1 (A) в оценке (1.5); он зависит от того, как близок объём получающейся подматрицы к подматрице максимального объёма [13]. Для некоторого класса матриц, этот коэффициент может составлять 2r вместо r +1. В любом случае, даже не смотря на потерю точности (которая компенсируется увеличением r ), существенное снижение вычислительной сложности делает частичный выбор ведущего элемента единственным практически эффективным алгоритмом для метода неполной крестовой аппроксимации.

Метод неполной крестовой аппроксимации может быть легко приспособлен для случая матриц с неизвестными элементами (чёрными точками). Мы только должны следить за тем, как чёрные точки распространяются на каждом шаге. Строки и столбцы, содержащие чёрные точки, составят наши чёрные списки (обновляемые на каждом шаге).

Адаптивный метод чёрных точек.

1. Инициализация:

2. Найти ведущий элемент:

3. Вычислить крест-скелетон:

4. Вычислить новую невязку: Rk+1 = Rk Ck.

5. Добавить элемент i0 к Lc и элемент j0 к Lr.

6. Вычислить ошибку:

7. Если k достаточно мало, выйти и возвратить как аппроксимации к соответствующим диагональным элементам Rii. Иначе, увеличить k на 1 и перейти к шагу Приведённый алгоритм на самом деле возращает вариант скелетного разложения с некоторой адаптивно выбранной подматрицей R, размер которой не был задан заранее и ожидается, что это хорошая подматрица для аппроксимации. Ошибка измеряется по элементам в известной части матрицы.

Если n достаточно велико, тогда метод чёрных точек возвращает приближения ко всем неизвестным элементам, за исключением 2r диагональных элементов. После этого мы должны запустить наш метод во второй раз, выбирая ведущие элементы только в полностью известных столбцах и строках. В конце концов, мы получаем скелетную аппроксимацию для матрицы R :

Выражение для диагональной части D + R аппроксимации имеет вид Как уже отмечалось выше, описанный адаптивный метод чёрных точек предназначен для задачи D + R аппроксимации. Однако, он может быть легко приспособлен для многих других проблем с другими шаблонами расположения чёрных точек.

Представленные выше алгоритмы требуют O(n2 (log n + r)) операций для построения C+R аппроксимации к неструктурированной матрице A. Эта цена складывается из двух частей: использование БПФ для вычисления FAF и из-за стратегии полного выбора ведущего элемента. В случае тёплицевых матриц это неприемлемая сложность. В дальнейшем мы покажем, как наши методы могут быть адаптированы для тёплицевого случая. В результате будет предложен алгоритм сложности O(n(log n + r2)).

1.5. Тёплицев случай дана тёплицева матрица T и A её образ Фурье:

Внедиагональные элементы матрицы A = [Akl ] можно параметризовать следующим простым и эффективным образом.

Лемма Доказательство. По определению образа Фурье A, Изменяя порядок суммирования, получаем:

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

Главная диагональ A также вычисляется за одно дискретное преобразование Фурье. Как только предварительный шаг завершён, любой элемент A может быть вычислен с помощью(1.6) очень быстро.

Вместо полного выбора ведущего элемента, мы должны использовать какую-нибудь стратегию частичного выбора. Мы предлагаем использовать ладейную схему:

1. На каждом шаге вычислить наддиагональ матрицы-невязки 2. Найти максимальный по модулю элемент в S и его индексы (i0, i0 + 1).

3. Найти максимальный по модулю элемент в i0 -ой строчке матрицы Rk и использовать его для вычисления следующего креста.

После этих модификаций суммарная сложность метода чёрных точек для тёплицевых матриц снижается до Множитель r2 появляется из-за того, что для вычисления какого-либо столбца или строчки матрицы Rk мы должны вычислить соответствующие элементы в предыдующих k 1 крестах, поэтому сложность ведёт себя как В итоге, для тёплицевого случая, C+R аппроксимация может быть вычислена быстро, если только требуемый ранг r n. Верхние оценки на ранг r, который мы называем циркулянтным рангом, будут представлены в следующем параграфе.

1.6. Существование C+R аппроксимации для некоторых классов тёплицевых матриц Как было заявлено ранее, некоторые широкие и практически важные классы символов приводят к тёплицевым матрицам, которые могут быть очень точно аппроксимированы суммой циркулянта и матрицы малого ранга.

Лемма 2 Пусть T нижнетреугольная тёплицева матрица с первым столбцом tk =, и при этом n = 1. Тогда T = C + R, где C циркулянтная матрица, а R ранг малоранговой добавки равен 2.

Доказательство. В качестве R достаточно взять одноранговую тёплицеву матрицу Напрямую проверяется, что матрица C = T R циркулянтная. Если же n = 1, то в качестве rk возьмём функцию вида:

Тогда (мы использовали, что n = 1 ). Нам нужно, чтобы ck ck n = 0.

Для этого выражение справа должно тождественно равняться нулю, т.е. = /n. Ранг же матрицы R будет равным 2, т.к. rij = i j(i Следствие 1 Если тёплицева матрица T порождена символом Доказательство. Возможны два случая: || < 1 и || > 1. В первом случае f(z) раскладывается в ряд Тейлора поэтому мы можем применить лемму 2. В случае || > 1 мы можем разложить f(z) в ряд Лорана. Единственное отличие от первого случая состоит в том, что мы получаем не нижнетреугольную, а верхнетреугольную матрицу.

Это был случай функции с простым полюсом. А что будет, если символ f будет иметь полюс более высокого порядка? Для этого нам потребуется доказать следующую лемму.

Лемма 3 Пусть T нижнетреугольная тёплицева матрица с первым столбцом tk = s(k), q натуральное число, а s(x) многочлен степени не выше q, некоторое число, Тогда T = C + R, где циркулянтная матрица, и при этом Доказательство. В качестве R возьмём тёплицеву матрицу специального вида:

где p многочлен, который ещё предстоит определить.

Нам нужно, чтобы матрица C = T R была циркулянтом. Это означает, что Отсюда получается следующее уравнение на неизвестный многочлен Возможны два случая. Пусть n = 1. Тогда полином p определяется с точностью до постоянной, и мы можем положить p(0) = 0.

т.е. значения многочлена p задаются в точках Построим теперь в точках 0, n,..., (q + 1)n многочлен (q + 1) степени c этими значениями.

Тогда многочлен p(x) p(x n) будет многочленом q -ой степени и будеть совпадать с s(x) (в силу выполнения интерполяционных условий).

Ранг матрицы R не больше q + 2. Действительно, В каждом слагаемом раскроем скобки:

подставим в выражение для rij и поменяем порядок суммирования:

Видно, что rij представимо в виде т.е. является матрицей ранга не выше q + 2.

Случай n = 1 разбирается аналогично. Единственное отличие состоит в том, что многочлен p можно выбирать не q + 1, а q -ой степени, так как добавляется одно уравнение на старший коэффициент p (который сокращался в выражении p(x) p(x n) ).

Матрица C T R будет циркулянтной в силу того, что Теорема 2 Пусть тёплицева матрица T порождена рациональным тригонометрическим символом где P, Q, L многочлены, L не имеет корней на единичной окружности, степень Q меньше степени L и они не имеют общих корней.

Тогда Доказательство. Представляем в виде суммы простых дробей:

где d степень многочлена L. Рассмотрим теперь отдельно кажc дую простую дробь вида (za). Элементы соответствующих такому символу тёплицевых матриц можно получить, разлагая (za) либо в ряд Тейлора, либо в ряд Лорана, в зависимости от того, где находится полюс a относительно единичной окружности. В одном случае это приведёт к нижнетреугольным тёплицевым матрицам, в другом к верхнетреугольным, поэтому можно ограничиться лишь изучением одного из двух случаев, второй рассматривается аналогично. Коc эффициенты ряда Тейлора (или ряда Лорана) функции (za) легко находятся. А именно, Видно, что эти коэффициенты удовлетворяют условиям леммы 3, применяя которую мы и получаем требуемый результат.

Перейдём теперь от рациональных символов к символам с логарифмическими особенностями.

Лемма 4 Пусть T нижнетреугольная тёплицева матрица с первым столбцом Тогда для любого существуют циркулянтная матрица C и матрица R ранга r такие, что причём где c0, c1, c2 зависят только от.

Доказательство. Для любого существуют fm, wm такие, что [44] Остается применить лемму 2.

Следствие 2 Пусть тёплицева матрица T порождена символом Тогда для любого существуют циркулянтная матрица C и матрица R ранга r такие, что причём Для доказательства следствия достаточно заметить, что ряд Лорана рассматриваемой функции f(z) имеет коэффициенты вида fk = k k. После этого мы можем применить лемму 3.

Аналогичным способом доказывается следующее утверждение:

Следствие 3 Пусть тёплицева матрица T порождена символом Тогда для любого существуют циркулянтная матрица C и матрица R ранга r такие, что причём Результаты данного раздела объединяет следующая Теорема 3 Пусть тёплицева матрица T порождена кусочно-аналитическим символом вида где функция g является аналитической в кольце, содержащем |z| = 1.

Тогда для любого существуют циркулянтная матрица C и матрица R такие, что причём а c0, c1, c2, c3 не зависят от n и.

Здесь, кроме слагаемых с логарифмическими особенностями, добавлена ещё и функция, аналитическая в кольце, содержащем |z| = 1. Для такой функции коэффициенты ряда Фурье убывают экспоненциально, и соответствующие тёплицевы матрицы могут быть аппроксимированы ленточными тёплицевыми матрицами. Ленточные тёплицевы матрицы могут быть приближены суммой циркулянта и матрицей малого ранга очень просто: достаточно добавить ненулевые элементы в двух противоположных углах матрицы.

Теорема 2 утверждает, что тёплицевы матрицы, порожденные произвольным рациональным тригонометрическим символом являются точной суммой циркулянта и матрицы фиксированного (не зависящего от n) ранга. Теорема 3 относится к случаю, когда символ является суммой аналитической функции и функции с логарифмическими особенностями. Соответствующие тёплицевы матрицы могут быть аппроксимированы матрицей вида C + R с очень высокой точностью. На первый взгляд может показаться, что класс таких символов достаточно узок. Однако, этот шаблон для символов содержит все примеры в литературе, посвящённые построению суперлинейных предобуславливателей. Действительно, функции вида (zk ) log(zk ) имеют скачок в производной с номером.

Для примера, рассмотрим символ f = x4, опредённый на интервале < x < и продолженный как 2 –периодическая функция на все остальные x. Он имеет скачки в первой и третьей производной.

Вычитая из f функции вида где точка, в которой происходит скачок, а A и B пропорциональны величине скачков, мы получаем аналитическую функцию, которая может быть аппроксимирована тригонометрическими полиномами (которые приводят к ленточным тёплицевым матрицам) с экспоненциально убывающей (по степени полинома) ошибкой.

Суммируя всё вышеизложенное, можно сказать, что все функции с конечным числом скачков конечного порядка могут быть аппроксимированы суммой циркулянта и матрицы малого ранга с оценкой вида (1.7).

1.7. Численные эксперименты Циркулянты, получающиеся в результате C + R аппроксимации тёплицевых матриц, естественно использовать в качестве предобуславливателей в методе сопряжённых градиентов (когда это возможно) или GMRES (во всех других случаях). Мы использовали эти циркулянты для тёплицевых матриц, порождённых следующими типичными символами (определёнными на интервале < x < и продолженными как 2 -периодические функции на все вещественные x ):

(A) Положительно определённые эрмитовы тёплицевы матрицы:

(B) Знаконеопределённые эрмитовы тёплицевы матрицы, взятые из (C) Неэрмитовы тёплицевы матрицы ( z = eix ):

В Таблице 1.3 приведены соответствующие времена (в единицах времени, равных времени вычисления одного умножения тёплицевой матрицы на вектор) для ладейной схемы. Эти матрицы аппроксимированы с относительной точностью = 107. Отметим, что на практике для различных символов могут требоваться различные точности.

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

Таблица 1.3. Время (в matvecs) для построения C + R аппроксимации ( = 107 ).

Тёплицевы матрицы, порождённые рассматриваемыми символами, удовлетворяют условиям Теорем 2 и 3, и поэтому могут быть хорошо аппроксимированы матрицами вида C + R. Для символов f10 и f оказалось, что невязка порядка машинной точности. На самом деле, это наблюдение привело нас к формулировке Теоремы 2, о которой нам не было известно до начала экспериментов. В Таблице 1.4 приведены ранги малоранговых поправок (матриц R в C + R представлении), вычисленные с помощью ладейной схемы. Для краткости будем называть эти ранги циркулянтными рангами.

Таблица 1.4. Циркулянтные ранги, ( = 107 ).

Определённый интерес представляет зависимость циркулянтных рангов от. Их типичное поведение приведено в Таблице 1.5.

Когда C+R аппроксимация построена, мы используем C1 в качестве явного предобуславливателя. Для случая симметричной и полоРанг 13 19 21 Таблица 1.5. Зависимость циркулянтного ранга от, n = 512, символ жительно определённой матрицы важно, чтобы C тоже была симметричной и положительно определённой матрицей. Численные эксперименты показывают, что симметрия нашим алгоритмом поддерживается.

Однако, наши циркулянты иногда имеют нулевые или отрицательные собственные значения, что мешает их использованию в качестве предобуславливателя. В этом случае, мы исправляем их, меняя неправильные (нулевые или отрицательные) собственные значения на 1. Это малоранговая коррекция, которая делает циркулянты положительно определёнными. В Таблице 1.6 показано число отрицательных и нулевых собственных значений циркулянтов C в C + R аппроксимации для семейства символов вида |x|k, k = 1, 2, 3, 4. Оказалось, что это число не зависит от n для других n результаты получились абсолютно такими же.

Таблица 1.6. Число отрицательных/нулевых собственных значений для построенных циркулянтов Наконец, в Таблице 1.7 показано число итераций, требуемых для решения предобусловленной системы. Порог остановки итерационного метода был установлен на 106.

1.8. Метод чёрных точек для произвольного шаблона Область применения метода чёрных точек отнюдь не ограничивается лишь диагональными матрицами. Вместо диагональной матрицы можно рассматривать произвольный шаблон разреженности S. Тогда получается следующая задача:

где S разреженная матрица, а R матрица малого ранга. Такие задачи часто встречаются в приложениях. Например, пусть A некоn 128 256 512 Таблица 1.7. Число итераций в методе сопряжённых градиентов торый большой массив данных с пропусками такая ситуация часто возникает в математической статистике при анализе рядов наблюдений. Требуется заполнить пропуски, исходя из того, что исходный массив был малого ранга. Для решения задач вида (1.8) предлагались различные методы, однако все они являются итерационными и имеют достаточно высокую вычислительную сложность. Мы же покажем, что метод чёрных точек позволяет получить приближение быстро и за конечное число операций (т.е. метод не является итерационным). Также очень интересена ситуация, когда нам неизвестен шаблон матрицы S. В этом случае не удаётся получить какие-либо теоретические результаты. Однако возможно построить эвристический алгоритм, который позволяет определить положение ненулевых элементов в матрице S (т.е. определить, какие элементы малоранговой матрицы были испорчены ). На практике это может, например, давать возможность обнаруживать ошибки экспериментаторов, которые измерили элементы матрицы A. Однако вначале подробно рассмотрим случай, когда нам известно положение ненулей в S.

Схема метода чёрных точек для случая произвольного шаблона S практически не отличается от случая диагонального шаблона. Если ранг матрицы R равен r, то мы выбираем r r подматрицу A и вычисляем скелетное разложение вида где L и U исходной матрицы. Матрицы C и R содержат пропуски ( чёрные точки ). Нетрудно понять, как будут расположены чёрные точки в матрице B. А именно, элемент (i, j) матрицы B будет известным, если в i -ой строке матрицы C нет чёрных точек и в j -ом столбце матрицы R тоже нет чёрных точек. Для выбора подматрицы опять воспользуемся методом неполной крестовой аппроксимации. На каждом шаге к шаблону разреженности добавляются новые элементы. Однако, как было объяснено выше, на каждом шаге фактически вырезаются некоторые строки и столбцы. Поэтому удобно завести два булевых массива длины n и m для того, чтобы помечать строки и столбцы, которые нельзя использовать при исключении. Полностью описание алгоритма выглядит так.

Алгоритм 1 Пусть дана матрица A, шаблон разреженности S и допустимая точность аппроксимации. Требуется построить матрицу R как можно меньшего ранга r и разреженную матрицу S c шаблоном S так, чтобы 1. Инициализация: I =, J =, r = 0.

2. Найти ведущий элемент: Если r = 0, то 4. Посчитать крест:

5. Пересчитать ранг:

6. Пересчитать чёрные списки :

7. Проверить критерий остановки, если он не выполнен, вернуться к шагу 1.

Что же получается в результате работы алгоритма? А в результате работы алгоритма получается матрица ранга (r + 1) в которой известна целая (как мы надеемся, большая!) прямоугольная подматрица размера ( card(·) число элементов в множестве).

Как же теперь найти оставшиеся чёрные точки ? Тут есть два пути. Один, более сложный в программировании, и, возможно, менее устойчивый, состоит в следующем. Мы запускаем алгоритм ещё раз, но накладываем на ведущие элементы дополнительные ограничения:

Если есть целиком известный столбец или строчка, выбирать ведущие элементы оттуда.

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

Отметим важное отличие от случая диагонального шаблона. Там, после первого шага, мы находили (n 2r) полных столбцов и строчек.

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

Другой, гораздо более простой в реализации способ состоит в следующем. Вспомним, что выполняются равенства Тогда можно посчитать V с помощью алгоритма, основанного на простой идее. Для каждого j J построим вектор и матрицу Фактически мы просто вырезали неудобные нам строки из j -ого столбца матрицы A. Теперь для нахождения j -ой строчки матрицы V необходимо решить линейную переопределённую задачу После этого нужно проделать аналогичную процедуру для нахождения матрицы U.

Какова цена этого упрощения ? Если применять обычные методы типа QR-разложения, то потребуется nr2 операций для нахождения одного столбца V. Несложно построить модификации, учитывающие происхождение матрицы U. Мы не будем обращать на это большого внимания, так как скорость вычислений метода в экспериментах оказалась достаточно высокой при n 100050000. Так как мы работаем с полной матрицей, содержащей n2 элементов, то главным ограничением становится не скорость работы алгоритма, а память, требуемая для хранения элементов матрицы. На данный момент размеры порядка десятков тысяч являются предельными. Достоинством нового метода является то, что он всегда сходится и является более устойчивым, то есть восстанавливает элементы матриц с меньшей погрешностью.

1.9. Неизвестный шаблон Рассмотрим теперь более сложный случай, когда нам неизвестен шаблон разреженности S. Что же делать? На самом деле понятно, что нам нужно. Необходимо найти такую матрицу малого ранга R, чтобы матрица A R была максимально разреженной. Однако у читателя может возникнуть естественный вопрос, а что же такое разреженность ? Как её померить? В литературе имеется большое количество вариантов количественного определения разреженности. Мы остановимся на одном из них. Пусть дан некоторый вектор x Rn.

Тогда определим разреженность S(x) по формуле где некоторый небольшой параметр. Теперь, когда у нас есть определение разреженности, мы можем сформулировать задачу:

Применим для решения задачи (1.11) метод Гаусса-Ньютона. Нетрудно показать, что на каждом итерационном шаге придётся решать минимизационную задачу вида где веса wij вычисляются по матрице с предыдущей итерации R(k) :

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

с некоторым положительными весами W = [wij ]. ( поэлементное произведение матриц). Очевидно, что задача приближения суммой разреженной матрицы и матрицы малого ранга может быть представлена в виде (1.14). Достаточно лишь положить Сделаем теперь предположение, что матрица является суммой разреженной матрицы и матрицы малого ранга:

и текущее приближение к решению достаточно близко к точному:

Тогда заменить использованную выше норму любой удобной нормой, например, фробениусовой:

Здесь опять развилка: возможно использование двух различных подходов. Первый состоит в том, что весовая матрица W фильтруется по некоторому порогу :

Здесь wmax максимальный по модулю элемент матрицы W. После этого решается задача S + R аппроксимации с известным шаблоном, задаваемым матрицей W с помощью метода чёрных точек. Задача состоит в том, чтобы не ошибиться с положением чёрных точек это контролируется параметром.

Другой подход состоит в использовании метода переменных направлений для минимизации функционала ||W (AR)||F. Вспоминим, что R = UV, U Rnr, V Rmr. Метод переменных направлений состоит в следующем.

Фиксируем U, находим V из минимизации квадратичного функционала ||W (A UV )||.

Фиксируем V, находим U из минимизации квадратичного функционала ||W (A UV )||.

Теперь определимся с тем, как мы будем пересчитывать матрицу весов W. Оказывается, что при использовании фробениусовой нормы формула (1.13) уже не является столь эффективной и требует замены.

Смысл этой формулы очевиден. Элементам с маленькими значениям невязки (A R)ij соответствуют большие веса, а элементам с большими значениями невязки маленькие. Оказывается, что функция является более эффективной.

1.10. Выводы В этой главе мы рассмотрели классическую задачу теории структрированных матриц построение циркулянтных предобуславливателей для общих и тёплицевых матриц. Эта задача была переформулирована как задача аппроксимации матрицы суммой циркулянта и матрицы малого ранга. Последняя задача была полностью решена с помощью построенного метода чёрных точек. Для тёплицева случая построен вариант метода линейной по размеру матрицы сложности, доказаны теоремы о существовании C + R аппроксимации для тёплицевых матриц широкого класса. Метод чёрных точек оказался применим не только для построения циркулянтных предобуславливателей, но и для задачи S + R аппроксимации (аппроксимации матрицы суммой разреженной матрицы плюс матрицы малого ранга). Такой метод был также построен в данной главе.

Глава 2.

Нестандартные вейвлет-преобразования 2.1. Введение В предыдущей главе мы рассмотрели задачу аппроксимации матрицы (тёплицевой или матрицы общего вида) суммой циркулянта и матрицы малого ранга. Циркулянтные матрицы диагонализуются с помощью преобразования Фурье. Однако в вычислениях оказывается полезным использовать другие быстрые преобразования для сжатия данных помимо преобразования Фурье, например так называемые вейвлет-преобразования (иногда вместо слова вейвлеты употребляют слова локальные волны, всплески ).

Классические вейвлеты и их применение в иерархическом анализе данных обычно связаны с равномерными сетками и использованием преобразования Фурье [51, 58]. Практический численный анализ приводит, как правило, к неравномерным сеткам. Построение функций и преобразований со свойствами классических вейвлетов в этом случае также возможно, но требует совершенно иной техники. В данной работе для построения нестандартных вейвлетов ( нестандартность связана с использованием неравномерных сеток) используются B-сплайны, построенные по неравномерной сетке на интервале.

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

2.2. Основные понятия и определения Напомним определение B-сплайна (см.[47]).

Определение Пусть заданы некоторые точки, среди которых по крайней мере различных:

y0.... yk+1, y0 < yk+1.

Тогда B-сплайном k-ого порядка, построенным по этим точкам, называется функция где через [y0 ;...; yk+1] обозначен оператор разделённой разности, взятой по точкам y0,..., yk+1. Функция (y x)k определяется как Носителем B-сплайна является множество supp(B) = [y0, yk+1]. Теперь рассмотрим произвольную сетку удовлетворяющую условию На этой сетке возможны совпадающие точки. Построим по этим точкам B-сплайны k-ого порядка.

Ввиду (2.1), носители этих сплайнов являются отрезками. Определим пространство V как линейную оболочку этих функций: V = Span(B1,..., Bn ). Размерность этого пространства равна n. Теперь введём более грубые В-сплайны. Для этого рассмотрим дополнительную сетку xi, i = 1,..., N + k + 1, которая содержится в сетке {xi } (т.е. для любого номера i существует номер s такой, что xi = xs ).

От этой сетки мы потребуем, чтобы её граничные точки совпадали с граничными точками исходной сетки:

Пример Если общее количество точек n+k+1 нечётно, в качестве xi можно взять все нечётные точки:

В этом случае N =. Построим по сетке {xi } неравномерные Bсплайны:

Рассмотрим пространство V = Span(B1,..., BN ). Размерность этого пространства равна N. Пространства V, V называются масштабирующими пространствами. Рассмотрим два примера масштабирующих пространств.

Пример Для любой функции f из V supp(f) n supp(Bi ) = [a, b]. Если все точки xi различны, то f непрерывна, и поэтому Пример Если мы хотим, чтобы функции из V не обращались в 0 на границе(например, в b ), нужно выбрать {xi }, i = 1,..., n + 1 различными, 2.3. Вейвлет-пространство. Масштабирующие и лифтинговые коэффициенты.

Рассмотрим функцию Bi. Её носителем является множество supp(Bi ) = [xi, xi+k+1 ]. Для простоты будем считать все эти точки различными. В силу определения дополнительной сетки cуществуют номера s0, s1 такие, что xi = xs0, xi+k+1 = xs1. В случае, когда точки xi различны, известно,что B-сплайны Bs0,..., Bs1 k1 образуют базис в пространстве (k 1) раз дифференцируемых функций, которые на каждом отрезке [xi, xi+1], i = s0,..., s1 1 являются полиномами k-ой степени. Bi принадлежит этому пространству, поэтому Коэффициенты ris называются масштабирующими коэффициентами (renement coecients). Это свойство неравномерных B-сплайнов является необходимым для того, чтобы использовать их при построении вейвлетов. Вейвлет-пространством W называется пространство, которое является дополнением (необязательно ортогональным) к пространству V в V.

(прямая сумма). Его размерность равна dim W = dim V dim V = n N. Базис в W мы будем обозначать через {i, i = 1,..., n N}.

Рассмотрим какое-нибудь пространство W с известным базисом.

Например, в примере (2.3) при k = 1 в качестве i можно выбрать функции i = B2i1, i = 1,..., (n+1)/2. Теперь построим новое вейвлетпространство с базисом Преобразование (2) является частным случаем общего преобразования, называемого лифтинговой схемой. Коэффициенты ij называются лифтинговыми коэффициентами. В данной работе эти коэффициенты выбираются так, чтобы {i } имели заданное количество нулевых моментов (p-ый момент функции f определяется как (f, xp ), где (, ) - скалярное произведение в L2(a, b) ). Подробно лифтинговая схема и различные способы выбора лифтинговых коэффициентов рассматриваются в [34].

2.4. Основная система Потребуем теперь, чтобы функции i имели m нулевых моментов.

Это означает, что где (, ) = (, )L2 (a,b) - скалярное произведение в L2 (a, b). Подставляя выражение для i, получим следующую систему на лифтинговые коэффициенты Для вычисления моментов B-сплайнов нам понадобится следующая Лемма 5 Пусть B(x) - B-сплайн k-ого порядка, построенный по произвольным (совпадающим или нет) точкам a = y0 y1... yk yk+1 = b. Тогда p-ый момент сплайна вычисляется по формуле Доказательство. Воспользуемся интегральным выражением для разделённых разностей(см. [45]): для любой (k + 1) раз дифференцируемой на [a, b] функции f Полагая в (2.7) f(x) = xk+p+1 получим, что Лемма доказана.

Используя лемму 5 и (2.5), запишем систему, которой удовлетворяют лифтинговые коэффициенты:

Через Li обозначен следующий оператор:

В этой системе i - фиксировано. Индекс j меняется от jmin до jmax, причём для того,чтобы число неизвестных совпадало с числом уравнений необходимо условие jmax = jmin + m.

2.5. Решение основной системы Займёмся теперь решением системы (2.8). Так как (2.8) должно выполняться при 0 p m, то оно равносильно тому, что для любого многочлена P(x) = m as xs+k+1. Но так как разделённая разность (k + 1) -го порядка, взятая от многочлена степени не выше k равна 0, то основная система эквивалентна уравнению(2.9). При этом P(x) – уже произвольный многочлен степени (m + k + 1).

Найдём теперь такие многочлены Pj(x), jmin j jmax, что Тогда, подставляя эти многочлены в (2.9), получим, что Для решения системы (2.8) достаточно указать какие-нибудь многочлены Pj (x). Пусть сначала все xi, i = jmin,..., jmax + k + 1 различны.

Тогда многочлен Pj однозначно определяется своими значениями в этих точках. Докажем следующую теорему.

Теорема 4 1. Многочлен Pj(x) такой, что удовлетворяет (2.10) при k = 0.

2. Многочлен Pj(x) такой, что удовлетворяет (2.10) при k 1.

Доказательство 1. То, что (2.11) удовлетворяет (2.10), проверяется непосредственной подстановкой.

2. Проверим, что (2.12) действительно даёт решение системы (2.10).

Возможны 3 случая:

а) Пусть j < r jmax. Тогда Pj (xs ) = 0 для всех s = r,..., r + k + 1. Поэтому Pj (xs ) = 0 = q(xs ) ). Поэтому [xr ;...; x(r+k+1) ]Pj(x) = [xr ;...; x(r+k+1) ]q(x) = 0, так как разделённая разность порядка k + 1 от многочлена степени не выше k равна 0.

в) Пусть r = j. Тогда Теорема доказана.

Запишем теперь многочлен Pj(x) в виде, подходящим и для случая совпадающих узлов. Pj (x) представим в виде где многочлен Pj(x) задаётся интерполяционными условиями где При k = Поэтому, если мы запишем Pj как интерполяционный многочлен Ньютона:

мы получим формулу, верную и для совпадающих узлов. Действительно, так как xj+k+1 > xj, все разделённые разности от функции f, входящие в выражение (2.13), определены и в случае совпадающих узлов. Поэтому выражение (2.13) для совпадающих узлов получается предельным переходом.

2.6. Нахождение масштабирующих коэффициентов Покажем, как находить масштабирующие коэффициенты. Напомним, что масштабирующие коэффициенты определяются как коэффициенты разложения сплайнов Bi по исходным сплайнам :

. Умножим это равенство на xp, 0 s1 s0 k 1 и проинтеp грируем от a до b. Ввиду (2.6) мы получим следующую систему на коэффициенты ris :

Эта система – система в точности того же вида, как и система (2.8).

Здесь m = s1 s0 k 1. Поэтому, используя полученные в предыдущем параграфе результаты, можно выписать формулы для масштабирующих коэффициентов.

2.7. Алгоритм вычисления вейвлет-преобразования Доказанная теорема даёт нам следующий алгоритм вычисления вейвлет-преобразования.

Алгоритм 2 Имея:

вычислить лифтинговые коэффициент при помощи следующего псевдо-кода :

Вычислить Pj (x), используя формулы интерполяции Ньютона.

Вычислить ij = s1 0 is [xi ;...; xi+k+1]Pj (x).

Для построения дискретного преобразования мы берём вектор и рассматриваем его как разложение некоторой функции f V :

где i, i = 1,..., n, представляют собой дуальный базис к базису Bi V. Обозначим через i и i дуальные базисы к Bi и i, соответственно. Тогда и требуемое преобразование выглядит как Алгоритм 3 (Один уровень вейвлет-преобразования.) Компоненты вектора ai, i = 1,..., n, Число нулевых моментов m, Порядок сплайна k, Массивы коэффициентов вычислить компоненты преобразованного вектора zi, i = 1,..., n, с помощью следующего кода :

Алгоритмы 2 и 3 реализуют один уровень вейвлет-преобразования.

Для вычисления l -уровневого преобразования нужно применить их рекурсивно l раз. Для этого необходимо лишь задать последовательность вложенных сеток.

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

Рассмотрим пока более подробно случай k = 1. Предположим для простоты, что общее число точек нечётно и xi = x2i1, is = (2i1)s, N = (n 1)/2. В этом случае, Матрица преобразования может быть представлена в виде где P матрица перестановки такая, что блочная матрица вида бидиагональная матрица размера N (n N) с элементами Матрица L блочная матрица вида где ненулевые элементы A являются лифтинговыми коэффициентами:

Обратное преобразование имеет следующий вид:

Обратное транспонированное преобразование имеет вид Отметим, что матрицы D, A, B являются ленточными. Поэтому, W, W 1, W T могут быть умножены на вектор за O(n) операций.

Нетрудно также отказаться от ограничения k = 1. В этом случае, для вычисления обратного преобразования, придётся решать линейную систему с ленточной матрицей.

В заключение определим матричное вейвлет-преобразование матрицы Z размера p p по формуле 2.8. Численные эксперименты В этом параграфе мы покажем, что нестандартные вейвлет-преобразования, адаптированные к заданной сетке превосходят преобразования Добеши. Будем сравнивать вейвлет-преобразование Добеши и нестандартное вейвлет-преобразование с парметрами k = 1 и m = 4. Для этого мы зануляем все элементы, которые меньше порога 106 и сравниваем число ненулевых элементов в каждой матрице. Вычислительная сложность нестандартных вейвлет-преобразований совпадает со сложностью вейвлет преобразований Добеши порядка 3.

Рис. 2.1. Матрица из примера 1 с p = 1. Слева вверху: исходная матрица; Справа вверху: D-3; слева снизу: NS-4 справа: D-4.

2.8.1. Пример (где p = 1 ) определённая на сетке xi = 1 cos(i/2n), вместе с её преобразованными версиями. Гладкость исходной матрицы (левый верхний угол рисунка) приводит к большому количеству маленьких элементов в преобразованной матрице. Нестандартные преобразования(левый нижний угол) приводят к гораздо большему количеству очень маленьких элементов чем и Добеши-3 (правый верхний угол) и Добеши 4 (правый нижний угол).

Рис. 2.2. Приближения к вейвлет-преобразованным матрицам из примера 1, с p = 1. Слева: новые, нестандартные преобразования с нулевыми моментами; центр: D-3; справа: D-4.

Картина становится более наглядной, если мы сравним портреты разреженности при использовании порога 106. На рисунке 2.2 показаны матрицы, аппроксимированные с помощью преобразования Добеши и с помощью нестандартных вейвлетов. Число ненулевых элементов при использовании нестандартных вейвлетов в 1, 5 раза меньше, чем при использовании преобразования Добеши порядка 4, и в два раза меньше числа ненулевых элементов в матрице, преобразованной с помощью преобразования Добеши порядка 3. Это типично для других функциональных матриц на этой сетке. В частности, мы протестировали матрицы вида (2.17) для p = 1/2, 1, 3/2, 2, 5/2. В каждом случае, нестандартные вейвлет-преобразования давали существенное уменьшение числа ненулевых элементов при заданной точности. Это продемонстрировано на рисунке 2.3. На нём показана фробениусова норма ошибки в зависимости от степени сжатия (доли ненулевых элементов в разреженной матрице). Как мы можем видеть, нестандартные вейвлеты требует на 40%. меньше памяти для хранения преобразованной матрицы.

Рис. 2.3. Фробениусова норма ошибки в зависимости от степени сжатия для нестандартных вейвлетов и для преобразования Добеши-3 для матрицы из примера 1 с p = 1/ 2.8.2. Пример определённую на сетке xi = ln(i)/ ln(n).

На рисунке 2.2 показаны результаты (с порогом 106 ) для Добешислева) и для нестандартных вейвлетов. Число ненулевых элементов при использовании нестандартных вейвлетов в два раза меньше, чем при использовании преобразований Добеши.

Эксперименты с другими функциями на этой сетке дают похожие результаты с выигрышем 50% в числе ненулевых элементов.

Даже при использовании фактора неортогональности, равного 10, (т.е. матрицы приближаются с порогом 107 ), выигрыш от использования нестандартных вейвлетов составляет от 30% до 40%.

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

Глава 3.

Тензорные аппроксимации матриц со структурированными факторами 3.1. Введение В предыдущей главе подробно рассматривались задачи нелинейной аппроксимации, возникающие при построении матричных предобуславливателей вида C + R, S + R и т.п. Перейдём теперь к абсолютно другой, на первый взгляд, задаче аппроксимации матриц задаче тензорной аппроксимации. Прежде, скажем несколько слов о тех практических задачах, к которым мы собираемся применять методы тензорной аппроксимации. Это, в первую очередь, многомерные интегральные уравнения. В диссертации речь пойдёт в основном о двумерных интегральных уравнениях, однако будет обозначено, как обобщить полученные результаты на случай большего количества измерений.

В качестве приложения рассматривается задача потенциального обтекания прямоугольного крыла (квадрат = [0, 1] [0, 1] в плоскости Oxy) установившимся потоком V(x, y, z). Для потенциала возмущённой скорости получаем краевую задачу (с условием непротекания на границе) При поиске решения в виде потенциала двойного слоя получается следующее уравнение Прандтля:

Интеграл в (3.2) понимается в смысле конечной части по Адамару.

Для численного решения уравнения (3.2) применяется метод дискретных вихрей [46],[52]. Введём сетки 0 x0... xp = 1, y0... yp = 1, и будем искать приближённое решение up (x, y) u(x, y) в виде кусочно-постоянной функции:

Выбрав точки коллокации (x0i, y0j ) (xi1, xi ) (yj1, yj), получаем систему линейных алгебраических уравнений вида где В матрично-векторной форме система (3.3) принимает вид где A = [akm ] - двухуровневая матрица (согласно терминологии [48]), в которой элемент akm находится на пересечении строки i + (j 1)p и столбца k + (m 1)p ; в векторах u = [uij] и f = [fkm ] элементы uij и fkm занимают позиции i + (i 1)p и k + (m 1)p соответственно.

Сходимость приближённого решения к единственному решению уравнения (3.2) была исследована в [53], но только для случая равномерных сеток. Было показано, что имеет место интегральная сходимость ||U Up ||L1 () 0, p, и равномерная сходимость на любом множестве точек на расстоянии не меньше от границы квадрата :

где 0 < < 1/4, а h = 1/p – шаг сетки. Оценка (3.4) указывает на то, что для достижения даже умеренной точности требуется решать линейные системы достаточно больших размеров.

Будем рассматривать два типа сеток:

(1) Равномерная сетка:

(2) Неравномерная чебышёвская сетка:

Что можно сказать о точности оценки (3.4)? Для изучения этого вопроса путем численного эксперимента необходимо уметь решать системы (3.3) для достаточно больших p. Соответствующие результаты и выводы представлены в разд. 4. Для случая неравномерных сеток результатов о сходимости вообще нет, поэтому результаты численного эксперимента особенно интересны.

В нестационарных задачах аэрогидродинамики представляет интерес решение систем вида (3.3) для большого числа различных правых частей. Именно для таких задач важно получать достаточно точные структурированные аппроксимации к A1 с малой памятью для хранения и быстрой процедурой умножения на вектор. Подобные аппроксимации могут использоваться и как явные предобусловливатели. Заметим, что в случае равномерной сетки матрица A оказывается дважды теплицевой [48] (блочно теплицевой с теплицевыми блоками).

В этом случае для A1 известны формулы Гохберга-Хайнига [12], но они содержат O(n3/2 ) параметров и поэтому малополезны. В случае неравномерных сеток вообще неизвестно какое-либо явное описание строения матриц A и A1. Поэтому поиск хороших аппроксимаций и развитие соответствующих вычислительных технологий представляются очень перспективным направлением с точки зрения практических вычислений. Аппроксимации плотных матриц, возникающих при решении интегральных уравнений [15, 16, 18, 27, 33, 40, 38] основываются на разбиении исходной матрицы на блоки и аппроксимации отдельных блоков матрицами малого ранга. Однако в случае областей, являющихся тензорными произведениями одномерных удаётся избежать использования сложных многоуровневых блочных структур с помощью тензорных аппроксимаций. Объяснить это можно следующим образом. Такие методы, как мультипольный метод Рохлина, мозаично-скелетонный метод и другие основаны на разбиении матрицы по принципу источник-приёмник. При построении тензорных аппроксимаций используется разделение по геометрическим координатам. А именно, для матрицы A мы будем искать аппроксимацию в виде где U V тензорное (кронекерово) произведение матриц, определяемое как блочная матрица вида [uijV]. Покажем связь между (3.5) и разделением переменных. Запишем приближённое равенство (3.5) в индексной форме (нумеруя, как и договаривались, элементы матрицы A 4 индексами):

Объединяя теперь пару (i1 j1 ) в один общий индекс i и пару (i1 j1) в общий индекс j получаем Индексы, объединённые в пары, соотвествуют номерам по направлению x и y соответственно для источников и приёмников. Нетрудно теперь увидеть в (3.6) задачу малоранговой аппроксимации матрицы A. Для этой задачи у нас уже есть алгоритм неполной крестовой аппроксимации, позволяющий находить малоранговую аппроксимацию за O(Nr) вычислений элементов матрицы и O(Nr2 ) дополнительных операций (здесь N = n1 n2 ).

После применения метода неполной крестовой аппроксимации мы можем сохранить матрицу в оперативной памяти. Однако этого недостаточно необходимо ещё уметь умножать матрицу на вектор. Использование одного лишь тензорного формата потребует O(N3/2 ) операций, что уже довольно существенно при N 1000. Поэтому нужно сжимать факторы! Один из возможных подходов состоит в построении C + R аппроксимации каждого фактора или структурированных представлений, основанных на так называемом малом ранге смещения (этот подход будет подробно описан в Главе 4). Но это не единственный вариант. Второй предлагаемый подход основан на использовании вейвлетов. К каждому фактору применяется вейвлет-преобразование W, после чего фактор становится псевдоразреженным. Чуть ниже мы подробнее опишем этот шаг, а пока обсудим другой вопрос. Если мы умеем быстро умножать матрицу на вектор, то можно запустить любой удобный итерационный процесс. Однако, как известно, без использования предобуславливателя обычно требуется большое количество итераций. В этой главе мы предложим несколько методов предобуславливания: тензорный предобуславливатель (тензорного ранга 1), предобуславливатель основанный на неполном LU -разложении и двухуровневый циркулянтный преобуславливатель. Сразу скажем, что при использовании неравномерных сеток циркулянтный преобуславливатель работает плохо, и нужно использовать масштабирование.

Образованная таким образом тройка тензоры + вейвлеты + преобуславливатель позволяет очень эффективно решать двумерные интегральные уравнения. Мы рассматриваем два варианта предобуславливателей: масштабированный циркулянтный предобуславливатель и предобуславливатель, основанный на использовании метода Ньютона.

Предлагаемый алгоритм состоит из следующих этапов:

(A) Приближаем A суммой тензорных произведений где Uk и Vk размеров n1 n1 and n2 n2 соответственно, и N = n1 n2 размер матрицы A. Для простоты изложения будем предполагать, что n1 = n2 = n = N1/2.

(B) Применяем вейвлет преобразования с заданным числом нулевых моментов m к каждому фактору:

Здесь W матрица вейвлет-преобразования, например преобразования Добеши или нестандартного преобразования, построенного нами. Pk и Qk псевдоразреженные матрицы. Выбирая подходящий порог = (, Pk, Qk ) и занулив все элементы в Pk и Qk меньшие по модулю, чем, мы приближаем матрицу B матрицей вида Приближение строится с заданной точностью :

(C) Построение предобуславливателя F1 для матрицы A. Есть несколько вариантов.

(D) Применяем GMRES для решения системы Возвращаем F1y как приближение к точному решению x.

Пусть теперь B = r Uk Vk. Применим вейвлет-преобразование к каждому тензорному фактору:

Затем мы выбираем порог и заменяем матрицы Pk и Qk разреженными матрицами Pk и Q соответственно, получая аппроксимацию матрицы B :

Легко проверить, что где легко вычисляется. Нестандартные преобразования неортогональны, поэтому мы должны писать где фактор неортогональности. Численные эксперименты показывают, что он порядка 10. Когда используются ортогональные преобразования, такие, как преобразования Добеши, = 1. Но W оказывается гораздо меньше в случае нестандартных преобразований, и финальная ощибка аппроксимации меньше.

Если B аппроксимация A такая, что то ошибка аппроксимации матрицы A матрицей C оценивается как Выбирая r и правильным образом, мы сможем поддерживать требуемую точность аппроксимации матрицы A. Важным результатом шага (B) является лучшее сжатие данных, но для нас это не самое главное. Основная цель этого шага уменьшение сложности матричновекторного произведения. Если обозначает число ненулей во всех Pk и Q, то C может быть умножено на вектор за O( n) операций.

3.2. Масштабированные циркулянтные предобуславливатели С помощью вышеописанных преобразований матрица A может быть умножена на вектор быстро и с необходимой точностью. Будем использовать итерационный метод типа GMRES или PCG для решения линейной системы с матрицей A. Каждая итерация выполняется быстро, но для плохообусловленной матрицы таких итераций может быть очень много, особенно если мы переходим на неравномерные сетки. Нужен хороший предобуславливатель.

В работе [11] было предложено два варианта предобуславливателей (названные ILUT и IKT). При построении обоих предобуславливателей используется представление (3.7).Для построения ILUT использовалась неполная факторизация с динамическим выбором порога в духе [28]. Для построения IKP разреженная обратная к первому (наибольшему по норме) слагаемому тензорного ранга 1 в сумме. Однако для неравномерных сеток IKP оказался неэффективным, а для ILUT потребовалось слишком много памяти.

Другая идея состоит в том, чтобы строить предобуславливатель прямо по матрице A, но используя лишь O(n) её элементов. Мы предлагаем воскресить хорошо известные конструкции многоуровневых циркулянтных предобуславливателей (см [37, 32, 41]). Однако применение именно циркулянтов эффективно лишь на равномерных сетках.

На неравномерных сетках мы предлагаем новый подход использовать масштабирование. Сначало мы масштабируем матрицу с подходящими диагональными матрицами D1 и D2. Если диагональные элементы A положительны, мы можем выбрать D1 = D2 так, чтобы все диагональные элементы A были равны 1 (другие возможности связаны с выравниванием норм столбцов и строчек A ).

Оптимальный двухуровневый циркулянтный предобуславливатель Q для A это блочный циркулянт с циркулянтными блоками, удовлетворяющий где T множество всех двухуровневых циркулянтов с такими же размерами, как и матрица A. Так как A является двухуровневой матрицей, элементы aij могут быть проиндексированы мультииндексами (i1, j1 ) и (i2, j2), где (i1, j1) определяет блок содержащий aij, а (i2, j2 ) определяет место этого элемента внутри блока. Элементы первого столбца q(i1,i2 ) оптимального циркулянта Q (он имеет длину p и его элементы можно естественным образом пронумеровать парой индексов) выражаются по формуле [37] где элементы A считаются периодическими во всех 4 индексах.

Основной трудностью при вычислении Q по формуле (3.19) является то, что получается алгоритм сложности O(n2 ), что неприемлемо. Мы предлагаем строить приближённую обратную матрицу.

В (3.19) мы вычисляем средние значения элементов, находящихся на i2 -ой периодический диагонали каждого блока вдоль i1 -ой периодической блочной диагонали. Естественно заменить среднее значение по всей последовательности на среднее значение по некоторой подвыборке с заданным шагом.

Когда Q построен, мы получаем предобуславливатель вида Двухуровневые циркулянты диагонализуются двумерным преобразованием Фурье:

где F матрица БПФ и диагональная матрица с собственными значениями. Следовательно, является явным предобуславливателем для матрицы A. Используя БПФ, мы можем умножать F на вектор за O(n log n) операций, и поэтому M1 может быть умножена на вектор O(n log(n)) операций.

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

В вычислительной алгебре известен метод уточнения обратной матрицы (см. [57]), описанный Хотеллингом [20] и Шульцем [30]. Это не что иное, как метод Ньютона с k -й итерацией вида (Xk1)(Xk Xk1) = (Xk1) для решения нелинейного уравнения Очевидно, (X + X) (X) = X1 (X + X)1 = X1 X (X + X)1, и поэтому (Xk1 )X = X1 X X1. Следовательно, Xk Xk1 = Xk1 (A X1 )Xk1 и, окончательно, k -е приближение к A1 имеет вид и сводится к двум операциям умножения матриц.

Будем применять метод (3.23) для обращения матриц малого тензорного ранга:

При перемножении матриц их тензорные ранги перемножаются.

Поэтому на каждой итерации метода Ньютона мы ищем аппроксимацию к Xk с меньшим тензорным рангом. В итоге метод оказывается эффективным в тех случаях, когда A1 допускает аппроксимацию малого тензорного ранга. Важно отметить, что в работе предлагается простая модификация метода Ньютона, намного более эффективная, чем (3.23), для таких случаев. Кроме того, для уменьшения вычислительной сложности кронекеровские сомножители аппроксимируются матрицами вида WSW, где W - вейвлет-преобразование, а S - разреженная матрица.