Алгоритмы и применения тензорных разложений для численного решения многомерных нестационарных задач тема автореферата и диссертации по математике, 01.01.07 ВАК РФ

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

Долгов Сергей Владимирович

Алгоритмы и применения тензорных

разложений для численного решения многомерных нестационарных задач

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

Автореферат

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

1 8 СЕН 2014

Москва 2014

005552621

005552621

Работа выполнена в федеральном государственном бюджетном учреждении науки Институте вычислительной математики Российской академии наук (ИВМ РАН).

Научный руководитель: член-корреспондент РАН,

доктор физико-математических наук, профессор Тыртышников Евгений Евгеньевич.

Официальные оппоненты: Истомин Яков Николаевич,

доктор физико-математических наук, профессор. Федеральное государственное бюджетное учреждение науки Физический институт им. П.Н. Лебедева Российской академии наук, зам. руководителя отдела теоретической физики.

Третьяков Алексей Анатольевич, доктор физико-математических наук. Федеральное государственное бюджетное учреждение науки Вычислительный центр им. A.A. Дородницына Российской академии наук, ведущий научный сотрудник.

Ведущая организация: Федеральное государственное бюджетное учреждение науки Институт прикладной математики им. М.В. Келдыша Российской академии наук.

Защита состоится _ 2014 г. в ^"час. ОС? мин. на заседании

диссертационного совета Д 002.045.01 при федеральном государственном бюджетном учреждении науки Институте вычислительной математики Российской академии наук (ИВМ РАН) по адресу: 119333, г. Москва, ул. Губкина, 8, ауд. 729.

С диссертацией можно ознакомиться в библиотеке и на сайте ИВМ РАН http://www.inm.ras.ru.

Автореферат разослан

2014 г.

Ученый секретарь ___, ^

диссертационного совета

д.ф.-м.н. ' Бочаров Геннадий Алексеевич

Общая характеристика работы

Объект исследования и актуальность работы

