ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2009, том 49, № 10, с. 1741-1756
УДК 519.61
ПРИБЛИЖЕННОЕ УМНОЖЕНИЕ ТЕНЗОРНЫХ МАТРИЦ НА ОСНОВЕ ИНДИВИДУАЛЬНОЙ ФИЛЬТРАЦИИ ФАКТОРОВ1)
© 2009 г. Д. В. Савостьянов, Е. Е. Тыртышников
(119333 Москва, Губкина 8, Институт вычислительной математики РАН) e-mail: dmitry.savostyanov@gmail.com; tee @inm.ras.ru Поступила в редакцию 10.03.2009 г.
Предлагаются алгоритмы приближенного вычисления произведения матриц C « C = A ■ B, где матрицы A и B заданы тензорным разложением в каноническом формате или в формате Таккера ранга r. Матрица C в виде полного массива не вычисляется. Вместо этого она представляется сначала аналогичным разложением с избыточным значением ранга, а затем переаппроксимируется (сжимается) с целью уменьшения ранга в рамках заданной точности. Известные алгоритмы переаппроксимации в данном случае требуют хранения массива из r2d элементов, где d — размерность пространства. Из-за ограничений по памяти и быстродействию они неприменимы уже для типичных значений d = 3 и r ~ 30. В данной работе предлагаются методы, основанные на аппроксимации модовых факторов для C по индивидуально выбранным критериям точности. В качестве приложения рассматривается вычисление трехмерного потенциала Кулона. Показано, что предложенные методы эффективны, когда значение r достигает нескольких сотен, а сложность операций по переаппроксимации (сжатию) C невелика по сравнению с предварительным вычислением факторов тензорного разложения C с избыточным значением ранга. Библ. 38. Табл. 4.
Ключевые слова: многомерные массивы, многомерные операторы, малопараметрические представления, каноническое разложение, разложение Таккера, скелетонная аппроксимация, малоранговые матрицы, сжатие данных, быстрая рекомпрессия, потенциал Кулона.
1. ВВЕДЕНИЕ
Поскольку d-мерный массив размера n по каждому направлению содержит nd элементов, для эффективного решения d-мерных задач необходимы специальные форматы представления данных и вычислительные схемы, основанные на данных форматах. Уход от экспоненциальной зависимости от d является серьезной открытой проблемой. Тем не менее имеется ряд работ, в которых предлагаются форматы данных, где число параметров зависит от n линейно или почти линейно (см. [1]—[9]).
В настоящей статье решаем следующую задачу. Для матриц A и B, заданных некоторым тензорным разложением, требуется найти приближенное тензорное разложение их произведения C ~ C = A ■ B, число определяющих параметров которого меньше, чем у точного разложения C. Это дает основу для построения приближенных пространств Крылова или приближенного обращения матрицы методом Ньютона (см. [1], [6], [8], [10], [11]).
В [12] мы рассмотрели эту задачу для основных тензорных форматов и предложили метод решения, обладающий сублинейной сложностью по полному числу данных. Однако при этом выполняется глобальная фильтрация массива, размер которого экспоненциально растет по d. В настоящей работе предлагаются методы, основанные на предварительной индивидуальной фильтрации модовых массивов, размер которых не зависит от d. Существенно то, что малая погрешность в отдельной моде может сильно изменить элементы всего тензора. Поэтому независимая фильтрация модовых факторов требует аккуратного выбора индивидуальных критериев точности аппроксимации.
Общая схема умножения с независимой фильтрацией описана в разд. 2. В разд. 3 приведен метод расчета индивидуальных критериев, гарантирующих требуемую точность ответа. В разд. 4 обсуждаются различные методы малоранговой аппроксимации. Там же предложен быстрый алго-
1) Работа выполнена при финансовой поддержке РФФИ (коды проектов 08-01-00115, 09-01-12058), гранта РФФИ/DFG 09-01-91332 и Программой приоритетных исследований Отделения матем. наук РАН.
1741
ритм сжатия на основе неполных крестовых аппроксимаций. В разд. 5 разработанные методы прилагаются к вычислению трехмерного потенциала Кулона. В Приложении приводятся алгоритмы вычисления индивидуальных критериев в трехмерном случае, выраженные через вызовы оптимизированных процедур библиотеки BLAS/LAPACK.
Матрицы и векторы рассматриваются как многомерные массивы (тензоры). Под тензором размерности d мы понимаем массив В = [Ьуу ^ ] с d индексами. Отдельные индексы и соответствующие им направления называются модами, а число значений индексов — модовыми размерами. В случае модовых размеров п1, ..., п11 говорят, что тензор имеет размер п1 х ... х пл.
Тензор может возникать, например, при дискретизации функции на тензорной сетке. Оператор, действующий на подобную функцию, будет иметь дискретный аналог в виде матрицы А = = [ а( ¡1... )(/! ■■■]й) ]. Для простоты будем считать, что все модовые размеры равны п, а также откажемся от двухуровневой индексации и будем писать В = [Ъу. к] и А = [«(¡у...к-)(у...к)] (предлагаемые методы выглядят единообразно в случае d = 3 и в пространствах большей размерности).
Для компактного представления тензоров используются канонический формат и формат Так-кера. Канонический формат (см. [13]) для вектора В имеет вид
К К
В = £и, ® V, ®...® ИЛИ Ь1] к = £Щ^^...^ (1)
* = 1 , = 1 где через ® обозначено прямое (внешнее, кронекерово) произведение векторов. Аналогичное представление для матрицы А записывается в виде
К
А = £ а, ® Ъ, ®...® с,. (2)
, = 1
При этом предполагается выполнение обозначенного ниже символом а/^ переформатирования
К
а(11...к)(у...к) ^ а({1)Ш)-(кк) = £ а(Н')8ЬШ)8 ■■■С(кк)з,
, = 1
которое переставляет элементы и приводит к d-мерному массиву с парными мультииндексами (/'/), (/У), ..., (кк), объединяющими в один "длинный индекс" индексы для соответствующих строчных и столбцовых мод.
Минимальное число слагаемых в возможных представлениях (1), (2) называют тензорным или каноническимрангом. Если указанные равенства справедливы с некоторой погрешностью, то говорят о канонической аппроксимации тензоров.
Разложение Таккера (см. [14]) для вектора В имеет вид
Г1 г2 гё
В = С Х1 и Х2 V. ха Ж = ££ ... £ gab.Ua ® ^ ®.® Wc, (3)
а = 1 Ь = 1 с = 1
где тензор коэффициентов О = [^аЪ.с] называется ядромразложения, а матрицы и = [иа], V = [уъ], ... ..., Ж = ^С] — модовыми факторами. Для матриц разложение Таккера выписывается аналогично (при той же схеме группировки индексов, как и в случае канонического разложения). Символом хк в (3) обозначено умножение тензора на матрицу по к-му направлению. В частности, запись С = = В х2 Мозначает Сук = ^»ту Ъу•. к. Пределы суммирования г1, ..., г11 называются модовымирангами тензора. Для простоты можно считать, что г1 = ... = г11 = г и Я = г, но предлагаемые ниже алгоритмы применимы, конечно, и в тех случаях, когда все модовые ранги и модовые размеры различны.
Для краткости канонический формат будем называть С-форматом, а формат Таккера — Т-форма-том. Для матриц А и В допустимы всевозможные сочетания форматов С и Т, но для С применяется только Т-формат. Мы налагаем это ограничение для того, чтобы уйти от сложной задачи переаппроксимации в С-формате. Среди множества предложенных для нее алгоритмов (см. [5], [13], [15]—[20]) нет ни одного достаточно надежного. Напротив, задача сжатия в Т-формате решается с помощью сингулярного разложения (см. [21], [22]), причем существуют и методы сублинейной по количеству данных сложности (см. [23], [24], об их развитии см. [10], [25]). Таким образом, мы будем изучать умножение матриц в случаях, которые символически можно обозначить как СТ —- Т и ТТ —- Т В каждом случае есть алгоритмы с одной и той же асимптотикой по п (см. [12]). Но практические вычисления требуют более тонкого анализа, включающего сравнение этапов, сложность которых определяется только рангом г.
Операция Модовое умножение Ортогонализация факторов Размер вектора, N
СУ 0(п 1о§ п)г2 0(пг 4) п
МУ 0(п2)г 2 0(пг4) п
ММ 0(п3)г 2 0(п2г4) п2
При записи формул числа обозначаются строчными буквами, векторы строчными жирными, матрицы заглавными, а тензоры жирными заглавными. Запись А =: ОЯ определяет правую часть как некоторое разложение матрицы А. Если индексы тензора сгруппированы в два мультииндек-са, то он естественным образом определяет матрицу и может применяться матричное обозначение. Пределы суммирования опускаются в предположении, что индексы /,у, ..., к пробегают значения от 1 до п, а индексы а, Ь, с, р, q, 5, I — от 1 до г. При суммировании по паре индексов полагаем, что каждый индекс независимо пробегает свой интервал. Фробениусова норма тензора и скалярное произведение тензоров определяются следующим образом:
у...к
Ц...к ^ < В С> = X ЬЧ-к£ч~к.
у...к
2. ОБЩАЯ СХЕМА УМНОЖЕНИЯ ТЕНЗОРНЫХ МАТРИЦ
Пусть A и B заданы в Т-формате (С-формат можно рассматривать как Т-формат с полностью диагональным ядром):
А = в иА) х2 х
3 • ■
ЖА) = х
рд...:
Жв) = X к
(А) и I
>рд...$ р
В = Н х, и{В) х2 VВ) х3. Тогда для C = A • B получаем представление в Т-формате С
.(В),
М) ,
,(В) ,
w,
(А)
(4)
аЬ ...с™ а
1 w.
(В)
X X
рд.. .5аЬ.. .с
... х,.. с (иА) • иав) )®с у^а) . уь в)).
/ (А) (В),
(w¡ • Wc ) =
(5)
- XX
рд...5аЬ ...с
с модовыми факторами и = ^ра], V = [^Ь], ..., Ж = и ядром F(G, H), имеющим размер г2 х ... ... х г2 и некоторую специальную структуру. Для получения и сжатия Т-формата для C в [12] предложен следующий алгоритм.
Алгоритм 1. Умножение с глобальной фильтрацией Шаг 1. Вычислить модовые факторы
2 кг и
Ьрд... 5 аЬ ...с ра
V
= в, Н) х, и х2 V хз ... х, Ж
(А) (В) и := и • и ,
ра р а
(А) (В)
V := Уд • УЬ ,
wгc := w
(А) (В)
wC
Сложность этого шага пропорциональна г2, зависит от п и выполняемой операции. В табл. 1 приведены оценки сложности для трех случаев: СУ-матрица A является ^-уровневой тёплицевой, B — вектор (свертка); МУ — умножение матрицы на вектор; ММ — умножение матрицы на матрицу. Этот шаг может быть очень дорогим, но обладает большим внутренним параллелизмом, так как все умножения могут быть произведены на разных процессорах независимо.
Шаг 2. Провести ортогонализацию модовых факторов (например, с помощью QR-разложения)
и =: ии', V=: ¡Г', ..., Ж=: (6)
где матрицы и, V', ..., Ж' имеют размер г2 х г2, а и, и,
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.