Tensors in modelling multi-particle interactions

TL;DR

This paper demonstrates how low-rank tensor decompositions can efficiently model multi-particle interactions in aggregation kinetics, significantly reducing computational complexity and enabling high-accuracy solutions on supercomputers.

Contribution

The work introduces a novel application of tensor decompositions to accelerate the evaluation of nonlinear differential equations in multi-particle kinetics.

Findings

Reduced computational complexity using tensor decompositions.

High-accuracy numerical solutions achieved with modest computational times.

Preliminary parallel scalability results indicate efficiency on supercomputers.

Abstract

In this work we present recent results on application of low-rank tensor decompositions to modelling of aggregation kinetics taking into account multi-particle collisions (for three and more particles). Such kinetics can be described by system of nonlinear differential equations with right-hand side requiring $N^D$ operations for its straight-forward evaluation, where $N$ is number of particles size classes and $D$ is number of particles colliding simultaneously. Such a complexity can be significantly reduced by application low rank tensor decompositions (either Tensor Train or Canonical Polyadic) to acceleration of evaluation of sums and convolutions from right-hand side. Basing on this drastic reduction of complexity for evaluation of right-hand side we further utilize standard second order Runge-Kutta time integration scheme and demonstrate that our approach allows to obtain numerical solutions of studied equations with very high accuracy in modest times. We also show preliminary results on parallel scalability of novel approach and conclude that it can be efficiently utilized with use of supercomputers.

Tables

Table 1. Table 1: Speedup of computations with use of ClusterFFT operation for pure ternary aggregation in case of 3-dimensional generalized Brownian kernel with N = 2 19 𝑁 superscript 2 19 N=2^{19} equations. Benchmark for 100 time-integration steps with use of second order Runge-Kutta method. “Zhores” supercomputer of Skolkovo Institute of Science and Technology tables.

Number of CPU-cores	time, sec	Speedup
1	257.80	1.00
2	147.62	1.75
4	80.21	3.21
8	43.65	5.91
16	22.63	11.39
32	14.83	17.38
64	13.15	19.60
128	12.22	21.09

Equations

$$\frac{d{n}_k}{dt}=\frac{1}{2}\sum _{i+j=k}{C}_{i,j}{n}_i{n}_j-n_k\sum _{i=1}^{\infty }{C}_{k,i}{n}_k.$$

$$\frac{\mathrm{d}\mathbf{n}}{\mathrm{d}t}=\sum _{d=2}^D{\mathcal{S}}^{(D)}\left[\mathbf{n}\right]=\sum _{d=2}^D\left[{\mathcal{P}}^{(D)}\left[\mathbf{n}\right]+{\mathcal{Q}}^{(D)}\left[\mathbf{n}\right]\right],$$

$$p_k^{(d)}\left[\mathbf{n}\right]=\frac{1}{d!}\sum _{|{\mathbf{i}}_d|=k}{C}_{{\mathbf{i}}_d}^{(d)}{n}_{i_1}{n}_{i_2}\dots n_{i_d},k\in \mathbb{N},$$

$$q_k^{(d)}\left[\mathbf{n}\right]=-\frac{n_k}{\left(d-1\right)!}\sum _{{\mathbf{i}}_{d-1}\in {\mathbb{N}}^{d-1}}{C}_{{\mathbf{i}}_{d-1},k}^{(d)}{n}_{i_1}{n}_{i_2}\dots n_{i_{d-1}},k\in \mathbb{N},$$

$${\mathbf{i}}_d=\left(i_1,i_2,\dots ,i_d\right),\left|{\mathbf{i}}_d\right|=i_1+i_2+\dots +i_d,2\le d\le D.$$

$$C_{{\mathbf{i}}_D}^{(D)}\equiv C_{{\mathbf{i}}_D}^{(D)}\left[{\mu }_1,{\mu }_2,\dots ,{\mu }_D\right]=\sum _{\sigma }{i}_{\sigma \left(1\right)}^{{\mu }_1}\cdot i_{\sigma \left(2\right)}^{{\mu }_2}\cdot \dots \cdot i_{\sigma \left(D\right)}^{{\mu }_D}$$

