Введение

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

В качестве примера приведем функцию f(x,y,z,u,v,w) определенную в области и соответствующую плотности распределения Больцмана. Для сетки, содержащей 100 узлов по каждой координате в шестимерном случае необходимо узлов памяти, содержащих числа. Для 64-значного числа это соответствует оперативной памяти 64 ТБ, что соответствует параметрам мощнейших современных суперкомпьютеров. Это значительно осложняет практические приложения, использующие уравнения Больцмана, как с точки стоимости расчетов, так и оперативности их применения.

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

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

Аппроксимация функций в многомерном пространстве представляет собой сложную задачу из-за экспоненциального роста требуемой оперативной памяти при увеличении размерности задачи (“проклятия размерности”). Одним из путей преодоления этих трудностей является использование тензорной формы многомерных задач и их аппроксимация с помощью тензорных разложений [1,2,3]. Этому вопросу и посвящена данная работа.

1. Каноническое разложение

Каноническое разложение [1,2] является самым фундаментальным тензорным разложением (на его основе определяется ранг тензора). Здесь - внешнее (тензорное) произведение векторов, в индексной форме имеющее вид .

В индексном виде каноническое разложение записывается как

.

(1)

Это выражение единственно с точностью до перестановки членов и масштабирования [4,5].

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

Задача определения набора ядер в вариационной форме имеет вид

.

(2)

Ключевой величиной при использовании канонического разложения является ранг тензора . По данным [1,4] она не вычислима из-за некорректности постановки задачи

.

(3)

В соответствии с [1, 4] использование канонического разложения страдает от неустойчивостей и нуждается в регуляризации.

2. Методы расчета тензорных разложений

Методы расчета тензорных разложений можно с некоторой условностью разделить на два подкласса: методы, основанные на матризации тензора (линейной алгебре), например [3] и вариационные методы [5] (которые делятся на прямые и итерационные).

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

Вариационные формулировки, как правило, опираются на метод наименьших квадратов по переменным направлениям (alternating least squares (ALS)) [6,7], и, как правило, тоже используют матризацию тензоров. Как правило, для канонического разложения ALS реализуется с помощью произведения Хатри-Рао (Khatri Rao product, ⊙) [6,7]). К сожалению, матризация тензора не приводит к его сжатию и не помогает против проклятия размерности.

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

3. Градиентная оптимизация на случайном зонтике

Формально качество аппроксимации функции каноническим разложением на каждом шаге глобальной итерации может быть оценено с помощью следующей невязки (записанной для 6D)

.

(4)

Здесь используются все узлы тензора, как в ALS с использованием произведения Хатри-Рао.

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

,

(5)

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

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

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

 

Рис. 1. Ансамбль точек, на котором производится оптимизация (“зонтик”)

 

Ядра получаются последовательным перебором.

Рассмотрим алгоритм на примере ядра разложения по в точке : и соответствующую невязку , (остальные элементы ядра определяются последовательно по )

.

(6)

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

Здесь точное значение функции в точке .

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

Возмутим это ядро на величину . Соответствующее возмущение невязки имеет вид:

.

(7)

Соответственно градиент невязки имеет вид:

.

(8)

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

.

(9)

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

Таким образом, вместо одной глобальной невязки (используемой в расчетах для контроля)

(10)

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

.

(11)

Компоненты ядер ищутся последовательной оптимизацией невязок по номеру ядра и номеру узла на ядре

(12)

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

.

(13)

Общая структура алгоритма имеет вид (цикл по ядрам и их узлам):

 

Алгоритм 1. Поиск ядер канонического разложения с помощью наискорейшего спуска

Input Initialize for (X,Y,Z,U,V,W) do ! цикл по координатам for i = do ! цикл по узлам вдоль ядра Поиск минимума реализуется методом наискорейшего спуска for do ! цикл по индексу ранга end for output

 

Такую примитивную структуру алгоритма можно реализовать из-за того, что градиенты приближенно разложимы:

.

(14)

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

3.1 Скользящая регуляризация

В работе использована скользящая по индексу ранга регуляризация, использующая Тихоновскую регуляризацию нулевого порядка [9], с коэффициентом регуляризации, увеличивающимся по индексу ранга в виде

(15)

Здесь - базовый коэффициент регуляризации.

Соответствующие итерации имеют вид:

(16)

Общая оптимизация сводится к локальной последовательной оптимизации по невязкам. Для контроля сходимости отслеживаются величины (5) и (10).

