Тензоры: логарифмическая сложность

Иван Оселедец, Сколковский институт науки и технологий

oseledets.github.io, i.oseledets@skoltech.ru

О чем эта лекция

• Краткое введение в тензорные разложения
• Tensor-train разложение и его свойства
• Quantized-tensor-train (QTT)-разложение
• Теоретические результаты по аппроксимации решения эллиптических/параболических уравнених в 2д/3д в QTT-формате
• Многомасштабные задачи
• Первые численные примеры
• Другие приложения: многомерные интегралы в финансовой математике
• Небольшая дискуссия о будущем

Что мы хотим

Мы хотим эффективно решать многомерные дифференциальные уравнения с сильноменяющимися коэффициентов в областях сложной формы, например:

$$\mathrm{div}~k~\mathrm{grad}~u = f, \quad u_{\partial \Omega} = 0,$$

где $\Omega$ - некоторая область в 2D/3D.

Коэффициенты могут:

• Иметь высокий контраст $\frac{k_{\max}}{k_{\min}}$
• Сильно осциллировать (многомасштабные задачи)

Многомасштабные дифференциальные уравнения

Рассмотрим простейшее одномерное дифференциальное уравнение вида:

$$\mathrm{div}~k(x, x/\epsilon)~\mathrm{grad}~u = f, \quad u_{\partial \Omega} = 0,$$

Пример:

<img width=80% src='multiscale-sol.png' /img>

Уравнения с сильноменяющимися коэффициентами/полигональные области

Решение имеет особенность в углах,

поэтому использование конечных элементов низкого порядка

на регулярных сетках может иметь

невысокий порядок сходимости, например, $O(h^{1/4})$.

Известные подходы

Конечно, существует большое количество подходов для решения подобных задач:

1. Резко меняющиеся коэффициенты: адаптивные сетки, hp-методы
2. Многомасштабные задачи: осреднение (Бахвалов, Панасенко), многомасштабные методы конечных элементов

hp-методы (адаптивный выбор как шага сетки, так и порядка элементов) имеют оптимальную сложность;

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

Наш подход

Наш подход основан на принципиально новой идее.

1. На заданной области вводим логически прямоугольную сетку с $2^d \times K$ узлами, которая позволяет разрешить нужный масштаб задачи

2. Вектор решения представляем в виде многомерного массива (тензора) $2 \times 2 \times \ldots \times 2 \times K$ и приближаем в виде "тензорного поезда" (tensor train, TT-format)

   $A(i_1, \ldots, i_d) \approx G_1(i_1) G_2(i_2) \ldots G_d(i_d), $

   где $G_k(i_k)$ имеет размер $r_{k-1} \times r_k$ при фиксированном $i_k$.

Теория Казеева-Шваба

В 2015 году В. Казеев и К. Шваб в работе Quantized tensor-structured finite elements for second-order elliptic PDEs in two dimensions получили фундаментальный результат:

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

Впервые указан функциональный класс для таких решений

Оценка

Сходимость по числу степенй свободы в QTT-представлении имеет экспоненциальный вид:

$$\Vert u - u^{QTT}_h \Vert_{\mathbf{H}_1} \leq e^{-N^{1/k}}$$

Логически прямоугольная сетка

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

Многомасштабные задачи

В работе Kazeev, Oseledets, Rakhuba, Schwab QTT-FE Approximation For Multiscale Problems мы получили аналогичный результат об экспоненциальной сходимости для многомасштабных задач с квазипериодическим коэффициентом, основывая на теории осреднения:

Численный пример для многомасштабной задачи

$$ (a u')' = 1, \quad u(1) = u(0) = 0, $$

где $$a(x) = a_0(x) a_1\left(x/\delta)\right),$$

where $a_0 = 1 + x$, $a_1 = 2/3 (1 + \cos^2(2 \pi \frac{x}{\delta})).$

Как это работает

Выше мы описали теорему существования:.

Однако, нам надо решить линейную систему

$$Ax = f,$$

где $x$ на самом деле тензор $X$, представленный в TT-формате:

$$X(i_1, \ldots, i_d) = G_1(i_1) \ldots G_d(i_d).$$

Как мы можем это сделать?

Напомним сначала базовые свойства TT-формата.

