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