$$C_{{\mathbf{i}}_D}^{(D)}=\sum _{r_0,r_1,\dots ,r_D}{H}_{r_0,i_1,r_1}^{(1)}\cdot H_{r_1,i_2,r_2}^{(2)}\cdot \dots \cdot H_{r_{D-1},i_D,r_D}^{(D)},$$

$$1\le r_{\lambda }\le R_{\lambda },0\le \lambda \le D;R_0=R_D=1.$$

$$\max \left\{R_{\lambda }:0\le \lambda \le D\right\}=\left(\genfrac{}{}{0ex}{}{D}{\lceil D/2\rceil }\right)\equiv O\left(2^D/\sqrt{D}\right).$$

$$R_{\lambda }\le \left(\genfrac{}{}{0ex}{}{D}{\lambda }\right)\equiv \frac{D!}{\lambda !\cdot \left(D-\lambda \right)!},0\le \lambda \le D.$$

$$\left\{r_{\lambda }\in \mathbb{N}:1\le r_{\lambda }\le R_{\lambda }\right\},$$

$${\mathcal{R}}_{\lambda }:=\left\{\left(r_{1,\lambda },r_{2,\lambda },\dots ,r_{\lambda ,\lambda }\right)\in {\mathbb{N}}^{\lambda }:1\le r_{1,\lambda }<r_{2,\lambda }<\dots <r_{\lambda ,\lambda }\le D\right\}.$$

$$C_{{\mathbf{i}}_{\lambda +1}}^{(\lambda +1)}\left[{\mu }_1,{\mu }_2,\dots ,{\mu }_{\lambda +1}\right]\equiv \sum _{\xi =1}^{\lambda +1}{C}_{{\mathbf{i}}_{\lambda }}^{(\lambda )}\left[{\mu }_1,\dots ,{\mu }_{\xi -1},{\mu }_{\xi +1},\dots ,{\mu }_{\lambda +1}\right]\cdot i_{\lambda +1}^{{\mu }_{\xi }}.$$

$$\sum _{r_0,r_1,\dots ,r_{\lambda -1}}{H}_{r_0,i_1,r_1}^{(1)}\cdot H_{r_1,i_2,r_2}^{(2)}\cdot \dots \cdot H_{r_{\lambda -1},i_{\lambda },r_{\lambda }}^{(\lambda )}=C_{{\mathbf{i}}_{\lambda }}^{(\lambda )}\left[{\mu }_{r_{1,\lambda }},{\mu }_{r_{2,\lambda }},\dots ,{\mu }_{r_{\lambda ,\lambda }}\right],$$

$$r_0=1,1\le i_{\tau }\le N,1\le r_{\tau }\le R_{\tau },1\le \tau \le \lambda .$$

$$H_{r_0,i_1,r_1}^{(1)}:=C_{i_1}^{(1)}\left[{\mu }_{r_1}\right]\equiv i_1^{{\mu }_{r_1}},1\le i_1\le N,1\le r_1\le R_1.$$

$$C_{i_1,i_2}^{(2)}\left[{\mu }_1,{\mu }_2\right]\equiv \left[\begin{array}{cc}{i}_1^{{\mu }_1}& i_1^{{\mu }_2}\end{array}\right]\cdot \left[\begin{array}{c}{i}_2^{{\mu }_2}\\ i_2^{{\mu }_1}\end{array}\right];$$

$$C_{i_1,i_2,i_3}^{(3)}\left[{\mu }_1,{\mu }_2,{\mu }_3\right]\equiv \left[\begin{array}{ccc}{i}_1^{{\mu }_1}& i_1^{{\mu }_2}& i_1^{{\mu }_3}\end{array}\right]\cdot \left[\begin{array}{ccc}0& i_2^{{\mu }_3}& i_2^{{\mu }_2}\\ i_2^{{\mu }_3}& 0& i_2^{{\mu }_1}\\ i_2^{{\mu }_2}& i_2^{{\mu }_1}& 0\end{array}\right]\cdot \left[\begin{array}{c}{i}_3^{{\mu }_1}\\ i_3^{{\mu }_2}\\ i_3^{{\mu }_3}\end{array}\right].$$

