Быстрая полилинейная аппроксимация матриц и интегральные уравнения тема автореферата и диссертации по математике, 01.01.07 ВАК РФ
|
Савостьянов, Дмитрий Валериевич
АВТОР
|
||||
|
кандидата физико-математических наук
УЧЕНАЯ СТЕПЕНЬ
|
||||
|
Москва
МЕСТО ЗАЩИТЫ
|
||||
|
2006
ГОД ЗАЩИТЫ
|
|
01.01.07
КОД ВАК РФ
|
||
|
|
||||
На правах рукописи
Савостьянов Дмитрий Валериевич
Быстрая полилинейная аппроксимация матриц и интегральные уравнения
01.01.07 — Вычислительная математика
АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата физико-математических наук
Москва 2006
Работа выполнена в Институте вычислительной математики РАН Научный руководитель: член-корреспондент РАН, профессор Е. Е. Тыртышников
Официальные оппоненты:
доктор физико-математических наук, профессор А. В. Самохин, кандидат физико-математических наук И. Е. Капорин
Ведущая организация: Пензенский государственный университет
Защита состоится «_»_ 200 г. в _часов на
заседании диссертационного совета Д 002.045.01 в Институте вычислительной математики РАН по адресу: 119333, г. Москва, ул. Губкина, д. 8.
С диссертацией можно ознакомиться в библиотеке Института вычислительной математики РАН.
Автореферат разослан «_»_ 200 г.
Ученый секретарь
диссертационного совета Д 002.045.01 доктор физико-математических наук
Г. А. Бочаров
Общая характеристика работы
Актуальность работы
Диссертация посвящена построению быстрых методов приближенного умножения на плотные неструктурированные матрицы, возникающие при численном решении интегральных уравнений. Сложность заключается в том, что возникающие плотные матрицы обычно слишком большие, и их нельзя разместить в оперативной памяти компьютера или на жестком диске, а часто и вычислить за разумное время. Снижение размерности возникающих систем за счет выбора хорошей адаптивной сетки и аккуратного способа построения дискретизации является теоретически и практически сложной проблемой, решение которой требует глубокого изучения свойств исходной задачи. Зачастую в инженерных расчетах детальный анализ задачи не может быть произведен (в силу нехватки времени или информации о свойствах задачи) и дискретизация происходит на несложной, но чрезмерно мелкой сетке с использованием простых базисных функций. Для сохранения хорошей точности моделирования приходится рассматривать и решать линейные системы больших и очень больших размеров. При этом в них содержится много избыточной информации, которую можно отбросить, значительно снизив число параметров, описывающих систему, без большого ущерба для точности. Автоматическое обнаружение специфики исходных задач и получение малопараметрической аппроксимации заданных линейных систем — интересная и практически важная задача.
В работах Е. Е. Тыртышникова было показано, что с алгебраической точки зрения многие быстрые подходы к решению интегральных уравнений (мультипольные разложения, кластеризация граничных элементов, интерполяция на регулярную сетку, локальные волны) явно или неявно базируются на одном и том же наблюдении: для весьма широкого класса ядер матрицы, возникающие при дискретизации, могут быть разбиты на блоки, большинство из которых близки к матрицам малого ранга. В предложенном им же методе неполной крестовой аппроксимации для каждого блока строится билинейная аппроксимация ац ~ ciy = UL«Vja> где ранг г мал по сравнению с размерами блока. В результате для матрицы А размера N x_N строится аппроксимация А, такая что ||А — A||f < e||A||f, и хранение А и умножение на нее требует C(N loga N logb е-1) ( a > О, b > 0) ячеек памяти и арифметических действий соответственно.
При решении трехмерных объемных интегральных уравнений на тензорных сетках элементы матрицы А определяются шестью индексами Q-ij = О-м^з)ш)з и могут быть аппроксимированы в виде трилинейной суммы ai,i2i3jlj2j3 « a(i,j,)(i2j2Ki3j3) = Z«=i u(i,ji)«v(i2h)aw(i3j3)a- Матрица при этом представляется суммой прямых (кронекеровых) произведений
А « А = 21«=! иа х V« х \Уа. Как метод эффективного представления данных трилинейные аппроксимации дают сверхлинейное сжатие, их использование при решении линейных систем может приводить к алгоритмам сублинейной сложности, что значительно расширяет практические возможности применения объемных интегральных уравнений при решении трехмерных задач.
Существующие подходы к построению трилинейного разложения не позволяют оперировать большими объемами данных, характерными для интегральных уравнений. Поэтому развитие быстрого метода вычисления трилинейной аппроксимации — актуальная задача. Кроме интегральных уравнений, трилинейная аппроксимация применяется также при обработке результатов эксперимента, сигналов, статистических данных, больших массивов графической и видео-информации и в других задачах.
Цели работы
Основная цель работы состоит в развитии эффективного метода неполной крестовой аппроксимации для тензоров (трехмерных массивов). Перед автором были поставлены три задачи:
1. выяснить, возможно ли построение крестовой аппроксимации для тензоров на основе неполной информации об их элементах;
2. предложить быстрый алгоритм малопараметрической аппроксимации тензоров по небольшому числу его элементов;
3. рассмотреть применение этого алгоритма к модельным задачам, используя малопараметрическое представление не только матрицы, но и векторов итерационного процесса.
Методы исследования
Предлагаемые в данной работе алгоритмы малопараметрического описания линейных систем основаны на дискретном аналоге метода разделения переменных. Билинейная аппроксимация — это дискретный аналог расщепления источников и приемников, трилинейная аппроксимация — аналог расщепления по физическим координатам. Отметим, что при построении быстрого алгоритма трехмерной крестовой аппроксимации основное внимание уделено аппроксимации массива [ац^] в виде разложения Таккера
Т1 Т2 гз
адк ~ йцк = 91'] 'к'иц^^кк'
г'=1 ¡'=1 к'=1
с как можно меньшими значениями модовых рангов ТьТг.гз. При этом задача поиска трилинейного разложения для исходного массива [а^*] сводится к поиску трилинейного разложения для ядра [д^'к']. Поскольку размеры ядра много меньше размеров исходного массива, для решения этой задачи можно применить известные алгоритмы поиска трилинейного разложения, такие как супер-обобщенное разложение Шура, метод переменных направлений или метод Гаусса-Ньютона.
Научная новизна работы
Предлагаемый в работе алгоритм трехмерной крестовой аппроксимации (алгоритм 3) является новым и принадлежит автору. Автору не известны другие алгоритмы решения той же задачи, обладающие сравнимой или более низкой асимптотикой сложности. Теоремы, дающие обоснование возможности построения данного метода и корректности его составных частей, являются новыми и принадлежат автору, за исключением теоремы 1, полученной в соавторстве с И. В. Оселедцем и Е. Е. Тыр-тышниковым. Метод билинейной аппроксимации матрицы (алгоритм 1) является одной из возможных модификаций метода неполной крестовой аппроксимации, предложенного Е. Е. Тыртышниковым. Алгоритм поиска подматрицы наибольшего объема, то есть модуля детерминанта (алгоритм 4) предложен С. А. Горейновым. Алгоритмы построения тензорных предобусловливателей (алгоритмы 5 и 6) являются новыми и получены в соавторстве с И. В. Оселедцем. Все прочие результаты являются новыми и принадлежат автору.
Постановка задачи принадлежит Е. Е. Тыртышникову, результаты данной работы следует рассматривать как обобщение для трехмерных тензоров матричных техник, предложенных им же.
Защищаемые положения
На защиту выносятся следующие результаты работы:
1. Теорема существования трехмерной крестовой аппроксимации, которая строится по небольшому числу столбцов, строк и распорок исходного массива.
2. Алгоритм трехмерной неполной крестовой аппроксимации. В качестве входного параметра алгоритм использует процедуру вычисления любого элемента массива [ацк], однако этот массив не вычисляется и не хранится полностью. Для построения аппроксимации [а^к] в виде разложения Таккера достаточно вычислить порядка 0{пта)
(г = max(ri,г2,гз), 1 < а ^ 2) элементов массива и совершить порядка CfnT3) дополнительных действий. Алгоритм следует основным идеям, используемым при построении билинейной аппроксимации в методе неполной крестовой аппроксимации для матриц, и сохраняет главные достоинства этого метода — надежность, высокую эффективность, гибкость в применении к различным типам задач и простоту в практическом использовании.
3. Результаты применения методов полилинейной аппроксимации к модельным и практическим задачам. Тестирование алгоритма трехмерной неполной крестовой аппроксимации на модельном объемном интегральном уравнении и демонстрация сублинейной сложности возникающего метода.
Практическая значимость работы
Методы полилинейной аппроксимации могут быть использованы при численном решении интегральных уравнений, заданных на контурах, поверхностях или в конечном объеме. Они особенно эффективны в случае, когда размерность возникающей линейной системы N слишком велика, чтобы матрицу можно было за разумное время вычислить полностью и разместить в оперативной памяти или на жестком диске. Использование билинейных аппроксимаций в большинстве случаев позволяет снизить сложность вычисления матрицы и умножения на нее до C?(N log" N logb е-1) с некоторыми а > О, Ъ > 0. Использование трилинейных аппроксимаций при решении объемных интегральных уравнений на тензорных сетках позволяет снизить сложность этих операций до 0(N2/3 loga N logb ), то есть до величины асимптотически меньшей числа неизвестных.
Необходимо заметить, что методы би- и трилилинейной аппроксимации позволяют быстро вычислять аппроксимацию матрицы и умножать ее на вектор. Для их успешного применения на практике, как правило, необходимо предобусловливание, то есть приближение обратной матрицы с такими же ограничениями на память и объем вычислений. В работе предлагается способ построения супероптимального тензорного предобусловливателя малого ранга, отвечающего этим требованиям. Для решения трехмерных объемных интегральных уравнений с сублинейной асимптотикой сложности необходимо кроме того использовать в итерационном процессе структурированные векторы. Свойства сходимости возникающих итерационных процессов теоретически пока не изучены, однако в численном эксперименте с модельным уравнением демонстрируется их сходимость и сублинейная сложность.
Хотя в диссертации обсуждается только применение трилинейной аппроксимации для решения интегральных уравнений, трилинейное разложение является вычислительным ядром многих других задач, возникая при обработке данных эксперимента в физике, в задаче о ядерном магнитном резонансе в химии, при обработке сигнала в инженерной практике, слепом разделении сигналов в MIMO (multiple-input multiple-output) системах передачи информации, томографии мозга и внутренних органов в медицине, обработке статистического материала в факторном анализе, социологии и экономике, в мультимедийных приложениях для обработки больших массивов графической и видео-информации, в задачах идентификации голоса или образа, в теории сложности при поиске оптимального алгоритма. Во всех этих задачах для эффективного приближенного представления больших плотных массивов данных могут применяться предложенные в диссертации алгоритмы.
Апробация работы
Основные результаты диссертации докладывались на различных конференциях, в т.ч. на международной конференции «Матричные методы и операторные уравнения» (Москва, ИВМ РАН, 2005), на «Ломоносовских чтениях» (Москва, ВМиК МГУ, 2006), на третьей международной конференции «Математические идеи П. Л. Чебышева и их приложение к современным проблемам естествознания» (Обнинск, ИАТЭ, 2006), на 13-ой конференции ILAS (Амстердам, 2006), на «Тихоновских чтениях» (Москва, ВМиК МГУ, 2006); а также на семинаре «Вычислительные и информационные технологии в математике» (научн. рук. ак. В. В. Воеводин, проф. В. И. Лебедев, чл.-корр. Е. E. Тыртышников, ИВМ РАН).
Публикации
Основные результаты диссертации отражены в публикациях [1-7].
Структура и объем работы
Диссертация состоит из введения, четырех глав основного текста и заключения. Объем диссертации — 143 страницы. Библиография включает в себя 81 наименование. Диссертация содержит 41 рисунок и 10 таблиц.
Краткое содержание работы
Введение. Диссертация посвящена построению быстрых методов приближенного умножения на большие плотные неструктурированные матрицы, возникающие при численном решении интегральных уравнений. Называл метод «быстрым», обычно имеют в виду, что его сложность составляет o(N2), где N — размерность линейной системы; в последнее время этот термин относят к алгоритмам почти линейной по N сложности (C?(NlogaN), a > 0). Для приближенных методов отдельно выделяют зависимость сложности от параметра точности е, обычно тоже логарифмическую, указывая C?(N log"N logb е-1), a > 0, b > 0. При этом для сохранения аппроксимации параметр е выбирается порядка точности дискретизации.
Основным результатом нашей работы является развитие метода сублинейной сложности, т.е. o(N). Аналогичные алгоритмы с такой асимптотикой сложности нам не известны.
Естественно, возможность построения быстрого метода основана в первую очередь на обнаружении и использовании некоторой особенности возникающих плотных матриц. В основном обсуждается два вида специфики.
Первый вид связан с наличием у ядра трансляционной симметрии. При использовании равномерных сеток это приводит к возникновению в матрице теплицевых и ганкелевых структур. При этом затраты памяти снижаются до C(N), а сложность получаемых методов составляет ü(NlogN). Эти методы сейчас достаточно хорошо изучены и в данной диссертации рассматриваются лишь для сравнения.
Второй вид специфики связан с выделением в плотной матрице блоков, для которых может быть построено малоранговое приближение. Эта матричная трактовка была предложена Е. Е. Тыртышниковым в 1993 г. для таких известных методов построения быстрых алгоритмов, как муль-типольные разложения (В. Рохлин, L. Greengard, С. Anderson), кластеризация граничных элементов (М. Hackbusch), интерполяция на регулярную (иерархическую) сетку (Ю. М. Нечепуренко, А. Brandt), использование локальных волн (wavelets, среди авторов — R. Coifman, В. Рохлин, С. Schwab). Эти методы развивались, однако, вовсе не на основе наблюдения за неявной структурой возникающей матрицы. Первые три из них являются вообще говоря безматричными, то есть исходная матрица в явном виде не участвует в операции умножения. Они формулируются на аналитическом языке для конкретного ядра интегрального оператора (log|x—1/|х—у|, etK|x_vl/|x—у|). При этом действие интегрального оператора рассматривается в терминах взаимодействия источников и приемников. Для каждого источника (или группы источников) определяются
зоны дальнодействия, в которой потенциал взаимодействия не содержит сингулярности и интегральный оператор может быть приближен суммой нескольких операторов с разделенными переменными на основе какого-либо разложения ядра. Погрешность разложения при этом определяет точность аппроксимации. Дискретный аналог зон дальнодействия — блоки матрицы, отвечающие геометрически хорошо отделенным друг от друга элементам сеток. Они допускают малоранговое приближение, что снижает затраты на хранение и умножение на матрицу в целом. В 1995 году Е. Е. Тыртышниковым было показано, что при определенных предположениях для аппроксимации матрицы достаточно вычислить лишь небольшое число элементов каждого блока. В 1998 году для решения этой задачи им же предложен метод неполной крестовой аппроксимации — конструктивный алгоритм с почти линейной сложностью. Матрица присутствует в этом алгоритме лишь виртуально и задана процедурой вычисления любого элемента; полностью она никогда не вычисляется и не хранится.
Основной результат нашей работы основан на третьем, не описанном выше, виде специфики плотных матриц. Речь идет о представлении многоуровневой матрицы, порожденную интегральным уравнением, в виде суммы тензорных (прямых, кронекеровых) произведений. Эта операция являются дискретным аналогом расщеплению ядра по геометрическим координатам. При этом для трехуровневых матриц при определенных предположениях касательно ядра (в работах Е. Е. Тыртышникова показано, что им удовлетворяют все основные типы интегральных уравнений) достигается сверхлинейное сжатие данных, а возникающие методы решения имеют сложность, асимптотически меньшую числа неизвестных, то есть о(1Ч). Предлагаемый в нашей работе алгоритм вычисления тензорных аппроксимаций может считаться развитием (нетривиальным) метода крестовой аппроксимации для матриц и следует основным его идеям: в частности, рассматривает матрицу системы как виртуальный объект, и строит аппроксимацию на основе лишь небольшого числа элементов. За счет этого сверхлинейной скоростью обладают как этап решения системы, так и этап аппроксимации.
Первая глава диссертации посвящена крестовой аппроксимации матриц и является своеобразной предысторией метода крестовой аппроксимации для тензоров. В разделе 1.1 приводится алгоритм построения крестовой (билинейной) аппроксимации блоков матрицы на основе информации о малой части его элементов. Алгоритм основан на последовательном вычислении строк и столбцов с выбором ведущего элемента в духе метода определения ранга с помощью разложения Гаусса и с контролем точности по случайному множеству элементов. Приводится способ эффективного вычисления критерия остановки, не требующий проверки точности
аппроксимации для всех элементов блока. Определенная эмпирика, присущая этому критерию, на практике приводит к тому, что полученная аппроксимация блока имеет несколько завышенный ранг. В том же разделе рассматриваются алгоритмы снижения ранга вычисленной аппроксимации, которые мы называем методами «дожимания» или переаппроксимации. Раздел завершают численные эксперименты, демонстрирующие выгоду от применения дожимателей. В разделе 1.2 предложена параллельная версия метода. Поскольку после вычисления матричного разбиения все блоки могут быть вычислены и аппроксимированы независимо, задача состоит лишь в сбалансированном их распределении по процессорам. Предложена оценочная функция сложности и методика распределения блоков по процессорам, позволяющая достичь хорошего баланса нагрузки по процессорам при генерации матрицы и решения системы.
Вторая глава посвящена методу трехмерной крестовой аппроксимации. В разделе 2.1 содержится введение в задачу, показана выгода от применения тензорных аппроксимаций, объяснена связь тензорного разложения матриц с каноническим разложением массива, дан краткий обзор существующих методов построения трилинейного разложения, указаны их основные достоинства, ограничения и недостатки. Определено разложение Таккера, объяснена его связь с задачей трилинейной аппроксимации, приведен классический алгоритм его вычисления на основе сингулярного разложения. Поставлена задача о быстром вычислении разложения Таккера.
Содержанием раздела 2.2 является теорема о существовании разложения Таккера, в котором факторы и, V и IV состоят из элементов массива.
Теорема 1. Пусть массив Л = [ацк] может быть с точностью е приближен разложением Таккера
П Г2 тз
аук = У д^Тс'ииЛ^'-^Ис' + Е;,к, |£цк| ^ £•
1'=1 )'=1 к'=1
Тэгда существует аппроксимация того же вида, в которой факторы и, V, УУ состоят соответственно из Г] столбцов, Г2 строк и Тз распорок массива Л, а для вычисления ядра требуются лишь элементы факторов. Погрешность такой аппроксимации не превосходит
е' < ((Т1 + 1)+Т1(Г2 + 1)+Т1Г2(ТЗ + 1))Е.
Этот факт дает обоснование адаптивному методу построения разложения Таккера, последовательно наращивающему опорные базисы и,У,\У за счет вычисления новых строк, столбцов и распорок массива.
В разделе 2.3 дано описание алгоритма трехмерной крестовой аппроксимации. В простейшем виде он реализуется применением крестового метода к матрице А'3' = [а^к] (развертке массива вдоль третьего индекса)
Рис. 1. Схема работы трехмерного крестового метода
и затем рекуррентному применению его для эффективного вычисления срезок Ля = [ац] при заданных к (см. рис. 1). Затем рассматриваются более эффективные версии алгоритма, использующие единые базисы строк и столбцов для представления всех вычисленных Аь Предлагаются методы приближенной, но более быстрой реализации отдельных компонент алгоритма. В заключение раздела приведена оценка сложности полученного метода.
Раздел 2.4 посвящен деталям, которые необходимы для эффективной реализации трехмерного крестового метода. Обсуждается, как с почти линейной сложностью найти подматрицу прямоугольной матрицы, имеющую максимальный, или близкий к нему, объем. Вводится понятие доминирующей подматрицы — так называется подматрица Ап прямоугольной матрицы А, если все элементы АА51 по модулю не превосходят единицы. Обсуждается связь подматрицы максимального объема и доминирующей подматрицы. Максимальный элемент всей матрицы оценивается через максимальный элемент в ее подматрице максимального объема.
Теорема 2. Пусть В — подматрица максимального объема среди всех подматриц размера т х г в А, тогда maxel(B) > maxel(A)/(2r2 +г), где maxel(A) — максимальное значение модуля элемента матрицы А. Если при этом матрица А имеет ранг г, то maxel(B) > maxel(A)/r2.
На основе этой теоремы эффективно решается задача определения максимального (или близкого к нему) элемента в малоранговой матрице, заданной произведением факторов.
Обсуждается технология добавления векторов в опорный базис, основанная на применении реортогонализации. Показано, как при этом эффективно пересчитывать доминирующую подматрицу.
Теорема 3. Пусть в матрице II известно положение доминирующей подматрицы и на ее основе в расширенной матрице и = [1ШР] построена подматрица-надстройка размера (г+гр) х (т+тр) следующим образом: ее г строк совпадают с положением доминирующей подматрицы в и, а остальные тр — с положением доминирующей подматрицы в дополнении Шура к 11п. Тогда
1. все элементы в матрице ии^1 по модулю не превосходят тр + 1,
2. объем подматрицы Ши связан с максимально возможным среди всех подматриц такого размера неравенством
| аеЪ ии| > | det итах|/(т + гр)^+г*>/2.
Указана необходимость контроля над точностью аппроксимации на основе дополнительного, например, случайного, набора элементов. Предложена стратегия построения «нулевого приближения» к очередной вычисляемой срезке, значительно ускоряющая сходимость алгоритма после совершения нескольких первых итераций.
Теорема 4. Пусть для матрицы А размера щ X Пг существует малоранговое приближение вида
А = иВУт + Е, № = с,
где унитарные факторы и и V имеют размеры щ ХТ] и т х тг соответственно. Рассмотрим «нулевое приближение»
А0 = иВ0Ут, В0 = И51 АоУц1,
где 11а и Уа — доминирующие подматрицы в и и V, а Ап — подматрица, стоящая в А на пересечении строк, отвечающих положению и столбцов, отвечающих положению УЬ. Точность такой аппроксимации
||А- А0||р < а/П1П2Г1Г2+ 1е.
Раздел 2.5 содержит численные эксперименты по сжатию модельных числовых массивов, подобных ядрам типовых интегральных операторов. Численные эксперименты проводятся с апостериорной проверкой, на основе которой продемонстрирована надежность метода. Приведены результаты сжатия данных, показана почти линейная скорость работы.
В разделе 2.6 приведено более систематическое экспериментальное исследование зависимости асимптотики полученного алгоритма от размера
массива, от требуемой точности и от вида ядра. Показано, что практическая скорость работы метода полностью согласуется с теоретической оценкой, а в некоторых случаях оказывается даже лучше нее.
Третья глава содержит различные приложения описанных методов аппроксимации матриц для быстрого решения интегральных уравнений. Раздел 3.1 посвящен применению мозаично-скелетонного метода для решения задачи Дирихле для уравнения Гельмгольца на гладком контуре. В нем сформулированы задачи и приведены классические определения и факты из теории потенциала. Затем обсуждается метод дискретизации и способ вычисления матричных элементов. Эксперименты, проведенные с использованием параллельной вычислительной платформы, демонстрируют хорошее ускорение метода, а также значимую выгоду от применения алгоритмов переаппроксимации. Раздел 3.2 посвящен задаче гидроакустики мелкого моря, сформулированную в терминах теории потенциала. Применение для ее решения малоранговых аппроксимации позволило снизить время генерации матрицы и полного решения задачи от дней до минут. В разделе 3.3 приведены эксперименты по применению метода тензорных аппроксимаций для решения простейшего объемного интегрального уравнения. Результаты сжатия матрицы, отвечающей дискретизации на сетке т х т х т показаны в таблице 1.
Таблица 1. Сжатие матрицы модельной задачи. Точность е = Ю-5
Размер сетки т 16 32 64 128 256 512
Полная матрица 128МВ 8СВ 512СВ 32ТВ 2РВ 128РВ
Тензорная аппр. 74КВ 320КВ 1.1 МВ 6МВ 25МВ 114МВ
Ранг 11 13 15 16 17 19
Время аппр. 0.6 с. 1.5 с. 8.4 с. 54 с. 5.5 мин. 30 мин.
Для построения эффективного итерационного метода предлагается использовать структурированные векторы, аппроксимируя их в формате разложения Таккера. Предлагается способ построения тензорного предо-бусловливателя заданного ранга. Демонстрируется сходимость полученного алгоритма и высокая скорость его работы. Однако исследование концепции структурированных векторов и полученных на ее основе вариантов итерационных методов, равно как и техники предобусловливания, конечно же, является предметом отдельной и весьма обширной работы.
В четвертой главе описывается одна из особенно интересных для нас практических задач, сводимых к решению интегрального уравнения. Речь идет о распространении электромагнитной волны в неоднородной трех-
мерной среде с затуханием. Для того, чтобы избежать сложностей с правильным выбором неотражающего краевого условия, задача формулируется в виде объемного интегрального уравнения. В силу специального вида его ядра, на равномерной сетке матрица задачи является суммой трехуровневой теплицевой и трехуровневой теплиц-ганкель-теплицевой матрицы. Использование этих структур позволило сократить затраты памяти на хранение этой плотной матрицы с п6 до 0{п3) элементов, но тем не менее, оперативная память персональной станции позволяет использовать только сетки размером до 64 х 64 х 64. Для практических вычислений этого было недостаточно, так как в принципиально важном случае — приближении источника к зоне неоднородности — возникающие дискретизационные шумы не позволяли достичь требуемого разрешения. Размер сетки был увеличен за счет использования многопроцессорных систем. Предложен алгоритм распределенного хранения и умножения на теплицеву матрицу, учитывающий ее структуру и особенность кластерных систем. Его применение позволило производить расчет на сетках вплоть до 256 х 256 х 256, используя 256 процессоров. Результаты этого расчета были использованы для решения обратной задачи, то есть зондирования неоднородности, и привели к хорошему качеству восстановления искомых параметров.
Используя метод трехмерной крестовой аппроксимации, мы можем эффективно решить ту же задачу на одном процессоре, причем в том числе и на неравномерной сетке. Эксперименты демонстрируют хорошее результаты сжатия матрицы обсуждаемой задачи, что служит основанием для развития еще более эффективных методов ее решения.
Основные результаты работы
1. Получена теорема, обосновывающая возможность поиска трилинейной аппроксимации по неполной информации об исходном массиве.
2. Предложен быстрый алгоритм для нахождения аппроксимации трехмерного массива (тензора) на основе информации о небольшом числе его элементов. Основное внимание при этом было уделено разработке метода, приближающего заданный трехмерный массив в виде разложения Таккера. Алгоритм реализован набором процедур на языке фортран.
3. Алгоритм протестирован на модельных и практических задачах и продемонстрировал свою надежность, высокую эффективность, гибкость в применении к различным типам задач и простоту в практическом использовании.
Публикации по теме диссертации
1. Тыртышников Е. Е., Горейнов С. А., Савостьянов Д. В. Математическое и программное обеспечение для задач электродинамики на вычислительном кластере. Задачи электродинамики для распределенных вычислительных систем. Отчет по ФЦП «Интеграция», инв. номер 02.20.0303023, 2002г.
2. Савостьянов Д. В. Параллельная реализация метода решения объемного интегрального уравнения электродинамики на основе многоуровневых теплицевых матриц /В сборнике «Методы и технологии решения больших задач» под ред. В. И. Агошкова и Е. Е. Тыртышникова — М.: ИВМ РАН. 2004. С. 139-170.
3. Оселедец И. В., Савостьянов Д. В., Ставцев С. Л. Применение нелинейных методов аппроксимации для быстрого решения задачи распространения звука в мелком море / В сборнике «Методы и технологии решения больших задач» под ред. В. И. Агошкова и Е. Е. Тыртышникова — М.: ИВМ РАН. 2004. С. 139-170.
4. Оселедец И. В., Савостьянов Д. В. Трехмерный аналог алгоритма крестовой аппроксимации и его эффективная реализация //В сборнике «Матричные методы и технологии решения больших задач» под ред. Е. Е. Тыртышникова — М.: ИВМ РАН, 2005. С. 65-100.
5. Оселедец И. В., Савостьянов Д. В. Тензорные ранги сверхбольших трехмерных матриц /В сборнике «Матричные методы и технологии решения больших задач» под ред. Е. Е. Тыртышникова — М.: ИВМ РАН, 2005. С. 147-174.
6. Савостьянов Д. В., Тыртышников Е. Е. Применение многоуровневых матриц специального вида для решения прямых и обратных задач электродинамики Ц Вычислительные методы и программирование. 2006. Т. 7. С. 1-16.
7. Оселедец И. В., Савостьянов Д. В. Минимизационные методы аппроксимации тензоров и их сравнение Ц ЖВМиМФ. 2006. Т. 46, №10. С. 16411650.
8. Oseledets I. V., Savostianov D. V., Tyrtyshnikov E. E. Tucker dimensionality reduction of three-dimensional arrays in linear time // SIAM J. Matr. Anal. Appl., принята к публикации.
Изд. лиц. ИД №03991 от 12.02.2001. Компьютерный набор. Подписано в печать 11.11.2006. Усл. печ. л. 0.8. Тираж 100 экз.
Институт вычислительной математики РАН. 119333, Москва, ул. Губкина, 8.
Введение
1.1 Быстрые методы решения интегральных уравнений
1.2 Историческая справка
3 Основной результат данной работы.
1.4 Содержание работы по
главам
Глава 1. Мозаично-скелетонный метод
1.1 Описание и развитие мозаично-скелетонного метода
1.1.1 Описание метода.
1.1.2 Методы дожимания (переаппроксимации).
1.1.3 Численные эксперименты.
1.2 Параллельная версия метода.
1.2.1 Особенности кластерных станций.
1.2.2 Модель параллелизации и распределение нагрузки.
1.3 Выводы.
Глава 2. Метод трехмерной крестовой аппроксимации
2.1 Введение.
2.2 Теорема существования.
2.3 Трехмерный крестовый метод.
2.3.1 Как нельзя построить этот алгоритм.
2.3.2 Как можно построить этот алгоритм.
2.3.3 Как добиться почти линейной сложности.
2.3.4 Как сделать алгоритм эффективным.
2.3.5 Сложность полученного метода.
2.4 Важные детали.
2.4.1 Как найти подматрицу максимального объема.
2.4.2 Как найти наибольший элемент в срезке.
2.4.3 Как добавить векторы в ортогональный набор.
2.4.4 Как проверять точность аппроксимации.
2.4.5 Какова точность «нулевого приближения»
2.5 Модельные численные эксперименты
2.6 Асимптотика предложенного метода.
2.6.1 Теоретические оценки.
2.6.2 Практические значения ранга.
2.6.3 Время работы алгоритма.
2.7 Выводы.
1. Быстрые методы решения интегральных уравнений
В последние годы при решении сложных инженерных задач (геологоразведки, акустики моря, ЯМР в квантовой химии) все чаще используются интегральные уравнения. Интерес к ним обусловлен несколькими причинами.
• Интегральные уравнения позволяют понизить размерность задачи, сводя, например, трехмерную задачу к интегральным уравнениям записанным на границе области — двумерном многообразии.
• Интегральные уравнения позволяют сводить краевые задачи в бесконечных областях к задачам на ограниченном носителе, как например метод объемного интегрального уравнения для задачи рассеяния.
Серьезным фактором, затрудняющим применение интегральных уравнений в практических расчетах, являются затраты на вычисление, хранение и умножение на матрицу возникающей линейной системы. В подавляющем большинстве случаев она плотная, а значит при числе неизвестных тг в общем случае на ее хранение требуется тг2 ячеек памяти. Долгое время именно за счет этого интегральные уравнения проигрывали в глазах инженеров-вычислителей подходам, основанным на дифференциальных уравнениях. Последние при дискретизации приводят к разреженным матрицам, на хранение которых требуется лишь 0(п) элементов памяти. Однако возникают другие сложности: например, прямое решение такой системы приводит к заполнению, а итерационное требует построения предобусловливателя, так как возникающая матрица обычно очень плохо обусловлена.
Если мы хотим использовать методы интегральных уравнений так же эффективно, как и дифференциальные постановки, требуется предложить быстрые алгоритмы решения возникающих линейных систем. Называя метод «быстрым», обычно имеют в виду, что его сложность составляет о(п2), где тг — размерность линейной системы; в последнее время этот термин все чаще относят к алгоритмам почти линейной по п сложности, то есть сложности 0(пк^ап), с некоторым а > 0. Для приближенных методов чаще всего отдельно выделяют зависимость сложности от параметра точности е, обычно тоже логарифмическую, указывая 0[п к^ап к^ь£~1). При этом для сохранения аппроксимации параметр точности е естественно выбрать порядка точности дискретизации. При этом он оказывается связан с п степенной зависимостью (например е - тг-1 или е ~ п~2), что не нарушает почти линейную оценку сложности.
Основным результатом нашей работы является развитие метода, сложность которого асимптотически меньше, чем размерность линейной системы, т.е. о(п). Аналогичные алгоритмы с такой асимптотикой сложности нам не известны.
Естественно, возможность построения быстрого метода основана в первую очередь на обнаружении и использовании некоторой особенности возникающих плотных матриц. В основном обсуждается два вида специфики.
Первый вид связан с наличием у ядра трансляционной симметрии. При использовании равномерных сеток это приводит к возникновению в матрице теплицевых и ганкелевых структур. При этом затраты памяти снижаются до 0[п)} а сложность получаемых методов составляет (9(пк^п). В одномерном случае для решения возникающей системы можно предложить прямые методы, однако их обобщение на случай большей размерности невозможно. Это, кстати, является частой проблемой предлагаемых «быстрых» методов: они эффективны лишь для одномерных задач. Для теплицевых структур обобщение на многомерный случай, впрочем, возможно — при использовании итерационных методов решения они приводят к значительному ускорению вычислений; однако дополнительно встает вопрос об эффективном пре-добусловливании. Однако в ряде случаев использование равномерных декартовых сеток представляется невозможным (например в областях сложной формы) или невыгодным (когда выгодно сгустить сетку к некоторой особенности). Подобный пример для решения гиперсингулярного интегрального уравнения Прандтля в квадрате предложен в работе [62].
Второй вид специфики связан с выделением в большой плотной матрице блоков, для которых может быть построено малоранговое приближение. Таким образом, речь идет о приближении матрицы А другой матрицей А, такой, что:
• матрично-векторное умножение у = Ах выполняется быстро;
• объем памяти, необходимый для хранения Л, мал;
• погрешность в решении, допускаемая при замене Л на Л, сравнима с погрешностью дискретизации.
Под «быстро» и «мал» мы имеем в виду (9(nloganlogb е-1) арифметических операций и ячеек памяти. В качестве критерия допустимой погрешности мы понимаем условие ||Л —Л||р ^ £||A||f. При таком подходе для решения линейной системы должен использоваться какой-либо итерационный метод (возможно, с предобусловливанием), основной операцией которого является умножение на матрицу А. Эта операция совершается быстро за счет выявления той или иной неявной структуры в матрице, позволяющей приближенно описывать ее почти линейным числом параметров.
Именно так с матричной точки зрения можно трактовать следующие известные методы построения быстрых алгоритмов:
• мультипольные разложения [30, 31, 32, 1];
• кластеризация граничных элементов [15];
• интерполяция на регулярную (иерархическую) сетку [60, 4];
• локальные волны (wavelets) [3, 29].
Исторически, однако, эти методы развивались вовсе не на основе наблюдения за неявной структурой возникающей матрицы. Первые три из них являются вообще говоря без матричными, или операторными, то есть исходная матрица в явном виде не участвует в операции умножения. Операторные методы (обзор см., например, в [66]), как правило формулируются на аналитическом языке для конкретного ядра интегрального оператора (log |х — у|, 1/|х — у|, егк'ху'/|х —у|). При этом действие интегрального оператора рассматривается в терминах взаимодействия источников и приемников, причем при описании такого взаимодействия для каждого источника (или группы источников) определяются зоны близкодействия, куда попадает небольшое число приемников и дальнодействия, куда попадают все остальные элементы. Поскольку в зоне дальнодействия потенциал взаимодействия уже не содержит сингулярности, интегральный оператор может быть приближен суммой нескольких операторов с разделенными переменными на основе какого-либо разложения ядра. Погрешность разложения при этом определяет точность аппроксимации. С матричной точки зрения это означает, что элементы матрицы разделяются на блоки, и для блоков, отвечающих геометрически хорошо отделенным друг от друга группам элементов, строится малоранговое приближение.
Поскольку операторные методы построения аппроксимаций основаны на специальном представлении исходного оператора, для их практического применения в инженерных расчетах требуется переработка всех этапов алгоритма, начиная от процедуры дискретизации и заканчивая обработкой результата решения линейной системы. Вычислительная процедура может быть очень сложной и многоэтапной, и при практической реализации эффективного метода инженер заинтересован подвергать модификации как можно меньшее количество разработанных блоков программы. Именно этим особенно привлекательны методы, использующие непосредственно элементы блоков аппроксимируемой матрицы, например, мозаично-скелетонный метод [37, 51, 38, И, 39] или же метод иерархических матриц [16]. Для их применения надо изменить лишь процедуру построения матрицы системы и операцию матрично-векторного умножения. Подчеркнем, что способ дискретизации системы и процедура вычисления элементов матрицы остаются прежними, меняется лишь техника вычисления блока матрицы — если в исходном алгоритме он строился полностью и хранится в плотном виде, то в быстром методе для него ищется малоранговое приближение. При определенных предположениях относительно свойств гладкости ядра интегрального оператора (точный вид которого знать не требуется) для построения такого приближения достаточно вычислить лишь небольшое число элементов блока. Это означает, что матрица присутствует при работе алгоритма лишь виртуально, в виде процедуры вычисления отдельных элементов. Полностью она не вычисляется и не хранится, что позволяет снизить сложность метода до почти линейной не только на этапе решения, но и на этапе генерации матрицы.
Основной результат нашей работы основан на еще одном, не обсужденном выше, виде специфики плотных матриц. Речь идет о возможности представить многоуровневую матрицу, порожденную интегральным уравнением, в виде суммы прямых или тензорных произведений [67, 40, 17]. На операторном уровне это отвечает геометрическому расщеплению ядра по координатам. При этом для трехуровневых матриц при определенных предположениях касательно ядра (которым удовлетворяют все основные типы интегральных уравнений) достигается сверх-линейное сжатие данных, а возникающие методы решения имеют сложность, асимптотически меньшую числа неизвестных. Предлагаемый в нашей работе алгоритм вычисления тензорных аппроксимаций может считаться развитием (нетривиальным) метода крестовой аппроксимации и следует основным его идеям: в частности, рассматривает матрицу системы как виртуальный объект, заданный процедурой вычисления его элементов, и строит аппроксимацию на основе лишь небольшого числа элементов матрицы. За счет этого сверх-линейной скоростью обладают как этап решения системы, так и этап генерации аппроксиманта.
Мы считаем, что до презентации основного результата и детального его обсуждения, стоит рассказать об истории развития алгоритмов аппроксимации матриц, этапом которой, мы надеемся, является и данная работа.
2. Историческая справка
Еще в 1969 году Н. С. Бахваловым была сформулирована следующая гипотеза [47], относящаяся к корректно поставленным краевым задачам математической физики:
Пусть известно, что правая часть f(x) принадлежит некоторому компактному множеству, а пе — минимальное число вычислений правой части f(xt), достаточное для того, чтобы при любой правой части f (х) по этим значениям можно было восстановить с точностью £ решение и(х). Тогда для решения корректной задачи математической физики достаточно использовать 0{nt) значений ft и дополнительно произвести 0{пс log0 ne logb е-"1) математических действий.
Возможно, первый широко известный алгоритм с такой асимптотикой сложности был предложен в 1985 году В. Рохлиным [30] и получил название мулътиполъного метода. Весьма близок ему появившийся в 1989 году метод кластеризации граничных элементов (panel clustering), предложенный В. Хакбушем и 3. П. Новаком [15]. В этих работах элементы исходной матрицы разделялись на две группы, отвечающие близкодействию и дальнодействию. Выбор этих зон определяется видом оператора исходной задачи и особенностями метода. Поскольку сингулярность содержится только в зоне близкодействия, интегральный оператор в дальней зоне является достаточно гладким и может быть приближен оператором с разделенными переменными с помощью разложения в ряд Тейлора (как в методе кластеризации) или разложения по сферическим функциям (как в мультипольном методе). Остаточный член разложения определяет погрешность приближения, а количество членов в разложении определяет ранг аппроксимации.
Приведенные методы формулировались на геометрическом и операторном уровне, в терминах взаимодействия источников и приемников. Возможность получить малопараметрическое описание линейной системы (и тем самым построить быстрый метод ее решения) описывалась через объекты, предшествующие самой системе: оператор взаимодействия исходной задачи и геометрию расчетной области. Однако с позиций инженерной практики куда более привлекательной выглядит возможность работать на матричном уровне, то есть строить малопараметрическую аппроксимацию на основе элементов имеющейся матрицы, а не исходного оператора. В 1992 году Е. Е. Тыртышников предложил матричную трактовку существующих методов. Понятие зон близко- и дальнодействия оператора выразилось на матричном языке в виде техники иерархического разбиения матрицы на блоки. При этом блоки, отвечающие хорошо разделенным узлам сетки (зоне дальнодействия), допускали малоранговое приближение. Чтобы еще больше сблизить операторное и матричное описание, аппроксимация интегрального оператора в области дальнодействия производилась на основе интерполяционного многочлена. Если для дискретизации интегрального уравнения применяется простейший метод коллокации, то значения интегрального оператора в точках сетки совпадают с элементами матрицы, а значит их интерполяция эквивалентна аппроксимации блока матрицы. Эти результаты, полученные совместно с М. В. Мягчиловым, были опубликованы в 1992 году в работе [27]. На почве этих исследований в это же время происходило сотрудничество группы Е. Е. Тыртышникова с компанией Elegant Mathematics, результаты которого были отражены в работе [37] (1993 год). В частности, было предложено аппроксимировать блоки с помощью алгоритма типа Ланцоша (это быстрее, чем применяемое традиционно сингулярное разложение) и закреплена идея иерархической структуры блоков, идентичная концепции Н-матриц, предложенной В. Хакбушем в 1999 году в работе [16].
Для предложенного метода быстро нашлось практическое применение. Уже в 1993 году он использовался при решении серьезных промышленных задач, например для исследования рассеяния электромагнитной волны на гладкой поверхности (для фирмы Cray Research). С тех пор моделирование электромагнитных процессов было постоянной (но совсем не единственной) областью приложения быстрых алгоритмов, развиваемых в группе Е. Е. Тыртышникова. В 1995 году в работах [20, 21] рассматривались трехмерные интегральные уравнения, возникающие в задачах распространения электромагнитной волны. В связи с плотной структурой матрицы в таких задачах, ее вычисление и хранение в полном виде не представлялось возможным. Тогда для снижения арифметических затрат была использована имеющаяся блочно-теплицевая структура матрицы — это оказалось выгоднее, чем строить малопараметрические аппроксимации. Стоит отметить, что в данной работе мы предлагаем алгоритм, превосходящий по эффективности метод использования теплицевых структур. Еще более важным является то, что наш метод может быть скомбинирован с имеющейся блочно-теплицевой (или любой другой свойственной матрице) структурой, что приведет к дополнительному сжатию данных.
Начиная с 1993 года одновременно с разработкой приложений происходило и развитие основного алгоритма. В 1995 году в Докладах РАН была опубликована работа [51], в которой показывалось возможность построения малоранговой аппроксимации для блоков матрицы лишь по малой части ее элементов — нескольким столбцам и строкам. Было доказано, что такие «опорные» столбцы и строки существуют, но оставался открытым вопрос о конструктивном критерии их выбора. В том же году для дальнейшего развития метода в фонд VolkswagenStiftung был подан совместный российско-германский проект, руководителями которого стали Е. Е. Тыртышников (Россия) и С. Рязанов (ФРГ). В состав проекта с российской стороны входили также С. А. Горейнов, Н. Л. Замарашкин, И. В. Ибрагимов и М. С. Мартынов, с немецкой стороны М. Вебендорф. Идеи представления матрицы в виде блоков, допускающих малоранговую аппроксимацию, представлялись на различных конференциях по вычислительной математике и приложениям: Болгария, 1995; Enumath, Париж, 1995; Руссэ, 1996; Картона, 1996. Результатом последних выступлений стала публикация в итальянском журнале «Calcolo» [38], в это же время вышли статьи о псевдоскелетном методе в журнале Linear Algebra and Applications [11, 12]. Все это способствовало популяризации метода.
В 1998 году состоялся заключительный рабочий семинар проекта Volkswagen-Stiftung, на котором Е. Е. Тыртышников впервые продемонстрировал численные результаты аппроксимации матриц на основе неполной информации о ее элементах. Алгоритм аппроксимации стал полностью автоматическим, в нем столбцы и строки аппроксимируемого блока вычислялись адаптивно, при этом на каждом шаге вычислялись те строки (столбцы), позиции которых совпадали с расположением подматрицы наибольшего объема (модуля детерминанта) в уже вычисленных столбцах (строках). Эта техника, названная неполной крестовой аппроксимацией, была затем опубликована в работе [39].
По всей видимости, этот результат стал наиболее полезным для дальнейшего развития алгоритма. Именно тогда стало понятно, что аппроксимация блока может быть вычислена без использования всех его элементов. Техника поиска при этом могла иметь различные варианты. Например, в 1999 году М. Вебендорф предложил в своей диссертации другой алгоритм поиска опорного креста, позже включенный в статью [2]. Алгоритм Бебендорфа можно охарактеризовать как «жадный» до вычислений матричных элементов — на каждом шаге к имеющемуся кресту добавлялись только одна строка и столбец, причем решение о том, какой крест должен быть вычислен на следующем шаге, принималось только по уже имеющейся информации. Конкретнее, если на р-том шаге метода аппроксимация блока А ~ Ар1 = uavj была уточнена за счет вычисления столбца Ир и строки vp, то следующий столбец выбирается среди еще не вычисленных на той позиции, где погрешность приближения строки vp аппроксимантом Ар1 оказалась максимальной. Аналогичный выбор происходит и для строки. Точность аппроксимации оценивалась методами теории интерполяции, основанными на результатах работы [27]. При этом в оценках точности сохранялась не до конца исследованная зависимость от выбранной сетки.
Позже Е. Е. Тыртышниковым было замечено, что фактически предложенный М. Бебендорфом метод есть ни что иное, как метод разложения Гаусса с выбором ведущего элемента по столбцу и с определенным образом зафиксированной техникой выбора столбцов в активной подматрице. В алгоритмах определения ранга метода Гаусса применялся уже давно. Если ранг матрицы равен г, то после г шагов метода Гаусса ведущая подматрица оказывается в точности нулевой, что позволяет определять ее ранг. Если же матрица А приближена матрицей ранга г с некоторой точностью £, число шагов метода Гаусса зависит от техники выбора ведущих элементов, выбора первого ведущего столбца и в значительной степени определяется свойствами самой матрицы. Таким образом, по существу в работе [2] для задачи определения ранга вычисляемых блоков предложен метод Гаусса с определенной стратегией выбора ведущих элементов и эмпирически показана его применимость для некоторых типовых матриц, возникающих при решении интегральных уравнений.
Следует отметить, что выбор ведущего элемента по столбцу является стандартным способом избежать роста погрешности в методе Гаусса, но выбор новых столбцов в ведущей подматрице, вообще говоря, может быть осуществлен произвольно, лишь бы получившийся набор сохранял линейную независимость. Трактовка метода крестовой аппроксимации как варианта метода Гаусса позволила окончательно прояснить его связь с концепцией поиска подматриц наибольшего объема. В 2000 году в работах [39, 13] было доказано, что среди всех возможных аппроксимаций блока, построенных по его строкам и столбцам, наилучшей является та, строки и столбцы в которой образуют на пересечении подматрицу максимального объема. Быстрого метода поиска этой подматрицы не существует, однако можно искать близкие к ней подматрицы на основе тех или иных адаптивных стратегий. Таким образом было определено направление дальнейшей работы — эвристические техники поиска подматрицы большого объема. При этом возникала определенная вариативность метода поиска такой подматрицы, что позволяло достаточно несложным образом улучшать его вычислительные свойства. Известны примеры, в которых метод М. Бебен-дорфа «не сходится», то есть неоправданно увеличивает размер вычисляемой аппроксимации без заметного улучшения ее точности. Этот недостаток устраняется, если при выборе ведущих столбцов несколько расширить список элементов, среди которых выбирается элемент с наибольшей погрешностью аппроксимации. Например, в работе [9] в этот список кроме элементов текущей строки добавлена еще диагональ аппроксимируемой матрицы, а в работе [69] рассматривается алгоритм, содержащий дополнительную проверку точности аппроксимации на некотором случайном множестве элементов. На численных примерах была показана высокая надежность и эффективность полученного метода. В работе [69] также была предложена версия скелетного метода для параллельных кластерных платформ. Эти результаты частично включены в данную диссертацию.
Связь «принципиальных» строк и столбцов матрицы с подматрицей максимального объема активно используется и при построении трехмерного аналога метода крестовой аппроксимации для тензоров.
Необходимо также отметить, что одновременно с развитием алгоритмической части метода постоянно продолжалось расширение класса задач, для которых теоретически обосновывалось существование малоранговых аппроксимаций. Для асимптотически гладких и осцил-ляционных ядер оценки были даны С. Горейновым в [52]. Они же позже вошли в его диссертационную работу [53]. Таким образом, в 2001 году метод был в высокой степени теоретически, алгоритмически и технологически развит.
3. Основной результат данной работы
Основные результаты предлагаемой работы относятся к быстрому построению малоранговой аппроксимации для трехмерных массивов. Задача формулируется так: для заданного массива Л = [ау^] найти с заданной точностью представление в виде г
Оцк ~ X. ^^аУ^ка. а=1 с возможно меньшим значением ранга г. В применении к решению линейных систем, возникающих при дискретизации интегральных уравнений, эта задача равносильна поиску аппроксимации матрицы системы суммой прямых (кронекеровых, тензорных) произведений. Как способ представления данных тензорные аппроксимации дают сверхлинейное сжатие, что делает их крайне привлекательными при решении интегральных уравнений, когда размер возникающих массивов очень велик. Именно эта область применения алгоритма обсуждается в данной работе. Однако сама по себе трилинейная аппроксимация имеет массу других приложений:
• при обработке данных эксперимента в физике,
• в задаче о ядерном магнитном резонансе в химии,
• при обработке сигнала в инженерной практике,
• при слепом разделении сигналов в MIMO (multiple-input multiple-output) системах передачи информации,
• при томографии мозга и внутренних органов в медицине,
• при обработке статистического материала в факторном анализе, социологии и экономике,
• в мультимедийных приложениях для обработки больших массивов графической и видео-информации,
• в задачах идентификации голоса или образа,
• в теории сложности при поиске оптимального алгоритма.
Читателю, интересующимися приложениями, мы рекомендуем замечательный обзор [7].
Вообще говоря, задача аппроксимации трехмерных массивов есть частный случай задачи приближенного построения канонического (или мультилинейного) разложения для тензора произвольной размерности. Однако, сужение рассматриваемой проблемы на трехмерный случай вполне оправдано, поскольку
• именно разложение трехмерных массивов является ключевым этапом практически во всех указанных приложениях,
• алгоритм, предлагаемый для трехмерного случая, тривиально обобщается для случая большей размерности. При этом следует отметить, что сам по себе он не является тривиальным обобщением двумерного алгоритма.
Первая работа по разложению многомерных массивов принадлежит Л. Р. Таккеру, в 1966 году предложившему использовать трилинейное разложение в факторном анализе [36]. Оказалось, что разложение многомерных массивов существенным образом отличается от малорангового разложения матриц. Например, ранг трилинейного разложения (тензорный ранг) может превосходить размеры массива, причем до сих пор не получены точные верхние оценки на его значение. Тензорный ранг зависит от числового поля, в котором строится аппроксимация. Если для матриц пхп множество матриц ранга меньше п имеет меру ноль, то для тензоров это не так: среди тензоров 2x2x2 ранг 2 имеет 79% тензоров, а ранг 3 — 21 %. Эти экспериментальные результаты были получены в 1977 Крускалом. Кроме указанных результатов, в работе [23] он сформулировал и условия единственности такого разложения.
Первоначально для конструктивного построения канонического разложения использовался только метод переменных направлений [18, 5]. Этот метод весьма просто запрограммировать, стоимость одной итерации крайне невысока, и невязка при итерациях не возрастает, однако он может сойтись в некоторый локальный минимум. Достоинства и недостатки этого метода в сравнении с другими минимизационны-ми методами поиска трилинейного разложения обсуждаются нами в работе [81]. Скачок в развитии этой задачи произошел в 1997 году, когда Ливен де Латауэр предл°жил в работе [24] использовать для задач обработки сигналов технику многомерного БУБ-разложения. Это фактически разделило рассматриваемую проблему на две независимые части:
• поиск промежуточного разложения, называемого разложением Таккера
Г1 Т2 г3
Ицк ~ цц ц (*) у=1 1с'=1 который осуществляется на основе сингулярного разложения, и
• поиск трилинейного разложения для ядра 0 = [д^.
Поскольку размеры ядра 0 в практических приложениях оказываются меньше размеров исходного массива (обычно они совпадают с его тензорным рангом), трилинейное разложение для 0 строится существенно быстрее, чем для Л. Кроме того, при Г]}Г2,тз п уже разложение Таккера дает значительное сжатие данных.
Впервые задача обобщенного сингулярного разложения серии матриц — формальный аналог задачи о трилинейном разложении трехмерного массива — рассматривалась группе Е. Е. Тыртышникова в 1999 году И. В. Ибрагимовым [56], в приложении к проблеме люминесцентной спектроскопии. При этом сразу же возникло желание построить быстрый алгоритм для разложения тензоров на основе имеющегося опыта работы с матрицами. Однако тривиальные обобщения матричных подходов не приводили к успеху. Стандартные методы, основанные на ЗУБ высокого порядка (НОЗУБ) [25, 26], обладают слишком большой асимптотикой сложности, и применимы поэтому только при небольших размерах тензоров, что характерно для задач обработки сигналов или статистических приложений. Но при решении интегральных уравнений возникающие тензоры очень большие. Например для уравнения на сетке пхпхп матрица имеет п6 ненулевых элементов, и сложность НОБУБ порядка 0[п8). Даже полное вычисление такого массива и хранение его в оперативной памяти при сколь-либо существенных размерах сетки невозможно — оперативной памяти размером 2вВ оказывается достаточно лишь для матрицы, отвечающей неравномерной сетке 25 х 25 х 25. При использовании равномерной сетки этот порог поднимается до 645, но время вычисления все равно слишком велико.
Первый известный нам быстрый алгоритм трилинейного разложения был предложен в 2002 году И. В. Ибрагимовым, работавшим тогда в университете земли Саар (Германия). В его работе [19] приведен метод, который строит для матрицы, порожденной интегральным уравнением на сетке пхпхп, тензорную аппроксимацию за 0[пА) действий. Матрично-векторное умножение выполняется с той же сложностью. Алгоритм использует информацию лишь о малой части элементов матрицы, однако в нем необходимо указать некоторый набор «важных элементов» матрицы, на множестве которых происходит оптимизация. Автоматической процедуры для их определения не предлагается, что является, на наш взгляд, недостатком метода.
Другой подход к построению быстрых алгоритмов для такого рода плотных систем в это же время был предложен в работе [9]: он состоит в комбинировании тензорного разложения с эффективным представлением факторов на основе вейвлет-спарсификации. Одновременно с этим Е. Е. Тыртышников изучал вопрос существования тензорной аппроксимации и показал [67, 40], что для матриц, порожденных асимптотически гладкими функциями (а в этот класс входят дискретные аналоги всех известных интегральных уравнений) существует тензорная аппроксимация, на ранг которой можно дать весьма оптимистичные верхние оценки: г ^ с loga п.
После того, как существование аппроксимации доказано, снова встает вопрос о максимально эффективном его вычислении. По аналогии с методом крестовой аппроксимации, хотелось бы отыскать способ нахождения аппроксиманта по малой части элементов исходного массива. Здесь принципиальны два вопроса.
• Существует ли разложение (*), в котором матрицы U, V и W состояли бы из элементов исходного массива?
• Если да, то есть ли способ найти эти матрицы «быстро» ?
На оба эти вопроса в нашей работе дан положительный ответ, что составляет ее основной результат. Мы показали, что если массив Л может быть приближен разложением Таккера с погрешностью порядка е, то существует другое разложение Таккера, в котором факторы U, V, W состоят из столбцов, строк и «распорок» массива Л, и оно приближает массив Л с точностью се, где коэффициент с зависит только от модовых рангов Г], Г2, гз. Для эффективного вычисления такого разложения мы предлагаем версию крестового алгоритма, обобщенную на многомерный случай. Это обобщение не тривиально. Оно базируется на рекуррентном применения крестового метода: сначала для разложения матрицы А = [a(ij)iJ, а затем для вычисления всех требуемых срезок А^ = [a-j]. Однако для того, чтобы сделать метод по-настоящему эффективным, необходимо было предложить новый, отличный от матричного случая, формат представления данных, а также решить целый ряд задач, касающихся эффективной работы с матрицами в малоранговом представлении. Алгоритм действует адаптивно, последовательно вычисляя столбцы, строки и распорки массива Л и за счет этого наращивая размер опорных базисов U,V, W. Принципиальную роль при выборе новых вычисляемых элементов играют позиции подматриц максимального объема в текущих факторах U, V, W. Если алгоритм применяется к массиву, у которого существует представление тензорного ранга г с точностью е, то после г шагов вносимая коррекция будет иметь порядок е. На этом факте основан критерий остановки для метода.
Для массива Л размером n2 х п2 х п2, отвечающему дискретизации интегрального оператора на сетке n х n х п, сложность нашего метода составит С^п2!3) операций. При г ~ logan это на два порядка лучше, чем сложность метода, предложенного И. Ибрагимовым, и на шесть порядков лучше, чем сложность метода, основанного на непосредственном вычислении НОЗУБ. Метод включает эвристические техники выбора ведущих элементов, однако его надежность не ниже, чем у крестового метода для матриц. В серии проведенных численных экспериментов не было замечено отклонений от требуемой точности. На тестовых массивах Л размера тп х т х т, размером от тп = 103 до тп = 105, скорость метода отвечает теоретической асимптотике, а на некоторых примерах быстрее нее. Применение этого метода к матрице, порожденной объемным интегральным уравнением, позволило представить данные, которые при полном хранении требовали бы 128РВ = 128 • 250В в объеме порядка 100МВ. Вычисление потребовало всего 30 минут. Применение метода к аппроксимации матриц, порожденных объемными интегральными уравнениями, позволяет достичь сверхлинейного сжатия данных. Применяя тот же метод для сжатия векторов итерационного процесса, мы получаем метод решения сложности порядка б>(п2 к^ьп). При использовании равномерных сеток или дополнительной технологии эффективного представления факторов тензорного разложения (например, вейвлет-спарсификации), мы можем получить алгоритм почти линейной по п сложности.
4. Содержание работы по главам
Первая глава посвящена мозаично-скелетонному методу аппроксимации матриц и является своеобразной предысторией метода крестовой аппроксимации для тензоров. Ее первый раздел содержит краткое описание основных понятий и структуры мозаично-скелетонного метода. Там же приводится алгоритм построения крестовой аппроксимации блока на основе информации о малой части его элементов. Алгоритм основан на последовательном вычислении строк и столбцов с выбором ведущего элемента в духе метода определения ранга с помощью разложения Гаусса и с контролем точности по случайному множеству элементов [69, 71]. Приводится способ эффективного вычисления критерия остановки, не требующий проверки точности аппроксимации для всех элементов блока. Определенная эмпирика, присущая этому критерию, на практике приводит к тому, что полученная аппроксимация блока имеет несколько завышенный ранг. В том же разделе рассматриваются алгоритмы снижения ранга вычисленной аппроксимации, которые мы называем методами «дожимания» или переаппроксимации. Раздел завершают численные эксперименты, демонстрирующие выгоду от применения дожимателей.
В разделе 1.2 предложена параллельная версия мозаично-скелетон-ного метода. Поскольку после вычисления матричного разбиения все блоки могут быть вычислены и аппроксимированы независимо, задача состоит лишь в сбалансированном их распределении по процессорам. Предложена оценочная функция сложности и основанная на ней методика распределения блоков по процессорам, позволяющая достичь практически полной загруженности процессоров на этапе генерации матрицы. На этапе решения также достигается линейное ускорение при числе процессоров до 16.
Вторая глава посвящена методу трехмерной крестовой аппроксимации. В первом разделе содержится введение в задачу, показана выгода от применения тензорных аппроксимаций, объяснена связь тензорного разложения матриц с каноническим разложением массива, дан краткий обзор существующих методов построения трилинейного разложения, указаны их основные достоинства, ограничения и недостатки. Определено разложение Таккера, приведен алгоритм его вычисления на основе ЗУБ. Поставлена задача о вычислении разложения Таккера с почти линейной сложностью.
Содержанием раздела 2.2 является теорема о существовании разложения Таккера (*), в котором факторы 11, V и состоят из элементов массива. Это дает обоснование адаптивному методу построения разложения Таккера, последовательно наращивающему опорные базисы за счет вычисления новых строк, столбцов и распорок массива.
В разделе 2.3 дано описание алгоритма трехмерной крестовой аппроксимации. При этом изложение ведется последовательно в виде цепочки постепенно усложняемых алгоритмов. Начиная с простейшего метода, мы добавляем в него ускоряющие вычисление процедуры и рассматриваем способы быстрой реализации отдельных компонент, приходя в итоге к эффективному, но довольно объемному и сложному алгоритму. Таким образом мы демонстрируем значение всех частей метода, стараясь одновременно с этим не затуманить его основной идеи. В заключение раздела приведена оценка сложности нашего метода.
Раздел 2.4 посвящен деталям, которые необходимы для эффективной реализации трехмерного крестового метода. Обсуждается, как с почти линейной сложностью найти подматрицу прямоугольной матрицы, имеющую максимальный, или близкий к нему, объем. Использованный для этой цели алгоритм заимствован из диссертационной работы С. А. Горейнова [53]. Показано, как эффективно пересчитывать эту подматрицу при добавлении векторов в базис. Доказывается теорема, оценивающая максимальный элемент всей матрицы через максимальный элемент в ее подматрице максимального объема. На основе этой теоремы с почти линейной сложностью решается задача определения максимального (или близкого к нему) элемента в малоранговой матрице, заданной произведением факторов. Обсуждается технология добавления векторов в опорный набор, основанная на применении реортогонализации. Указана необходимость контроля над точностью аппроксимации на основе дополнительного, например, случайного, набора элементов. Предложена стратегия построения «нулевого приближения» к очередной вычисляемой срезке, значительно ускоряющая сходимость алгоритма после совершения нескольких первых итераций. Дана оценка точности полученного приближения.
Раздел 2.5 содержит численные эксперименты по сжатию модельных числовых массивов, подобных ядрам типовых интегральных операторов. Здесь мы приводим точную формулировку вычислительного метода (метод трилинейной аппроксимации является его ключевым этапом, но некоторый простор для реализации все же остается). Численные эксперименты проводятся с апостериорной проверкой, на основе которой продемонстрирована надежность метода. Приведены результаты сжатия данных, показана почти линейная скорость работы.
В разделе 2.6 приведено более систематическое экспериментальное исследование зависимости асимптотики полученного алгоритма от размера массива, от требуемой точности и от вида ядра. Показано, что практическая скорость работы метода полностью согласуется с теоретической оценкой, а в некоторых случаях оказывается даже лучше нее.
Третья глава содержит различные приложения описанных методов аппроксимации матриц для быстрого решения интегральных уравнений. Раздел 3.1 посвящен применению мозаично-скелетонного метода для решения задачи Дирихле для уравнения Гельмгольца на гладком контуре. В нем сформулированы задачи и приведены классические определения и факты из теории потенциала. Затем обсуждается метод дискретизации и способ вычисления матричных элементов. Эксперименты, проведенные с использованием параллельной вычислительной платформы, демонстрируют хорошее ускорение метода, а также значимую выгоду от применения алгоритмов переаппроксимации. Раздел 3.2 посвящен задаче гидроакустики мелкого моря, сформулированную в терминах теории потенциала. Применение для ее решения малоранговых аппроксимации позволило снизить время генерации матрицы и полного решения задачи от дней до минут. В разделе 3.3 приведены эксперименты по применению метода тензорных аппроксимаций для решения простейшего объемного интегрального уравнения. Для построения эффективного итерационного метода предлагается использовать структурированные векторы, аппроксимируя их в формате разложения Таккера. Предлагается способ построения тензорного пре-добусловливателя заданного ранга. Демонстрируется сходимость полученного алгоритма и высокая скорость его работы. Однако исследование концепции структурированных векторов и полученных на ее основе вариантов итерационных методов, равно как и техники предо-бусловливания, конечно же, является предметом отдельной и весьма обширной работы.
В четвертой главе описывается одна из особенно интересных для нас практических задач, сводимых к решению интегрального уравнения. Речь идет о распространении электромагнитной волны в неоднородной трехмерной среде с затуханием. Для того, чтобы избежать сложностей с правильным выбором неотражающего краевого условия, задача формируется в виде объемного интегрального уравнения. В силу специального вида его ядра, на равномерной сетке матрица задачи является суммой трехуровневой теплицевой и трехуровневой теплиц-ганкель-теплицевой матрицы. Использование этих структур позволило сократить затраты памяти на хранение этой плотной матрицы с п6 до (9(п3) элементов, но тем не менее, оперативная память персональной станции позволяет использовать только сетки размером до 64 х 64 х 64. Для практических вычислений этого было недостаточно, так как в принципиально важном случае — приближении источника к зоне неоднородности — возникающие дискретизационные шумы не позволяли достичь требуемого разрешения. Размер сетки был увеличен за счет использования многопроцессорных систем. Предложен алгоритм распределенного хранения и умножения на теплицеву матрицу, учитывающий ее структуру и особенность кластерных систем. Его применение позволило производить расчет на сетках вплоть до 256 х 256 х 256, используя 256 процессоров. Результаты этого расчета были использованы для решения обратной задачи, то есть зондирования неоднородности, и привели к хорошему качеству восстановления искомых параметров.
Используя метод трилинейной крестовой аппроксимации мы можем эффективно решить ту же задачу на одном процессоре, причем в том числе и на неравномерной сетке. Предварительные эксперименты показывают, что для матрицы обсуждаемой задачи эффективно строится малоранговое тензорное приближение. Это дает основание для еще более эффективных методов решения электродинамической задачи.
4.8. Выводы
В этой главе описана одна из особенно интересных для нас практических задач, сводимых к решению интегрального уравнения — задача о распространении электромагнитной волны в неоднородной трехмерной среде с затуханием. В силу специального вида ядра возникающего интегрального уравнения на равномерной сетке матрица задачи является суммой трехуровневой теплицевой и трехуровневой теплиц-ганкель-теплицевой матрицы. Использование этих структур позволило нам сократить затраты памяти на хранение этой плотной матрицы с п6 до 0{пъ) элементов, но тем не менее, оперативная память персональной станции позволяет использовать только сетки размером до 64 х 64 х 64. Для практических вычислений этого было недостаточно, так как в принципиально важном случае — приближении источника к зоне неоднородности — возникающие дискретизационные шумы не позволяли достичь требуемого разрешения. Размер сетки был увеличен за счет использования многопроцессорных систем. Предложен алгоритм распределенного хранения и умножения на теплицеву мат
Заключение
В заключение диссертации сформулируем ее основные результаты.
1. Доказана теорема существования трехмерной крестовой аппроксимации, которая строится по небольшому числу столбцов, строк и распорок исходного массива.
2. Предложен алгоритм трехмерной неполной крестовой аппроксимации. В качестве входного параметра алгоритм использует процедуру вычисления любого элемента массива [а^, однако этот массив не вычисляется и не хранится полностью. Для построения аппроксимации [а^ в виде разложения Таккера с модо-выми рангами (г1,г2,гз) достаточно вычислить порядка 0{пта) (г = тах(г1,г2,гз), 1 ^ а ^ 2) элементов массива и совершить порядка (^(пг3) дополнительных действий.
3. Проведено тестирование алгоритма трехмерной неполной крестовой аппроксимации на модельном объемном интегральном уравнении и показано, что на основе этого метода можно построить алгоритмы сублинейной по числу неизвестных сложности.
Кроме указанных результатов в диссертации предложен метод дожимания для блоков, возникающих в ходе работы мозаично-скелетон-ного метода, предложена параллельная версия мозаично-скелетонного алгоритма и проведено ее тестирование на практических примерах, продемонстрировано применение мозаично-скелетонного метода к решению задачи Дирихле для двумерного уравнения Гельмгольца и к задаче гидроакустики. Рассмотрена задача распространения электромагнитной волны в трехмерной среде с затуханием, сформулированная в виде интегрального уравнения. Предложен параллельный алгоритм решения этой задачи, использующий специфику возникающих структурированных матриц. Алгоритм применен для решения прямой и обратной задач электродинамики. Показано, что применение трехмерного крестового метода может приводить к значительному сжатию матрицы системы и в перспективе позволит решать ту же задачу на неравномерных сетках без использования многопроцессорных систем.
1. Bebendorf M. Approximation of boundary element matrices // Numer. Math. 2000. V. 86, No. 4. P. 565-589.
2. Beylkin G., Coifman R., Rokhlin V. Fast Wavelet Transforms and Numerical Algorithm I // Yale University Preprint. 1990.
3. Brandt A., Lubrecht A. Multilevel matrix multiplication and fast solution of integral equations. // J. Comp. Phys. 1990. V. 90. P. 348-370.
4. Caroll J. D., Chang J. J. Analysis of individual differences in multidimensional scaling via n-way generalization of Eckart-Young decomposition // Psychometrica. 1970. V.35 p. 283-319.
5. Chan R. H., Tyrtyshnikov E. E. Spectral Equivalence and Proper Clusters for Boundary Element Method Matrices // Internat. J. Numer. Methods Engrg. 2000. V. 49. No. 9. P. 1211-1224.
6. Comon P. Tensor decomposition: State in the Art and Applications / In IMA Conf. mathematics in Signal Processing, Warwick, UK, Dec. 18-20,2000.
7. Dennis J. E., Schabel R. B. Numerical methods for unconstrained optimization and nonlinear equations — Plentice-Hall: Englewood Clis, 1983.
8. Ford J. M., Tyrtyshnikov E. E. Combining Kronecker product approximation with discrete wavelet transforms to solve dense, function-related systems //SIAM J. Sci. Comp. 2003. V. 25, No. 3. P. 961-981
9. Gao G., Fang S., Torres-Verdin C. A new approximation for 3D electromagnetic scattering in the presence of anisotropic conductive media / 3DEMIII Workshop, 2003, Adelaide.
10. Goreinov S. A., Tyrtyshnikov E. E., Yeremin A. Y. Matrix-free iteration solution strategies for large dense linear systems // Numer. Linear Algebra Appl. 1996. V. 4 (5). P. 1-22.
11. Goreinov S. A., Tyrtyshnikov E. E., Zamarashkin N. L. A theory of Pseudo-Skeleton Approximations // Linear Alebra Appl. 1997. V. 261. P. 1-21.
12. Goreinov S. A., Tyrtyshnikov B. B. The maximal-volume concept in approximation by low-rank matrices / Contemporary Mathematics. 2001. V. 208. P. 47-51.
13. Greengard L., Rokhlin V. The rapid evaluation of potential fields in three dimentions / Lecture Notes in Mathematics. 1988. V. 1360. P. 121-141.
14. Hackbusch W., Nowak Z. P. On the fast matrix multiplication in the boundary elements method by panel clustering / Numer. math. 1989. V. 54 (4). P. 463-491.
15. Hackbusch W. A sparce matrix arithmetic based on H-matrices. Part I: Introduction to H-matrices / Computing. 1999. V. 62. P. 89-108.
16. Hackbush W., Khoromskij B. N., Tyrtyshnikov B. B. Hierarchical Kronecker tensor-product approximations. / J. Numer. Math. 2005. V. 13. P. 119-156.
17. Kress R. Linear Integral Equations. / Appl. Math. Sci., 1982.
18. Kruskal J. B. Three-way arrays: Rank and uniqueness for 3-way and n-way arrays / Linear Algebra Appl. 1977. V. 18. P. 95-138.
19. De Lathauwer L. Signal processing based on multilinear algebra. Ph.D. Thesis. — Katholieke Univ.: Leuven, 1997.
20. De Lathauwer L., de Moor B., Vandewalle J. A multilinear singular value decomposition // SIAM J. Matrix Anal. Appl 2000. V.21. P. 1324-1342.
21. De Lathauwer L., de Moor B., Vandewalle J. On best rank-1 and rank-(Rl,R2,.,RN) approximation of high-order tensors //SIAM J. Matrix Anal Appl 2000. V. 21. P. 1324-1342.
22. Myagchilov M. V., Tyrtyshnikov B. B. A fast matrix-vector multiplier in discrete vortex method // Russian Journal of Numer. Anal and Math. Modelling. 1992. V. 7 (4). P. 325-342.
23. Olshevsky V., Oseledets I. V., Tyrtyshnikov B.B. Tensor properties of multilevel Toeplitz and related matrices. // Lin. Algebra Appl. 2006. V. 1. P. 1-21.
24. Petersdorff T., Schwab C., Schneider R. Multiwavelets for second kind integral equations. // Preprint of University of Maryland. 1991.
25. Rokhlin V. Rapid solution of integral equations of classical potential theory // J. Comput. Physics. 1985. V. 60. P. 187-207.
26. Rokhlin V. A fast algorithm for particle simulations. // J. Comput. Physics. 1987. V. 73. P. 325-348.
27. Rokhlin V. Rapid solution of integral equations of scattering theory in two dimensions // J. Comput. Physics. 1990. V. 86. P. 414-439.
28. Saad Y., Schultz M. H. GMRBS: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM F. Scientific and Stat. Comp. 7:856-869, 1986.
29. Smirnov Yu. G., Tsupack A. A. Volume singular integral equations for solving diffraction problem of electromagnetic waves in mocrowave oven. // Proc. of European Symp. on Numer. Meth. in Electromagnetics March 6-8, 2002, Tolouse, Prance. P. 172-176.
30. Sun X., Pitsianis N. P. A matrix version of the fast multipole methos //SIAM Review. 2001. V. 43, No. 2. P. 289-300.
31. Tucker L. R. Some mathematical notes on three-mode factor analysis. //Psychometrika. 1966. V. 31. P. 279-311.
32. Tyrtyshnikov B. B. Matrix approximations and cost-effective matrix-vector multiplication — M.: HBM PAH, 1993.
33. Tyrtyshnikov E. E. Mosaic-skeleton approximations Ц Calcolo. 1996. V. 33 (1-2). P. 47-57.
34. Tyrtyshnikov E. E. Incomplete cross approximation in the mosaic-skeleton method // Computing. 2000. V. 4. P. 367-380.
35. Tyrtyshnikov E. E. Kronecker-product approximations for some function-related matrices Ц Linear Algebra Appl. 2004. V. 379. P. 423-437.
36. Tyrstyshnikov E. E. Fast computation of Toeplitz forms and some multidimentional integrals. Rus. J. Numer. Anal 2005. V. 20, No. 4. P. 383-380
37. Tyrtyshnikov E. E. Optimal and Super-optimal Circulant Preconditioned // SI AM J. Matrix Anal. Appl, 1992. V. 13. №2, p. 459-473.
38. Watson G. N. A Threatise on the Theory of Bessel Functions. — Cambridge University Press, 1996.
39. Zhang Т., Golub G.H. Rank-one approximation to high-order tensors //SIAM J. Matrix Anal Appl 2001. V. 23. P. 534-550.
40. Zwamborn A. P. M., Van der Berg. The three-dimensional weak form of the conjugate gradient FFT method for solving scattering problems. IEEE Trans. Microwave Theory Tech., 1992, MTT-40, 9:1757-1765.
41. Ахмед H., Pao К. Ортогональные преобразования при обработке цифровых сигналов. М.: Связь, 1980.
42. Бахвалов Н. С. Об оптимальных методах решения задач // Appl. Mat. 1969. V. 13. P. 27-43.
43. Веккенбах Э., Веллман Р. Неравенства. М.: Мир, 1965.
44. Воеводин В. В., Тыртышников Е. Е., Вычислительные процессы с теплицевыми матрицами. — М.: Наука, 1987.
45. Голуб Дж., Ван Лоун Ч. Матричные вычисления / Пер. с. англ. М.: Мир, 1999.
46. Горейнов С. А., Замарашкин Н. Л., Тыртышников Е. Е. Псевдоскелетные аппроксимации матриц //Доклады Российской академии наук. 1995. V. 343 (2). Р. 151-152.
47. Горейнов С. А. Мозаично-скелетонные аппроксимации матриц, порожденных асимптотически гладкими и осцилляционными ядрами //Матричные методы и вычисления — М.: ИВМ РАН, 1999. С. 42-76.
48. Горейнов С. А. Псевдоскелетные аппроксимации для блочных матриц, порожденных асимптотически гладкими ядрами. Диссертация на соискание ученой степени кандидата физико-математических наук. — М.: ИВМ РАН, 2001.
49. Довгий С.А., Лифанов И.К. Методы решения интегральных уравнений. Киев: Наукова думка, 2002. - 344 с.
50. Еремин Ю. А., Ивахненко В. И. Строгие и приближенные модели царапины на основе метода интегральных уравнений //Дифф. уравнения. Т. 37, №10. 2001. С. 1386-1394.
51. Ибрагимов И. В. Новый подход к решению задачи обобщённого сингулярного разложения // Матричные методы и вычисления М.: ИВМ РАН, 1999. С. 193-201.
52. Колтон Д., Кресс Р. Методы интегральных уравнений в теории рассеивания. М.: Мир, 1987. - 311с.
53. Лифанов И.К. Метод сингулярных интегральных уравнений и численный эксперимент. М.: ТОО Янус, 1995 520с.
54. Мартынов М. С. Результаты применения мозаично-скелетонного метода для решения интегральных уравнений на параллельном компьютере МВС-100 Ц Численный анализ и математическое моделирование М.: ИВМ РАН, 1999. С. 109-115.
55. Нечепуренко Ю. М. Быстрые численно устойчивые алгоритмы для широкого класса линейных дискретных преобразований. — М.: 1985. Препринт ОВМ АН СССР №98.
56. Никифоров А. Ф., Уваров В. Б. Специальные функции математической физики — М.: Наука, 1984.
57. Оселедец И. В., Тыртышников Е. Е. Приближенное обращение матриц при решении гиперсингулярного интегрального уравнения // ЖВМиМФ, 2005. Т. 45, №2. С. 315-326.
58. Парлетт Б. Симметричная проблема собственных значений. Численные методы / Пер. с англ. — М.: Мир, 1983.
59. Самохин А. Б. Интегральные уравнения и итерационные методы в электромагнитном рассеянии — М.: Радио и связь, 1998.
60. Самохин А. Б. Исследование задач дифракции электромагнитных волн в локально-неоднородных средах Ц ЖВМиМФ. 1990. Т. 30. С. 107-121.
61. Тыртышников Е. Е. Методы быстрого умножения и решение уравнений Ц Матричные методы и вычисления — М.: ИВМ РАН, 1999. С. 4-41.
62. Тыртышников Е. Е. Тензорные аппроксимации матриц, порожденных асимптотически гладкими функциями / Матем. сб. 2003. Т. 194, 6. С. 147-160.
63. Публикации по теме диссертации:
64. Oseledets I. V., Savostianov D. V., Tyrtyshnikov E. E. Tucker dimensionality reduction of three-dimensional arrays in linear time Ц SIAM J. Matr. Anal. Appl., принято к публикации.
65. Савостьянов Д. В. Мозаично-скелетонные аппроксимации. Выпускная квалификационная работа на степень бакалавра — М.: ИВМ РАН, 2001.
66. Тыртышников Е. Е., Горейнов С. А., Чугунов В. Н., Святский Д. А., Савостьянов Д. В. Математическое обеспечение для решения задач на многопроцессорных вычислительных кластерах. Отчет по ФЦП «Интеграция», инв. номер 02.200.203461, 2001г.
67. Савостьянов Д. В. Параллельная реализация метода решения объемного интегрального уравнения электродинамики на основе многоуровневых теплицевых матриц Ц Методы и технологии решения больших задач — ИВМ РАН. 2004. С. 139-170.
68. Оселедец И. В., Савостьянов Д. В., Ставцев С. Л. Применение нелинейных методов аппроксимации для быстрого решения задачи распространения звука в мелком море Ц Методы и технологии решения больших задач — ИВМ РАН. 2004. С. 139-170.
69. Савостьянов Д. В., Тыртышников Е. Е. Применение многоуровневых матриц специального вида для решения прямых и обратных задач электродинамики / Вычислительные методы и программирование. 2006. Т. 7. С. 1-16.
70. Оселедец И. В., Савостьянов Д. В. Методы разложения тензора / Матричные методы и технологии решения больших задач- М.: ИВМ РАН, 2005. С. 51-64.
71. Оселедец И. В., Савостьянов Д. В. Трехмерный аналог алгоритма крестовой аппроксимации и его эффективная реализация / Матричные методы и технологии решения больших задач- М.: ИВМ РАН, 2005. С. 65-100.
72. Оселедец И. В., Савостьянов Д. В. Быстрый алгоритм для одновременного приведения матриц к треугольному виду и аппроксимации тензоров Ц Матричные методы и технологии решения больших задач М.: ИВМ РАН, 2005. С. 101-116.
73. Оселедец И. В., Савостьянов Д. В. Об одном алгоритме построения трилинейного разложения Ц Матричные методы и технологии решения больших задач — М.: ИВМ РАН, 2005. С. 117130.
74. Оселедец И. В., Савостьянов Д. В. Минимизационные методы аппроксимации тензоров и их сравнение Ц Матричные методы и технологии решения больших задач — М.: ИВМ РАН, 2005. С. 131-146.
75. Оселедец И. В., Савостьянов Д. В. Тензорные ранги сверхбольших трехмерных матриц / Матричные методы и технологии решения больших задач — М.: ИВМ РАН, 2005. С. 147-174.
76. Оселедец И. В., Савостьянов Д. В. Минимизационные методы аппроксимации тензоров и их сравнение Ц ЖВМиМФ. 2006. Т. 46, МО. С. 1641-1650.