Оптимизация элементов ядер заканчивалась при , . Процесс построения канонического разложения заканчивался при , .

4. Результаты численных экспериментов по интерполяции многомерных функций

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

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

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

Результаты расчетов дают достаточно устойчивые и воспроизводимые оценки.

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

(17)

Это предельно простая функция, но ее априорный ранг не известен. Использовано 30 итераций, ранг вариировался от одного до 10. Зависимость величины невязки (5) от величины ранга представлена в Табл. 1.

Таблица 1. Зависимость невязки (5) от ранга для функции (17).

ранг	1	2	3	4	5	10
невязка

Характерно быстрое падение невязки при достижении истинного значения ранга (здесь- 2), и медленный рост при дальнейшем увеличении ранга. Из этих результатов видно, что, в целом, величина невязки может служить индикатором истинного ранга тензора. По мере увеличения ранга увеличивается зашумленность результата, по всей видимости, это отражает неустойчивости, возникающие вследствие некорректности для канонического разложения [1,4].

На Рис. 2 и 3 представлена точная функция (17) и ее аппроксимация при ранге 2.

 

Рис. 2. Точная функция	Рис. 3. Аппроксимация с помощью канонического разложения

 

Функция (17) является двумерной, рассмотрим более сложный истинно шестимерный случай, задаваемый гауссианой и суммой синусов по всем направлениям в виде

.

(18)

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

 

Рис. 4. Точная функция (18)	Рис. 5. Аппроксимация (18) с помощью канонического разложения

 

На Рис. 6 представлено поведение разных критериев сходимости в зависимости от числа итераций. В логарифмической форме представлены невязки (10) (Eps), глобальная невязка, оцененная с помощью Монте-Карло (5) (Eps_MC) и норма градиента невязки (grad norm).

 

Рис. 6. Разные критерии сходимости в зависимости от числа итераций.

 

Невязки (10) (Eps) и глобальная невязка, оцененная с помощью Монте-Карло (5) (Eps_MC), демонстрируют близкое поведение, пригодное для контроля сходимости. Норма градиента невязки (grad norm) на начальных итерациях ведет себя немонотонно, что ограничивает ее применимость в качестве критерия останова.

Строго говоря, используемый метод не предусматривает монотонной минимизации глобальных критериев качества, так как минимизация осуществляется поточечно. Однако, численные эксперименты демонстрируют монотонную сходимость глобальных критериев. Данный алгоритм обеспечивает заметно более быструю сходимость в сравнении с алгоритмом, использованным в [10].

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

 

Рис. 7. Зависимость невязки от величины ранга для функции (18)

5. Результаты численных экспериментов по аппроксимации результатов параметрических расчетов

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

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

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

Используя структуру данных типа “зонтик” будем строить каноническое разложение. Один наш расчет поля течения (элементарный “строительный блок”) - это полностью заполненный тензор порядка 3.

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

На Рис. 8 и 9 представлены результаты расчета поля течения (плотности) для (отклонение верхнего потока), (отклонение нижнего потока), и соответствующей аппроксимации с помощью канонического разложения (ранг 15, число точек ансамбля-2000). При этом использован тензор порядка три . Направление течения справа, на переднем плане – зона повышенной плотности за пересекающимися скачками уплотнения.

 

Рис. 8. Расчет поля плотности	Рис. 9. Аппроксимация поля плотности с помощью канонического разложения (тензор порядка три)

 

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

На Рис. 10 и 11 представлены результаты расчета поля течения (плотности) (, , ) и соответствующей аппроксимации с помощью канонического разложения (ранг 3, число точек ансамбля-200). При этом использован тензор порядка четыре: , в котором использовано пять слоев по :

 

Рис. 10. Расчет поля плотности	Рис. 11. Аппроксимация поля плотности с помощью канонического разложения (тензор порядка четыре)

 

Результаты, представленные на Рис. 10 и 11 демонстрируют работоспособность аппроксимации (сжатия) однопараметрического поля течения с помощью канонического разложения. Наблюдаемые осцилляции аппроксимации вероятно связаны с некорректностью задачи определения канонического разложения. Есть определенная надежда, что при переходе к тензорному поезду качество решения улучшится.

6. Обсуждение