$$H_{r_{\lambda },i_{\lambda +1},r_{\lambda +1}}^{(\lambda +1)}:=\{\begin{array}{ll}{i}_{\lambda +1}^{{\mu }_{r_{\xi ,\lambda +1}}},& \left\{r_{\xi ,\lambda +1}\right\}\cup {\mathcal{R}}_{\lambda }={\mathcal{R}}_{\lambda +1},\\ 0,& \text{otherwise};\end{array}$$

$$1\le r_{\lambda }\le R_{\lambda },1\le i_{\lambda +1}\le N,1\le r_{\lambda +1}\le R_{\lambda +1}.$$

$$\sum _{r_{\lambda }=1}^{R_{\lambda }}{C}_{{\mathbf{i}}_{\lambda }}^{(\lambda )}\left[{\mu }_{r_{1,\lambda }},{\mu }_{r_{2,\lambda }},\dots ,{\mu }_{r_{\lambda ,\lambda }}\right]\cdot H_{r_{\lambda },i_{\lambda +1},r_{\lambda +1}}^{(\lambda +1)}=$$

$$=\sum _{\xi =1}^{\lambda +1}{C}_{{\mathbf{i}}_{\lambda }}^{(\lambda )}\left[{\mu }_{r_{1,\lambda +1}},\dots ,{\mu }_{r_{\xi -1,\lambda +1}},{\mu }_{r_{\xi +1,\lambda +1}},\dots ,{\mu }_{r_{\lambda +1,\lambda +1}}\right]\cdot i_{\lambda +1}^{{\mu }_{r_{\xi ,\lambda +1}}}.$$

$${\left\{{\mathbf{a}}_N\right\}}_p:={\left[a_{\left(p-1\right)N/P+1},a_{\left(p-1\right)N/P+2},\dots ,a_{pN/P}\right]}^T,1⩽p⩽P.$$

$${\mathbf{h}}_{r_{\lambda -1},r_{\lambda }}^{(\lambda )}:={\left[H_{r_{\lambda -1},1,r_{\lambda }}^{(\lambda )},H_{r_{\lambda -1},2,r_{\lambda }}^{(\lambda )},\dots ,H_{r_{\lambda -1},N,r_{\lambda }}^{(\lambda )}\right]}^T\in {\mathbb{R}}^N,$$

$${\mathbf{H}}_{r_{\lambda }}^{(\lambda ,p)}:={\left[{\left\{{\mathbf{h}}_{1,r_{\lambda }}^{(\lambda )}\right\}}_p,{\left\{{\mathbf{h}}_{2,r_{\lambda }}^{(\lambda )}\right\}}_p,\dots ,{\left\{{\mathbf{h}}_{R_{\lambda -1},r_{\lambda }}^{(\lambda )}\right\}}_p\right]}^T\in {\mathbb{R}}^{R_{\lambda -1}\times N/P},$$

$$\displaystyle\mathscr{H}^{(\lambda,\hskip 0.85358ptp)}:=\left[\left\{\mathbf{h}^{(\lambda)}_{r_{\lambda-1},\hskip 0.85358ptr_{\lambda}}\right\}_{\!p}\right]\equiv\Bigl{[}\mathbf{H}^{(\lambda,\hskip 0.85358ptp)}_{\hskip 0.85358ptr_{\lambda}}\Bigr{]}\in\mathbb{R}^{R_{\lambda-1}\times N/P\times R_{\lambda}}.$$

Full text

1 Introduction

Aggregation of inelastically colliding particles plays important role in many technological and natural phenomena. In case of spatially homogeneous systems aggregation process can be described by famous Smoluchowski kinetic equations [1]. These equations describe time-evolution of mean concentrations $n_{k}$ of particles of size $k$ per unit volume of media:

[TABLE]

Such a model is well-studied by lots of analytical and numerical methods but allows to take into account only pairwise particles’ collisions. In this work we consider a generalization of aggregation equations for case of multi-particle interactions

[TABLE]

where operators  $\mathcal{P}^{(d\hskip 0.85358pt)}=\left[\hskip 0.85358ptp^{(d\hskip 0.85358pt)}_{1}\!,\hskip 0.85358ptp^{(d\hskip 0.85358pt)}_{2}\!,\hskip 2.84526pt\ldots\hskip 2.84526pt\right]^{T}$ and  $\mathcal{Q}^{(d\hskip 0.85358pt)}=\left[q^{(d\hskip 0.85358pt)}_{1}\!,\hskip 0.85358ptd^{(d\hskip 0.85358pt)}_{2}\!,\hskip 2.84526pt\ldots\hskip 2.84526pt\right]^{T}$ are defined as

[TABLE]

Those operators $p^{(d\hskip 0.85358pt)}$ and $q^{(d\hskip 0.85358pt)}$ correspond to description of simultaneous collisional aggregation of $d$ particles. We represent such a model informally in Fig. 1.

In case of defined initial conditions $n_{k}(t=0)$ we obtain a Cauchy problem which is known to be well-posed under assumption of bounded non-negative symmetric coefficients [2]. There also exist a very few examples of known analytical solutions for the Cauchy problem for multi-particle aggregation kinetic equations [3]. Unfortunately, such class of mathematical models is much less studied numerically than class of aggregation equations accounting only binary collisions. The reason lies in higher order of non-linearity and exponential growth of complexity of evaluation of the right-hand side with respect to number $D$ of simultaneously colliding particles. It makes numerical treatment of such systems extremely time-consuming.

In our recent work we proposed a novel approach based on application of low-rank tensor train (TT) decomposition  [5, 6] to acceleration of computations and presented its high accuracy  [7, 8]. This allowed us to reduce the complexity of evaluation of sums from the right-hand side from $O(N^{D})$ operations to $O(NDR^{2}\log N)$ operations, where $R\ll N$ is the maximal TT-rank of kinetic coefficients. In case of pre-defined low-rank Canonical Polyadic (CP) decomposition of kernel coefficients the complexity can be reduced even to $O(NDR\log N)$ but we do not know a robust way finding CP-decomposition even for 3-dimensional tensors. Nevertheless, assumption that $R\ll N$ is crucial for efficiency of proposed approach.

In the current work we prove that TT-ranks of a wide class of generalized Brownian kernels are low and do not depend on number $N$ of accounted kinetic equations for arbitrary dimension $D$. We also present an efficient parallel implementation of our TT-based approach and present preliminary tests of scalability our method.

2 Estimates of TT-ranks for generalized Brownian kernels

In this section we present estimates of tensor ranks for generalized Brownian kernel $\mathscr{C}^{(D)}=\left[C^{(D)}_{\mathbf{i}_{D}}\right]$ (see example of exact Brownian coefficients e.g. at [4]) with the elements of the following form:

[TABLE]

Here we assume the sum to be over all permutations $\sigma$ of the set $\left\{1,2,\ldots,D\right\}$.

Recall the definition of a TT-decomposition for a kernel $\mathscr{C}^{(D)}$ which is of the form

[TABLE]

For the kernel $\mathscr{C}^{(D)}$ of dimension $D$ with the elements (1) there holds an estimate on the TT-ranks $R_{\lambda}$:

[TABLE]

Such estimate can be verified by the means of the following

Theorem 2.1

Let the parameters $\mu_{1},\mu_{2},\ldots,\mu_{D}$ be fixed. For a given tensor $\mathscr{C}^{(D)}$ in $D$ dimensions of sizes $N\times N\times\ldots\times N$ with the elements (1) one can prove the following estimates on its TT-ranks $R_{\lambda}$:

[TABLE]

Proof

Let us put $R_{\lambda}:=\dbinom{D}{\lambda}$ for $0\leq\lambda\leq D$. To prove the theorem we simply need to construct a tensor train decomposition of the tensor $\mathscr{C}^{(D)}$ with these predefined ranks $R_{\lambda}$.

Further we assume that for each number $1\leq\lambda\leq D$ it is chosen a bijection $r_{\lambda}\to\left(r_{1,\hskip 0.85358pt\lambda},\hskip 0.85358ptr_{2,\hskip 0.85358pt\lambda},\,\ldots,\,r_{\lambda,\hskip 0.85358pt\lambda}\right)$ between the sets

[TABLE]

One can specify such mappings for sure due to the coincidence of the cardinalities of the considered sets.

Note, that for all numbers $1\leq\lambda\leq D-1$ one can check the correctness of the identity

[TABLE]

With the given identity we show by induction on $1\leq\lambda\leq D$, that it is always possible to choose the values $H^{(\tau)}_{r_{\tau-1},\hskip 0.85358pti_{\tau},\hskip 0.85358ptr_{\tau}}$ in order to satisfy the constraints

[TABLE]

Thus, the equality (4) with $\lambda=D$ gives us the required TT-decomposition for the tensor with the elements of the form (1).

The base of induction is trivial, if choose

[TABLE]

Next we rewrite the equation (1) for $D=2,3$:

[TABLE]

This representation allows us to describe the structure of factors $\mathscr{H}^{(\tau)}$ for each $\tau$. If the equality (4) is already proven for a given $\lambda$, then it is sufficient to choose

[TABLE]

Now it is not hard to prove, that by the virtue of the proposed choice for all numbers $1\leq i_{1},\hskip 0.85358pti_{2},\hskip 0.85358pt\ldots,\hskip 0.85358pti_{\lambda+1}\leq N$ and $1\leq r_{\lambda+1}\leq R_{\lambda+1}$ we have the identity

[TABLE]

Therefore, to verify an induction step we just need to use the identity (3), where it is necessary to use parameters $\mu_{\hskip 0.85358ptr_{1,\hskip 0.28453pt\lambda+1}},\hskip 0.28453pt\mu_{\hskip 0.85358ptr_{2,\hskip 0.28453pt\lambda+1}},\,\ldots,\,\mu_{\hskip 0.85358ptr_{\lambda+1,\hskip 0.28453pt\lambda+1}}$ instead of parameters $\mu_{1},\hskip 0.85358pt\mu_{2},\hskip 0.85358pt\ldots,\hskip 0.85358pt\mu_{\lambda+1}$ respectively. The last statement proves the theorem.

3 Parallel algorithm and numerical experiments

In our work we exploit organization of parallel computations along particle size coordinate with dimension $N$. It is worth to note that alternative way of parallelization of our approach along TT-ranks leads to dramatic overheads in terms of data exchanges and collective operations and does not lead to speedup of computations.

Let us assume that we have $P$ processors and number $N$ of studied kinetic equations is divisible by $P$. Thus, we introduce the following notation for getting $p$-th block of taken vector  $\mathbf{a}_{N}=\left[sa_{1},\hskip 0.85358pta_{2},\hskip 0.85358pt\ldots,\hskip 0.85358pta_{N}\right]^{T}$:

[TABLE]

With use of those notations we can denote blocks и $\mathscr{H}^{(\lambda,\hskip 0.85358ptp)}$ of cores $\mathscr{H}^{(\lambda)}$ of TT-decomposition of kernel kinetic coefficients  $\mathscr{C}^{(D)}$​ which will be used at processor with number  $p$ ( $1\leq p\leq P$) :

[TABLE]

The algorithm for operator $\mathcal{P}^{(d\hskip 0.85358pt)}$ is presented in Fig. 2 and for $\mathcal{Q}^{(d\hskip 0.85358pt)}$ in Fig. 3. On the input algorithms require to have blocks $\mathscr{H}^{(1,i)}$ of TT-decomposition for kinetic coefficients and vector of concentrations. For time-integration of the Cauchy problem we utilize standard explicit second order Runge-Kutta method, hence, each time-step requires two evaluations of $\mathcal{P}^{(d\hskip 0.85358pt)}$ and $\mathcal{Q}^{(d\hskip 0.85358pt)}$.

We present results of benchmarks of presented algorithm for the generalized Brownian coefficients in Table 1. In our experiments we used ClusterFFT library included into Intel MKL™. As soon as FFT is a dominating operation in our algorithm in terms of complexity, we obtain similar performance of our code to performance of ClusterFFT library. From these experiments we obtain acceleration of calculations by order of magnitude. This allows us to consider a broader class of problems of potential interest which can be studied in modest computational time.

4 Conclusions

In this paper we present recent developments of tensor based methods for modelling of multi-particle aggregation. We prove estimates of TT-ranks for generalized Brownian kinetic coefficients depending only on dimensionality $D$ but not mode-sizes $N$ of used arrays. We also propose an efficient way of parallel implementation of TT-based approach and demonstrate preliminary results of its parallel scalability.

In our work we used “Zhores” supercomputer installed at Skolkovo Institute of Science and Technology [9]. We also would like to acknowledge Talgat Daulbaev for an idea of representation of generalized Brownian coefficients in TT-format in case of $D=3$. The work was supported by the Russian Science Foundation, grant 19-11-00338.