ТТ-ранги - ранги матриц

Определим развертки:

$A_k = A(i_1 \ldots i_k; i_{k+1} \ldots i_d)$, $n^k \times n^{d-k}$

матрица

Теорема: существует TT-разложение с TT-рангами

$$r_k = \mathrm{rank} A_k$$

Доказательство

Доказательство очень простое и дает TT-SVD алгоритм

Устойчивость

На практике, точных малых рангов не бывает, поэтому есть оценка устойчивости:

If $A_k = R_k + E_k$, $||E_k|| = \varepsilon_k$ $$||\mathbf{A}-\mathbf{TT}||_F \leq \sqrt{\sum_{k=1}^{d-1} \varepsilon^2_k}.$$

Вычисления в TT-формате

Базовые операции выполняются за O(N).

Например, поэлементное умножение: $$C(i_1,\ldots,i_d) = A(i_1,\ldots,i_d) B(i_1,\ldots,i_d)$$

$$C_k(i_k) = A_k(i_k) \otimes B_k(i_k),$$ ранги умножаются, поэтому нужно уметь их уменьшать.

Округление

Округление решает следующую задачу: пусть дано TT-представление

$$A(i_1, \ldots, i_d) = G_1(i_1) G_2(i_2) \ldots G_d(i_d)$$

нужно найти другое с меньшей памятью, которое приближает $A$ с нужной точностью $\epsilon$.

Округление можно точно сделать за $\mathcal{O}(dnr^3)$ операций, используя линейную структуру формата .

Программная реализация

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

Для решения линейной системы:

$$x_{k+1} = R(x_k + B (f - A x_k)), $$

где $R$ - операция округления, а $B$ предобуславливатель.

А нужно ли предобуславливание?

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

Базовая идея: Пусть $A = A^* > 0.$

$$(Ax, x) - 2 (f, x) \rightarrow \min,$$

и минимум считаем по множеству всех тензоров ограниченного ранга (многообразию).

Метод переменных направлений

Простейшим оптимизационным методом является метод переменных направлений

$$X = G_1(i_1) \ldots G_d(i_d),$$

если мы фиксируем все $G$ кроме $1$, функционал становится квадратичным, и минимизация сводится к решению (небольшой) линейной системы.

Оптимизационные методы на многообразиях

Более сложные методы основаны на римановых методах минимизации на многообразиях.

Важна кривизна многообразия, однако мы установили, что множество TT-тензоров ведет себя как линейное многообразие и кривизна не играет особой роли.

Типичная картина сходимости

Здесь синия линия - сходимость без проекции на многообразие, зеленая - со сходимостью, ступенчатая - "нормальная компонента" решения.

Общая идея

Созданы (и постоянно улучшаются) оптимизационные алгоритмы типа "черный ящик" для минимизации функционалов на TT-многообразии Поэтому, объединяя:

• теорию о существовании решения на логически прямоугольных сетках
• методы построения таких сеток
• пакеты по TT-оптимизации,

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

Grand Challenges

Есть четыре больших направления исследований на ближайший год:

• Разработка специальных дискретизаций на сверхмелких сетках: обычные схемы становятся неустойчивыми при $h \sim 10^{-4}$. Простое решение - использовать вейвлет-конечные элементы, так как нет необходимости использовать локальные базисные функции. Для одномерных задач мы успешно решаем задачи на сетках порядка $2^{120}$.

• Оптимизация алгоритмов решения задач оптимизации на TT-многообразия (теория сходимости, методы ускорения)

• Создание пакета, объединяющего построения логически прямоугольных сеток (вычислительная геометрия, сотрудничество с ETH и NYU Courant)

Приложения

Создаваемый подходв позволяет решать огромное количество прикладных задач "прямым" моделированием:

• Композитные материалы с дефектами
• Фотонные кристалы, оптические волноводы (уравнения Максвелла).
• Акустика, задачи рассеяния

Во всех случаях можно обобщить теоретические результаты Казеева-Шваба.

Выводы

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

Ссылки

• oseledets.github.io : статьи, ссылки
• http://github.com/oseledets/TT-Toolbox - пакет на матлабе
• http://github.com/oseledets/ttpy - пакет на Питоне