Диссертация посвящена численному решению многомерных задач методами тензорных разложений. Например, мы хотим одновременно описать несколько взаимодействующих под действием случайных сил тел. Совместная функция плотности вероятности задается в ¿-мерном пространстве, и, как правило, требуют пл значений после дискретизации, где (1 составляет количество всех координат всех частиц, а п — число степеней свободы по каждой размерности. Под "размерностью" мы и будем иметь в виду количество конфигурационных координат с1 в пространстве состояний системы. В принципе, даже случаи й = 3 или д, = 2 можно рассматривать как "многомерные", если п очень велико. Ситуация становится еще более драматичной, если физическая или математическая модель предусматривает работу с размерностями порядка десятков, сотен и более. Если исключить тривиальный случай, когда п = 1 для большинства координат (в этом случае нет смысла рассматривать соответствующие переменные вообще), экспоненциальный рост вычислительной сложности с размерностью делает невозможными непосредственные расчеты на любой суперкомпьютере. Например, в квантовом мире, частицы со спином 1/2 (в определенных условиях, например в магнитно-резонансных экспериментах, рассмотрение электронов и ядер может быть ограничено спиновыми эффектами) обладают только п — 2 состояниями для каждой частицы, "спин вверх " и "спин вниз". Однако, простая линейная цепочка из (I = 100 взаимодействующих спинов (что рассматривается как модельная задача в квантовой физике) описывается уже волновой функцией с 2100 ~ Ю30 неизвестными.

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

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

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

ты (например, сингулярное разложение) для выделения редуцированного набора параметров.

Ключевым моментом в разделении переменных является представление многомерной функции (или ее дискретного аналога, тензора) в виде произведения одномерных объектов, т.е. x(ii,..., id) = a;'1'(¿i)a;'2'(z2) ■ ■ ■ x^{i¿¡). Это разложение в прямое произведение редко выполняется точно. Обычно используют определенную комбинацию тензорных произведений, но конкретные численные алгоритмы нахождения такой комбинации существенно отличаются по своим свойствам.

Так, "жадные" методы тензорных аппроксимаций часто стагнируют на каком-то (неудовлетворительно большом) уровне ошибки. Методы оптимизации общего вида страдают от некорректности задачи аппроксимации тензора размерности больше двух суммой R прямых произведений. В качестве альтернативы были предложены рекуррентные двумерные разложения: Tensor Tree Networks, Hierarchical Tucker, Matrix Product States/Tensor Train. Последний представляет данные в виде полилинейной комбинации трехмерных массивов, и требует объема памяти 0(dnr2), что дает возможность избежать проклятия размерности, если ранг г не очень большой.

Любая тензорная структура требует численных методов для аппроксимации данных и других операций, таких как решение линейных систем. Одними из самых эффективных явились методы переменных направлений, а точнее семейство алгоритмов Density Matrix Renormalization Group, разработанных в сообществе квантовой физики для моделирования волновых функций спиновых систем в форматах тензорных произведений. Тем не менее, даже функции ошибки вида Ца; — £*||2, по отношению к элементам тензорного разложения, могут иметь многочисленные локальные минимумы, и алгоритмы DMRG с высокой вероятностью возвращают локальный, но не глобальный минимизатор ошибки.

Другим подходом является адаптация классических итерационных алгоритмов (например, GMRES) к хранению и вычислению векторов в виде тензорных произведений. Этот путь, однако, также не является достаточно надежным, поскольку чем ближе приближенное решение к точному, тем более сложной и осциллирующей является структура ошибки и других вспомогательных (в т.ч. Крыловских) векторов, которые требуют либо грубого приближения, либо большого объема памяти для хранения их тензорных форматов.

Цели диссертационной работы

1. Разработка нового вычислительного метода для решения больших систем линейных уравнений с представлением данных в формате тензорных произведений. Теоретический анализ его свойств сходимости.

2. Применение и проверка нового метода на задачах расчета стохастической химической кинетики (основное кинетическое уравнение) и задаче моделирования Фарлей-Вунемановской неустойчивости в плазме ионосферы Земли (уравнение Власова).

3. Анализ результатов, сравнение нового алгоритма с ранее существующими техниками, сравнение модельных численных данных с имеющимися в литературе экспериментами.

Основные положения диссертации

Главным результатом данной диссертации является вычислительный метод, который сочетает в себе сильные стороны как оптимизационных тензорных алгоритмов переменных направлений, так и классических приближенных итерационных схем. В процессе DMRG итерации, мы расширяем тензорный формат решения с помощью тензорного формата приближенной невязки. Это обеспечивает замечательную взаимную поддержку классических методов и методов переменных направлений (DMRG). Так как решение адаптируется в процессе DMRG оптимизации по элементам формата, даже очень грубые приближения невязки (причем без других Крыловских векторов) дают высокую точность решения. С другой стороны, подключение информации о глобальной невязке в локальных шагах процесса переменных направлений обеспечивает последний правильными градиентными направлениями, помогая ему избегать локальных минимумов. Новый AMEn (Alternating Minimal Energy) метод обладает доказанной геометрической сходимостью с точки зрения глобальных элементов тензора, аналогично методу градиентного спуска. Примечательно, что практическая скорость сходимости оказывается намного быстрое, чем теоретические оценки на основе градиентного спуска, что делает этот метод надежным даже для несимметричных линейных систем, наподобие метода полной ортогонализации.

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

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

Обе упомянутые модели описывают эволюцию функции вероятности во времени. Таким образом, в качестве третьего приложения, мы покажем, как схемы интегрирования по времени, например, Эйлера или Кранка-Николсон, могут быть ускорены методами тензорных произведений. Мы рассматриваем переменную времени в качестве дополнительного измерения. Затем используем тензорный формат QTT, что дает логарифмическую вычислительную сложность в зависимости от количества временных шагов.

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

мента для различных задач, которые в принципе подходят для концепции разделения переменных.

Достоверность научных положений и методология

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

Научная новизна

Предложенный метод АМЕп является первым алгоритмом для решения системы уравнений в представлении тензорными произведениями, который обладает теоретически доказанной оценкой глобальной сходимости, и при этом является эффективным на практике. В основе АМЕп алгоритма лежит сочетание метода градиентного спуска и метода переменных направлений (ВМГЮ). Однако, в отдельности метод градиентного спуска сходится весьма медленно и требует достаточно точной аппроксимации невязки (градиента), что приводит к большой вычислительной сложности. Отдельные шаги методов ВМГЮ весьма эффективны, но итоговый результат может давать очень плохое приближение для точного решения.

С помощью нового метода было впервые получено полное решение основного кинетического уравнения для Л-фага с высокой точностью. Классический способ аппроксимации решения основного кинетического уравнения состоит в усреднении большого набора случайных реализаций. Он хорошо применим в случаях, когда требуется невысокая точность. С ростом числа реализаций М метод сходится с достаточно медленной скоростью 0(М~1/2), что отрицательно сказывается на вычислительной сложности в высокоточных расчетах. В качестве альтернативы, были предложены и методы тензорных представлений. Однако ранее существовавшие алгоритмы позволяли моделировать только небольшие промежутки времени, что непригодно для реальных задач в системной биологии.

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

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

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

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

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

Апробация результатов диссертации и публикации

Основные результаты работы апробированы в следующих докладах па конференциях и семинарах.

[1] Solution of the chemical master equation by the separation of variables and alternating optimization methods. European Conference on Mathematical and Theoretical Biology, Gothenburg, June 16, 2014.

[2] A new family of the alternating linear optimization schemes in tensor product representations. Seminar of the group Computational Methods in Systems and Control Theory, Max Planck Institute Magdeburg, March 04, 2014.

[3] Alternating minimal energy methods for linear systems in higher dimensions. Part II: Faster algorithm and application to nonsymmetric systems. ENUMATH Conference, EPFL Lausanne, August 26-30, 2013.

[4] Alternating minimal energy methods for linear systems in higher dimensions. Part II: Faster algorithm and application to nonsymmetric systems. MAFELAP, Brunei University, London, June 11, 2013.

[5] Fast adaptive alternating linear schemes in higher dimensions. Part 2: eigenproblems. NASCA Conference, University of Calais, June 24, 2013.

[6] Fast adaptive tensor product approach to eigenvalue problems in higher dimensions. Seminar of the Department of Chemistry, University of Southampton, June 27, 2013.

[7] Alternating minimal energy methods for linear systems in higher dimensions. Part II: Faster algorithm and application to nonsymmetric systems. Smss numerics colloquium, EPFL, Lausanne, April 05, 2013.

[8] Fast adaptive alternating linear solvers. Implementation hints and application to Fokker-Planck and master equations. Workshop on algorithms for high-dimensional problems in quantum chemistry, University Southampton, February 26, 2013.

[9] Alternating minimal residual methods for linear systems in higher dimensions. Part II: heuristics and experiments. 29 GAMM Seminar on Uncertainty Quantification, MPI MIS, Leipzig, January 22, 2013.

[10] Advanced tensor representation and solution techniques with application to Fokker-Planck and master equations. CMAM-5, Humbolt University, Berlin, August 16-18, 2012.

[11] Advantages and difficulties of use of tensor methods in solution to the Fokker-Planck equation. 28 GAMM Seminar on Analysis and Numerical Methods in Higher Dimensions, MPI MIS, Leipzig, January 16-18, 2012.

[12] A gray-box DMRG algorithm for tensor structured solution to linear systems. 17th Conference of the International Linear Algebra Society, University Braunschweig, August 26, 2011.

[13] On a solution to a parabolic equation in the QTT format using the DMRG approach. 4£/i Workshop on High-Dimensional Approximation, University Bonn, June 27, 2011.

[14] Use of the DMRG scheme for structured linear system solution. 3rd International Conference on Matrix Methods in Mathematics and Applications, ИВМ PAH, Москва, 24 июня, 2011.

[15] Linear solvers in TT formats. 27th GAMM-Seminar on Approximation of Multiparametric functions, MPI MIS, Leipzig, January 26-28, 2011.

[16] TT-GMRES: о решении систем линейных уравнений в тензорных форматах. 53 паучпая конференция МФТИ, Москва, Ноябрь 26-28, 2010.

По результатам работы опубликовано 6 статей в международных рецензируемых журналах, и 1 статья в материалах конференций [3].

Структура диссертации

Диссертация состоит из введения, пяти глав, заключения, списка обозначений и списка литературы. Текст изложен на 161 странице, содержит 25 рисунков. Список литературы включает 253 наименования.

Благодарности

Эта диссертация содержит основные результаты моей научной работы с 2011 по 2014 год. Прежде всего, я хотел бы выразить благодарность моему научному руководителю, профессору Евгению Евгеньевичу Тыртышникову за внимательное руководство и неизменную поддержку в институте вычислительной математики.

Я глубоко признателен моим коллегам Ивану Валерьевичу Оселедцу и Дмитрию Валерьевичу Савостьянову за очень плодотворное сотрудничество и обмен

идеями. Их высокий профессионализм и поддержка были неоценимы в процессе исследования многомерных областей математики.

Я особо благодарен гостеприимству математического института им. Макса Планка в Лейпциге, и лично профессору Борису Николаевичу Хоромскому за вдохновляющую атмосферу и уникальные условия, великолепно способствующие научной деятельности.

Я высоко ценю замечательное и дружелюбное сотрудничество, установившееся с нашими коллегами в различных группах, занимающихся как прикладными задачами, так и разработками вычислительных методов, в особенности профессору Илье Купрову, University of Southampton, профессору Александру Павловичу Смирнову, факультет вычислительной математики и кибернетики МГУ, профессору Геннадию Алексеевичу Бочарову, институт вычислительной математики РАН, профессорам Томасу Шульте-Хербрюггену и Томасу Хукле, Technical University Munich, профессорам Удо Райхлю и Питеру Беннеру, институт сложных динамических систем им. Макса Планка, Магдебург. Объективная оценка текущей работы и дальнейших направлений исследований были бы невозможны без мотивирующих дискуссий о последних математических инструментах и интригующих практических приложениях.

Содержание работы

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

В первой главе мы описываем три основных приложения: уравнение Власова для Фарлей-Бунемановской неустойчивости, основное кинетическое уравнение для стохастической химической кинетики, и схема одновременной пространственно-временной дискретизации динамических уравнений.

Фарлей-Бунемановская (ФБ) неустойчивость возникает в слабо ионизированной плазме Е-области ионосферы Земли. Эта неустойчивость порождается в плазме с замагниченными электронами и незамагниченными ионами в электрическом поле, направленном перпендикулярно геомагнитному полю. В Е-области электроны подвержены заметному влиянию геомагнитного поля, в отличие от ионов, которые не замагничиваются из-за частых столкновений с нейтральными частицами газа. В результате, посредством скорости дрейфа, обусловленной электрическим полем, распределение электронов по скоростям сдвигается относительно ионного распределения. Соответствующие условия для возникновения неустойчивости возникают в экваториальных и полярных зонах Е-области ионосферы Земли на высотах порядка 90-130 км, где нестабильность проявляется как низкочастотные колебания плазмы с длинами волн порядка метра.

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

полю). Но поведение ионов требует кинетического описания, которое дается четырехмерным уравнением Власова (две пространственные, две скоростные координаты):

дф дф дф дфдф ( еЕ0 дф\ дф от ах ду ах ои ау) оти

где ф = ф(х, у, V, ш, ¿) является функцией распределения ионов по скоростям v,w в

каждой пространственной точке х, у во время t, ф обозначает электростатический

потенциал, возникающий из-за неравновесности концентраций ионов и электронов, е2

Аф =-— "1)1 Е0 это внешнее электрическое поле, е - заряд электрона, £о

£отп^п

- диэлектрическая постоянная, пц - масса иона, !/,-„ - частота столкновений ионов с нейтральными частицами, = \/Т^/гщ - характерная тепловая скорость ионов с температурой Т.В правой части уравнения (1) стоит упрощенный оператор столкновений Bhatnagar-Gross-Kгook, где

7Т" ( \ [

= ~ схр I----I , тц = ф(х, у, V, (2)

определяет правило вычисления нейтрального распределения: оно соответствует распределению плотности ионов в пространстве, но предполагается максвеллов-ским по скоростям.

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

= ф^,лф^(к)ф^{т) = £ (3)

с разумными ТТ рангами, и сокращения затрат памяти в десятки раз.

Для дискретизации по пространственным переменным х, у, а также по скоростям V, и> в уравнении (1) мы используем разностные схемы, принимая во внимание периодические граничные условия. Хотя областью определения для скоростей является вся плоскость, на практике функция распределения ф в (1), будучи возмущенным распределением Максвелла, быстро убывает с ростом V и ш, и мы можем ограничить область до квадрата (и, -ш) е [—г>„шх, vrлaxY, накладывая периодические граничные условия.

Для решения нелинейных уравнений модели, мы используем расщепления по физическим процессам и координатам:

1. диффузия электронной плотности;

2. конвекция электронной плотности;

3. Уравнение Пуассона для электростатического потенциала.

4. конвекция ионного распределения;

5. релаксация ионного распределения за счет столкновений.

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

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

Динамическую эволюцию стохастической системы реакций можно рассматривать как дискретный марковский процесс. Предположим, что (1 различных активных химических соединений 5],..., в хорошо перемешиваемой среде могут реагировать посредством М реакционных каналов. Каждый канал задается сте-хиометрическим вектороль гт е Ъл, и функцией скорости и!т(г) : —> К+, т = 1,..., А/, Е+ = {х € К : х ^ 0}, так что т-я реакция может быть записана в классическом виде

Чтобы ввести стохастическое описание, обозначим состояния мультииндексом г = (¿1,... , г^), и всегда будем иметь в вид}' числа молекул, так что является неотрицательным целым числом, ¿д.. 6 ({0} и М). Вероятностная роль функции скорости следующая: для бесконечно малого интервала времени <й, 1К"'(г,<,гй) = ыт{г)сИ является вероятностью того, что при числе молекул г в момент времени ¿, в последующем интервале , Ь + <й) в системе произойдет одна реакция по каналу т. Количественной характеристикой состояния г является вероятность того, что количества молекул веществ ¿>1,..., Ба в момент времени 4 принимают значения ¿1,..., i¡¡, Ф(г,() : ({0} иМ)<ги [0,Г] -»• Теперь, принимая (й достаточно малым, так что вероятность того, что больше, чем одна реакция будет происходить в интервале [Ь, < + (й) пренебрежимо мала, можно написать распределение в конце интервала, + в£), используя законы сложения и умножения вероятностей для независимых и взаимоисключающих событий. Переходя к пределу ей —> 0, получаем основное кинетическое уравнение:

а,"5, + • • ■ + аГЗ, К + + • ■ • + (ау +

Потенциально возможны любые числа молекул, т.е. уравнение (4) является бесконечномерным ОДУ. Разумеется, для проведения численного моделирования нам нужно ограничить его до конечной задачи. Так называемый алгоритм проекции в конечное пространство (Finite State Projection, FSP) использует тот факт, что очень большие количества молекул имеют очень малую вероятность возникнуть за конечное время: Ф(г, £) —> 0, |г| -> оо. Таким образом, мы полагаем что каждый ik лежит в конечном диапазоне, ¿ц- = 0,... ,пц — 1, и выбираем пдостаточно большими, так что, например, при г* > пвероятность Ф(г,4) меньше машинной точности, и можно пренебречь ошибкой, внесенной при ограничении пространства состояний. Однако, даже если каждый Пк = 0(п) не превышает десятков, общее число степеней свободы растет как nd, и сжатое хранение и обработка для Ф принципиально необходимы.

Схемы интегрирования по времени, например классические методы Кранка-Николсон и Эйлера для линейных ОДУ, могут быть переписаны в виде блочных систем линейных уравнений. Пусть в качестве базисной используется схема Кранка-Николсон,

(l+^A)Xpn = (/-|л)хр + |(/г, + /р+1), Р = 0.....JVt - 1, (5)

где хр ~ x(tp), fp = /(ip), tp = pSt задает сетку по времени, тогда мы можем собрать все временные шаги в одну глобальную линейную систему:

-1 + ЦА 1 + %А

-1 + ЦА 1 + ЦА

где д = [/о + fi fi + /2 ••• +/л',]Т. Конечно, эта система не несет смысла,

если мы будем держать ее в стандартном виде - пошаговое решение (5) является оптимальным методом для двухдиагональной матрицы. Тем не менее, предположим что можно вычислить ее решение более эффективно в структурированном представлении. Тогда становится заманчиво получить сразу всю историю динамики, и использовать очень мелкие шаги по времени, чтобы гарантировать достаточную точность для всех спектральных компонент матрицы А. С использованием приближений тензорными произведениями, мы можем ввести дополнительные виртуальные переменные, и в конечном итоге аппроксимировать и пространственную, и временную части системы (6) с логарифмической редукцией вычислительной сложности, 0(log(Ar) log(Ar()), которая действительно подтверждается в численных экспериментах.

Глава 2 посвящена описанию различных методов разделения переменных. Дается обзор существующих тензорных разложений: канонического (CP), формата Таккера, и рекуррентных представлений иерархический Таккер и MPS/TT, а также их свойств и связанных с ними алгоритмов. Рекуррентные тензорные представления являются эффективным инструментом сжатия данных. На первом шаге построения такого формата мы используем группировку индексов, чтобы уменьшить число переменных формально (но получить большее число степеней свободы

XI Хо i^Ax 0

Х2 = 0

0

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

В особенности, мы сфокусируемся на MPS/TT представлении, в котором размерности группируются в соответствии с линейным деревом. Практический опыт позволяет полагать, что наиболее важным является порядок индексов (а он может быть определен посредством физических соображений, или с использованием адаптивных алгоритмов), тогда как для группировки размерностей в большинстве случаев достаточно линейного дерева. Его простота обеспечивает определенную элегантность обозначений, что особенно важно в многомерных задачах из-за обилия индексов.

Говорят, что тензор х представлен (или аппроксимирован) в формате Matrix Product States (MPS), или Tensor Train (TT), если выполняется

*(«!,...,»,)« £ (7)

«1 - 1

где a;'*' 6 ct_1xntxrt назьшаются TT ядрами (блоками, или факторами), диапазоны Гц, = rk{x) ^ г(х) ранговых индексов ак = 1,..., г к называются ТТ рангами. Для однообразия записи можно полагать го = rd = 1. Любопытно, что это представление было независимо открыто несколько раз: как MPS в квантовой физике с 1980'х, и как ТТ в 2009 в вычислительной линейной алгебре. Асимптотические затраты памяти для ТТ формата составляют О(dnr2). Для небольших диапазонов исходных индексов ik (т.н. модовых размеров п) и более высоких рангов, ТТ формат содержит меньше неизвестных, чем более сложные древовидные разложения (последние имеют асимптотическую сложность не ниже 0(í-3)), что делает его особенно эффективным для описания, например, спиновых систем.

ТТ формат для матриц, также известный как Matrix Product Operator, вводится с использованием перестановки индексов,

A{4,...,id,ju...,jd)= Y, (8)

С использованием дистрибутивности, произведение матрицы на вектор пишется независимо для каждого ТТ блока: результат получается в аналогичном ТТ формате у(г) = 2/(1)(¿i) • • • yw(id), где

т.е. = (9)

h

причем pk = = 1,... ,rk(A)rk(x), к = 0,... ,d.

Для единого описания форматов (7) и (8), удобно использовать следующее Определение. Пусть даны ТТ форматы {i(t)(ú)} или {A(t,(¿t.jt)}. Векторное ТТ отображение раскрывает ТТ представление в полный тензор следующим образом:

т(х<р),.. .,х<«>) : {хЮ}чк=р -у х.....*> е cr"-lXn"-"'><r', где

при 1 < р ^ 9 < Для граничных случаев мы можем дополнительно обозначать _ _ х(<«+1) £ <£"1-",хг11 т(р,-.<') _ хОР) = е ^-^»г -»^

Матричное ТТ отображение т{А^р\..., Л'®') аналогично соответствует (8).

Другие полилинейные алгебраические операции, например, г — х + у, могут быть также перенесены с элементов исходных тензоров на элементы их тензорных форматов,

,("(¿0 = [*<»&) »(«(м)] , z«\id) =

xV\id)

*w(u) -

x<"Hik)

для к = 2,..., d — 1. Все операции обладают линейной сложностью по размерности d, и полиномиальными сложностями по п и г.

Главным преимуществом тензорных деревьев по сравнению с простой суммой R прямых произведений является процедура рекомпрессии или округления. Алгебраические операции могут увеличивать ТТ ранги, при этом они могут получаться завышенными для данных тензоров и точностей. ТТ округление позволяет понижать ранги до квазиоптимальных значений при заданной точности. Эта процедура использует известные QR и сингулярное разложения матриц, и может быть выполнена за 0(dnr3) операций.

Отдельно стоит отметить понятие тепзоризации, или квантизации, предложенной Оселедцем и Хоромским. Рассмотрим для краткости одномерный вектор, х = [х(г)]. Пусть число допустимых значений для г составляет степень 2, т.е. п = 2L. Рассмотрим позиционную запись г в двоичной системе исчисления, т.е.

L

< = $>■ 2'"\ г, €{0,1}. (И)

1=1

Это соответствует перегруппировке вектора в тензор с L виртуальными (кванти-зованнъши) размерностями. Теперь можно применить ТТ разложение:

Если ТТ ранги этого тензора малы, требуемый объем памяти составляет примерно логарифм от исходного, О(L • 2 • г2) = O(logra). Это позволяет добиться большей эффективности сжатия данных, в частности, в "маломерных", двух- и трехмерных задачах. Такая разновидность ТТ формата была названа Quantized Tensor Train (QTT) формат. Название Quantized связано с термином "квант", минимально возможный долей информации, получаемой при задании каждой цифры ц.

При работе с малопараметрическими представлениями данных, неизбежно задается вопрос о конкретных значения ТТ (или QTT) рангов, г. В конце главы предлагается новое комбинированное тензорное представление (QTT-Tucker), объединяющее аналитические и практические преимущества Таккеровского, ТТ и QTT форматов. В основе нового формата лежит разложение Таккера:

Х(н.....id) = J2 (12)

Проклятие размерности снимается за счет хранения Таккеровского ядра в ТТ формате:

«(е)(7ь-,7-)= Е ^^О^ЦЫ-'-х^Ы. (13)

Наконец, если Таккеровские факторы х^({к) из-за больших размеров пк занимают слишком много памяти, для индекса можно ввести тензоризацию:

= Е (н)

На практике, этот формат может приносить дополнительное сокращение сложности, поскольку С^ТТ ранги в Таккеровских факторах требуют меньших значений для поддержания той же точности аппроксимации, чем соответствующие С^ТТ ранги в глобальном линейном формате. Еще одной интересной особенностью нового формата является то, что он может наследовать оценки рангов из формата Так-кера, которые могут быть гораздо легче доказуемы, чем оценки в ТТ формате, за счет использования полиномиальной интерполяции. Вычислительные алгоритмы для формата С^ТТ-ТЬскег удобно описывать с точки зрения двухуровневого ТТ представления. Это позволяет использовать ТТ версии многих методов (например, переменных направлений).

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

В качестве наиболее примечательных примеров, следует отметить следующие.

Матрицы конечных разностей в С^ТТ формате: пусть Зп есть Жорданов блок размером п х п, б,, = /„ — ^ есть матрица жесткости в схеме Кранка-Николсон (6), и Мп = /„ + ,7,7 есть матрица масс. Если п = 2то можно записать следующее С^ТТ представление с рангами 2:

Мп{ьз)

где в =

1

-1

и М =

М,

"1 0

1 1

Jik.ll

ТТ

1 г

суть элементарные (2 х 2) матрицы конечных

разностей и масс.

Матрицы перехода для ионного уравнения Фарлей-Бунемановской неустойчивости. Например, одномерная матрица интегрирования на один шаг по времени конвекционной части уравнения Власова (1) пишется в виде

Мс = а_25_2 + + а050 + а^ + о252,

где 5_2, ...,52 суть матрицы периодического сдвига, а коэффициенты а_2,..., а2 определяют веса интерполяции с характеристики. Однако, в нашем случае каждый вес а; параметризовал значением скорости и. То есть, действие а; задается

диагональной матрицей

Л; = [I ® /] ® diag(a;(u)) ® I,

и в итоге для матрицы перехода в исходном ионном уравнение получаем ТТ представление ранга 5:

Mx(St) = [S_2 ® /] ® diag(o!_2(ti)) ® / + [S_i ® /] ® diag(a_1(î))) ® I + [50 ® J] ® diag(a0(u)) ® I

+ [S'i ® /] ® diag(ai(t))) ® / 4- [S2 ® I] ® diag(a2M) ® I,

где Кронекеровские произведения nom y (в квадратных скобках) вычисляются явно (мы не разделяем пространственные переменные), тогда как остальные подразумеваются неявно, и сомножители хранятся отдельно, как ТТ блоки.

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

Лемма. Пусть даны элементарные векторы или матрицы Ek(ik), F£(ik), F*+1(ù.) для k = 1,..., cl. Сумма попарных произведений

)i2La®iî® ( (g) Ej

d \ d /k-2 \j=2 1 k=2 \i=l

(15)

обладает точным ТТ разложением Я(г) = #(1)(ii) • • • Я(<г)(г^) ранга 3, вида

H(I1(i,) = [Wi) Fîdi)], H«\id) =

гад Ft+1(h)

0

Fdd{id) ESd)

H(k\ik) =

*î(ik)

Ek{iK),

(16)

если k = 2,..., d — 1.

Разложение Таккера (12) также выполняется с рангами не более 3.

В главе 4 изложены методы решения линейных уравнений в тензорных форматах: классические итерационные методы с тензорными округлениями, итерационные методы оптимизации переменных направлений, разработанные в квантовой физике (ОМГЮ) и вычислительной математике (метод наименьших квадратов в переменных направлениях), в также новый улучшенный алгоритм (АМЕп), которая дополняет схему переменных направлений классическим градиентным шагом с использованием глобальной невязки системы. Мы анализируем его сходимость, доказываем геометрическую скорость сходимости, и обсуждаем некоторые практические аспекты.

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

1. ошибка Ех,(.т) = - х||2 = J7iI.,

2. энергия Зд,ь(х) = ||.т, - х\\2А = (х,Ах) - 2Re(x,b) + const, где b = Axt, и

3. невязка ^¡.(х) = \\Ах - Ь||2 = Лд-л,л-ь(х) + const.

Эти три функции связаны с решением некоторой (может быть, с единичной матрицей) линейной системы. Функция энергии определяется для симметричной положительно определенной (SPD) матрицы, А = А" > 0, с А-скалярным произведением, и А-нормой, введенными обычным образом: (х,у)А = (х,Ау), |[ж||д = \/(х,х)А.

Одним из наиболее простых тензорных оптимизационных методов является метод наименьших квадратов в переменных направлениях (ALS), так же названный линейная схема в переменных направлениях, и совпадающий с (одноблочным) Density Matrix Renormalization Group (DMRG) подходом из квантовой физики. Оригинальная DMRG схема была предложена для вычисления основного состояния системы путем минимизации отношения Рэлея. Позднее она была применена для решения SPD линейной системы Ах = Ь, где А и 6 заданы в ТТ формате, путем минимизации функции энергии.

Итак, мы ограничиваем оптимизацию Лл,ь(х) на векторы х — т(.г'(1),... ,x(d)), представленные в ТТ формате с фиксированными ТТ рангами г = (п,..., rj_i), и выполняем фактические вычисления посредством последовательности микрошагов, т.е. последовательных оптимизаций по ТТ блокам х'*'. Каждая такая локальная задача ставится следующим образом:

■u(k) = argminJA,,,(T(x(1),...,j;(d))) по е cri-lX""<ri. (17)

jCO

ТТ блок ж'*' затем заменяется на и(к>, и процесс переходит к следующему ядру: обычно блоки пересчитываются в последовательном проходе по размерностям, например к = l,...,d (прямой полу-проход), или к = d,..., 1 (обратный полупроход), и так далее вплоть до сходимости х.

Используя определение ТТ отображения (10) для частичного ТТ представления, мы можем определить фрейм-матрпцу:

Хфк = х«*> ® 1Пк ® (х(>ч)Т е C",-"J><r'-'"trt. (18)

Она задает соотношение линейности ТТ формата, х = Х^х^К Это позволяет переписать (17) в виде решения линейной системы

„(*) _ л~1Ь,

к к' ГШ

Ак = Х^кАХФк Ък = X^kb е СТк-1"кТк. К '

Недостатком одноблочного алгоритма DMRG является то, что ТТ ранги остаются постоянными в процессе расчета. Мы должны угадывать правильные значения рангов решения априори, что может быть трудно; Кроме того, сходимость даже к квази-оптималыюму решению при данных ТТ рангах может быть очень медленной. Для адаптации ТТ рангов в процессе расчета был разработан двухблочный DMRG. Данный метод работает с вектором в следующем виде:

х = т(х<",... .¡Е«*-»), х<*+2\ ..., хЮ), (20)

где = = t(xw,x<-m)) пишется в соответствии с определением

(10). Шаг локальной оптимизации проводится аналогично (17), но по элементам суперядра ¡rt*-*"1-1) следующим образом:

u("+4 = arg min J^,,бk.k+Ax(hM1)) n° xik,k+1) (21)

Решение системы (21) затем разделяется обратно на ии и'*"1"11 для вос-

становления первоначальной ТТ структуры с помощью матричного малорангового разложения. Так как элементы и<-к-к+1) отличаются от элементов предыдущего приближения, ранг г^ может быть изменен. Тем не менее, из-за локального характера оптимизации, методы DMRG (даже двухблочные) могут терять важную часть информации о направлении к точному решению, и стагнировать в каком-то локальном минимуме, далеко от желаемого порога.

Одноблочный DMRG, аналогично методу градиентного спуска, является вариационным для функции энергии: на каждом шаги, эти методы ищут минимизатор функции энергии с ограничениями. В этом смысле DMRG метод аналогичен классическим проекционным методам, тому же алгоритму градиентного спуска или GMRES. Разница состоит в том, что фрейм-матрицы текущего решения не приближают Крыловские базисы, и, следовательно, не гарантируют сходимость метода. В данной диссертации, мы предлагаем новый метод (АМЕп), который сочетает классические алгоритмы и схемы переменных направлений, путем обогащения фрейм-матриц информацией Крыловского типа, а точнее, фрейм-матрицей невязки. По аналогии с алгоритмом минимальных невязок (Minimal Residual, MR), мы можем назвать метод градиентного спуска также методом минимальных энергий (Minimal Energy, MEn). Новый метод решает линейную систему в ТТ формате методом DMRG, дополненном шагом вида градиентного спуска. Это и послужило причиной названия АМЕп: Alternating Minimal Energy.

Для анализа алгоритмов переменных направлений (а именно, мы рассматриваем АМЕп) удобно переписать их в следующем рекуррентном виде.

1. Построить Ai = X!¡iiAX¿i, bi = X^Jj.

2. Решить Ai«(1) = bu принять и = r(u(1),х(2),...

г®

3. Объединить х(1) := [и'1' г'1'],

.(2) .

, где г'1' является ТТ блоком

приближенной невязки, г = ..., яз 6 — Аи.

4. Перезапустить алгоритм для приближенного решения А^я'5*2' = 2 (рекурсия). Редуцированная матрица А¡¡2 и правая часть собираются так:

А>2 = Х^{2.....Л)АХФ{ 2.....6>2 = АГ;{2.....где

= = ж(1) ® 1п,...па е с(вд-

Третий шаг принципиально отличает алгоритм АМЕп и алгоритмы DMRG. Действительно, редуцированная система на шаге 4 является системой Галеркина в базисе столбцов матрицы .....л}: полное решение удовлетворяет

х = Х1х(>2), где « (Х1АХхУ1 {Х{Ъ) ,

т.е. приближенно выполняет условия Галеркина, Х{ (b — Ах) яа 0 (приближенность возникает из-за того, что система A>2xi>2> = 6S2 в свою очередь решается АМЕп методом, т.е. неточно).

Однако, заметим, что Xi = [u(1) z(1>] ® In2...ni = [СД Zj], где Zi = г(1) ® является фрейм-матрицей приближенной невязки, г = Z\z(1). Следовательно, Га-леркннский базис содержит вектор неточной невязки, z е span(Xi). Это дает возможность доказать глобальную сходимость метода АМЕп на основе теории сходимости неточного градиентного спуска.

Утверждение. Пусть дана SPD линейная система Ах = b и начальное приближение i, рассмотрим невязку z = 6 — At и вектор г, такой, что ¿.{z\ z) ^ в < 7г/2. Тогда неточный метод градиентного спуска

x = t + zh, h = arg min JAib(t + zh') = (23)

сходится, и скорость сходимости оценивается следующим образом (х* =

II®* - ®IU к -1 а , - ,,..l + sin6

= -¡1-^ гг—г = П < 1, к = cond(X)--—. (24)

-t||.4 к+1 1 — sinö 4 '

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

Теорема. АМЕп алгоритм является сходящимся, если ошибки аппроксимации невязок Zk f» Zk = Ь»t — вносимые на всех шагах рекурсии, удовлетворяют условию 9 = max/t=i.....d-1 ¿{Zk\h) < т/2. Скорость сходимости одной АМЕп

итерации ограничена сверху так, что имеет место неравенство

_<1-(1-П) , П = —, K = cond И)г—(25)

где t является начальным приближением, х - приближением после АМЕп итерации, а х„ = A~lb - точное решение.

В конце четвертой главы разбираются практические аспекты эффективной реализации алгоритма АМЕп, и приводится сравнение его с другом типом адаптивной одноблочной схемы DMRG, разработанным в сообществе квантовой физики. Показаны основные преимущества АМЕп: отсутствие возмущения решения в процессе внесения блока невязки, как следствие отсутствие эвристических параметров для масштабирования поправочных данных, и в итоге возможность доказательства глобальной сходимости. Приведена ссылка на публикацию автора [3], в которой также подтверждается и более высокая эффективность АМЕп в практических расчетах основных состояний спиновых цепей.

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

Заключение: основные результаты работы

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

• Центральным результатом диссертации является новый вычислительный алгоритм (АМЕп) для решения больших систем линейных уравнений и аппроксимации данных в форматах тензорных произведений. Этот метод был получен путем объединения БМИС схемы оптимизации в переменных направлениях, и классического метода градиентного спуска, и позволяет вычислять приближенные решения систем уравнений значительно быстрее и точнее, чем каждая из используемых схем в отдельности.

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

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

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

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

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

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

Публикации автора по теме диссертации

[1] Computation of extreme eigenvalues in higher dimensions using block tensor train format / S. V. Dolgov, B. N. Khoromskij, I. V. Oseledets, D. V. Savostyanov // Computer Phys. Comm. — 2014. — Vol. 185, no. 4. — P. 1207-1216.

[2] Dolgov S. V., Smirnov A. P., Tyrtyshnikov E. E. Low-rank approximation in the numerical modeling of the Farley-Buneman instability in ionospheric plasma // J. Сотр. Phys. - 2014. - Vol. 263. - P. 268-282.

[3] Dolgov S. V., Savostyanov D. V. Corrected one-site density matrix renormalization group and alternating minimal energy algorithm // Proc. ENUMATH 2013, accepted.- 2014,— URL: http://axxiv.org/abs/1312. 6542.

[4] Dolgov S., Khoromskij B. Two-Level QTT-Tucker Format for Optimized Tensor Calculus // SIAM J. on Matrix An. Appl. — 2013. — Vol. 34, no. 2. — P. 593-623.

[5] Dolgov S. V. TT-GMRES: solution to a linear system in the structured tensor format // Russ. J. Numer. Anal. Math. Modelling. — 2013. — Vol. 28, no. 2. — P. 149-172.

[6] Dolgov S. V., Khoromskij Boris N., Oseledets Ivan V. Fast solution of multidimensional parabolic problems in the tensor train/quantized tensor train-format with initial application to the Fokker-Planck equation // SIAM J. Sci. Comput. — 2012. - Vol. 34, no. 6. - P. A3016-A3038.

[7] Dolgov S. V., Oseledets I. V. Solution of linear systems and matrix inversion in the TT-format // SIAM J. Sci. Comput. - 2012. - Vol. 34, no. 5. - P. A2718-A2739.

Заказ № 35-Р/09/2014 Подписано в печать 08.09.14 Тираж 100 экз. Усл. пл. 1,0

ООО "Цифровичок", г. Москва, Большой Чудов пер., д.5

о, тел. (495)649-83-30

1/) www. cfr. ru; e-mail: :akpark@cfr. ru