В связи с некорректностью задачи [4] при определении ядер канонического разложения предусмотрена возможность квадратичной регуляризации нулевого порядка по Тихонову [9] со скользящей регуляризацией, позволяющей сжимать аппроксимацию ранга , обеспечивающую невязку (получать оценку ранга сверху).

Тензорные разложения достаточно активно используются для целей визуализации [11], поскольку позволяют построить легко вычислимую модель, аппроксимирующую труднообрабатываемый набор данных в параметрическом пространстве. Тензорные разложения (каноническое разложение, тензорный поезд, иерархический Таккер) используются для экономичного решения многомерных задач типа уравнения Больцмана [5,13,14]. Часто для этих целей используется формат тензорного поезда и кросс-аппроксимация [13]. Рассматриваемый здесь подход легко расширяется на формат тензорного поезда и лишен ограничений на вычислительные ресурсы, характерных для кросс-аппроксимации (связанных с матризацией тензора). Используемое здесь каноническое разложение позволяет эффективно аппроксимировать и хранить многомерные функции. Затраты компьютерного времени на работу с функциями в шестимерном пространстве (при использовании 100 узлов по каждой координате, что формально требует хранение и работу с чисел) составляют 1-3 минуты на ПК (процессор Intel I5, 2.66 ГГц) при затратах памяти на хранение ядер в максимальном случае около чисел. Есть надежда, что переход к формату тензорного поезда (планируемый в дальнейших работах) позволит преодолеть часть недостатков канонического разложения, связанных с некорректностью и неустойчивостью.

Заключение

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

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

Библиографический список

1. W. Hackbusch. Tensor Spaces and Numerical Tensor Calculus. Springer, 2012.

2. H. Yorick, S. Willi-Hans, Matrix Calculus, Kronecker Product and Tensor Product: A Practical Approach to Linear Algebra, Multilinear Algebra and Tensor Calculus with Software I, WOS 2019

3. Oseledets I. V., Tensor- train decomposition, SIAM J. Sci. Comput., 33 (2011), pp. 2295–2317

4. V. D. Silva and L -H. Lim, Tensor rank and the ill-posedness of the best low rank approximation problem, SIAM J. Matrix Anal. Appl., 30(3) 1084-1127, 2008.

5. A. M. P. Boelens, D. Venturi, D. M. Tartakovsky, Parallel tensor methods for high-dimensional linear PDEs, J. Computat. Phys. 375 (2018) 519-539.

6. Shuangzhe Liu, Gaotz Trenkler, Hadamard, Khatri-Rao, Kronecker and other matrix products, Int. J. Information and system Scienses, 2008, Volume 4, Number 1, Pages 160-177

7. P. Comon, X. Luciani, and A.L.F. De Almeida, Tensor decompositions, alternating least squares and other tales, Journal of Chemometrics, vol. 23, no. 7-8, pp. 393–405, 2009. A. Uschmajew, Local convergence of the alternating least squares algorithm for canonical tensor approximation, SIAM J. Matrix Anal. Appl. 33 (2) (2012) 639–652,

8. S. Rabanser, O. Shchur, S. Gunnemann, Introduction to Tensor Decompositions and their Applications in Machine Learning, arXiv:1711.10781v1, 2017

9. Тихонов А.Н., Арсенин В.Я., Методы решения некорректных задач, М., Наука, 1979.

10. Алексеев А.К., Бондарев А.Е., Пятакова Ю.С., О визуализации функций в многомерном пространстве с помощью канонического разложения, Научная визуализация, т. 14, N 3, c. 73-91, 2022

11. Yannis Panagakis, Jean Kossaifi, Grigorios G. Chrysos, James Oldfield, Mihalis A. Nicolaou, Anima Anandkumar and Stefanos Zafeiriou, Tensor Methods in Computer Vision and Deep Learning, arXiv:2107.03436v1, 2021

12. Oseledets I., Tyrtyshnikov E., TT-cross approximation for multidimensional arrays, Linear Algebra Appl., 432 (2010), pp. 70–88.

13. Arnout M. P. Boelens, Daniele Venturi, Daniel M. Tartakovsky, Tensor methods for the Boltzmann-BGK equation, arXiv:1911.04904v2 2020

14. A.V. Chikitkin, E.K. Kornev, V.A. Titarev, Numerical solution of the Boltzmann equation with S-model collision integral using tensor decompositions, arXiv:1912.04582v1 2019

---

On Application of Canonical Decomposition for the Visualization of Results of Multiparameter Computations

Authors: A.K.  Alekseev^1,A,B, A.E.  Bondarev^2,A, Yu. S.  Pyatakova^3,B

^A Keldysh Institute of Applied Mathematics RAS

^B RSC Energia, Korolev, Russia

¹ ORCID: 0000-0001-8317-8688, aleksey.k.alekseev@gmail.com

² ORCID: 0000-0003-3681-5212, bond@keldysh.ru

³ ORCID: 0000-0002-8055-7807 , yuliya.pyatakova@rsce.ru

 

Abstract

The approximation of the tensor appearing at a discretization of the multidimensional function is considered from the viewpoint of storing and treating of the results of parametric computations obtained in computational aerogasdynamics. The new algorithm for the computation of the canonical decomposition using gradient descent and approximately decomposable goal functional is described.

This algorithm applies the random set of points on the hyperplane orthogonal to the computed core of the canonical decomposition (“umbrella”) that ensures its flexible application for an approximation of the tensors with a priori unknown rank and may be naturally transferred on such tensor decomposition as the tensor train. The results of the numerical tests are presented for the model six-dimensional functions and for an ensemble of the numerical solutions for the two-dimensional Euler equations. These equations describe the flow of the compressible gas with two crossing shock waves. The Mach number and angles of the flow deflection serve as the flow parameters. The results are provided for the dimensionality 3 (simple numerical solution) and 4 (the ensemble of the numerical solutions in dependence on the Mach number).

 

Keywords: canonical decomposition, parametric calculations, computational aerogas dynamics, gradient optimization, sliding regularization.

 

References

1. W. Hackbusch. Tensor Spaces and Numerical Tensor Calculus. Springer, 2012.

2. H. Yorick, S. Willi-Hans, Matrix Calculus, Kronecker Product and Tensor Product: A Practical Approach to Linear Algebra, Multilinear Algebra and Tensor Calculus with Software I, WOS 2019

3. Oseledets I. V., Tensor- train decomposition, SIAM J. Sci. Comput., 33 (2011), pp. 2295–2317

4. V. D. Silva and L -H. Lim, Tensor rank and the ill-posedness of the best low rank approximation problem, SIAM J. Matrix Anal. Appl., 30(3) 1084-1127, 2008.

5. A. M. P. Boelens, D. Venturi, D. M. Tartakovsky, Parallel tensor methods for high-dimensional linear PDEs, J. Computat. Phys. 375 (2018) 519-539.

6. Shuangzhe Liu, Gaotz Trenkler, Hadamard, Khatri-Rao, Kronecker and other matrix products, Int. J. Information and system Scienses, 2008, Volume 4, Number 1, Pages 160-177

7. P. Comon, X. Luciani, and A.L.F. De Almeida, Tensor decompositions, alternating least squares and other tales, Journal of Chemometrics, vol. 23, no. 7-8, pp. 393–405, 2009. A. Uschmajew, Local convergence of the alternating least squares algorithm for canonical tensor approximation, SIAM J. Matrix Anal. Appl. 33 (2) (2012) 639–652,

8. S. Rabanser, O. Shchur, S. Gunnemann, Introduction to Tensor Decompositions and their Applications in Machine Learning, arXiv:1711.10781v1, 2017

9. Tikhonov, A.N., Arsenin, V.Y.: Solutions of Ill-Posed Problems. Winston and Sons, Washington DC (1977).

10. A.K. Alekseev, A.E. Bondarev, Yu.S. Pyatakova, On the visualization of multidimensional functions using canonical decomposition, Scientific Visualization, 2022, volume 14, number 3, pages 73 - 91

11. Yannis Panagakis, Jean Kossaifi, Grigorios G. Chrysos, James Oldfield, Mihalis A. Nicolaou, Anima Anandkumar and Stefanos Zafeiriou, Tensor Methods in Computer Vision and Deep Learning, arXiv:2107.03436v1, 2021

12. Oseledets I., Tyrtyshnikov E., TT-cross approximation for multidimensional arrays, Linear Algebra Appl., 432 (2010), pp. 70–88.

13. Arnout M. P. Boelens, Daniele Venturi, Daniel M. Tartakovsky, Tensor methods for the Boltzmann-BGK equation, arXiv:1911.04904v2 2020

14. A.V. Chikitkin, E.K. Kornev, V.A. Titarev, Numerical solution of the Boltzmann equation with S-model collision integral using tensor decompositions, arXiv:1912.04582v1 2019