{'fraction_non_alphanumeric': 0.033807231613119895, 'fraction_numerical': 0.014216948393596777, 'mean_word_length': 4.572988147223955, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Model reduction in Smoluchowski-type equations

In this paper we utilize the Proper Orthogonal Decomposition (POD) method for model order reduction in application to Smoluchowski aggregation equations with source and sink terms. In particular, we show in practice that there exists a low-dimensional space allowing to approximate the solutions of aggregation equations. We also demonstrate that it is possible to model the aggregation process with the complexity depending only on dimension of such a space but not on the original problem size. In addition, we propose a method for reconstruction of the necessary space without solving of the full evolutionary problem, which can lead to significant acceleration of computations, examples of which are also presented.


Introduction
A classical model of aggregation kinetics is based on the Smoluchowski equations, dating back to the original work by Marian von Smoluchowski [1]. In the original form, these equations describe an evolution of a spatially uniform system of agglomerates of different sizes, via an infinite system of ordinary differential equations for the concentrations n k of particles of size k each. The original formulation has been later amended by Hans Muller [2] to model continuous particle size distribution or additional phenomena, such as particle fragmentation [3] and others.
The range of phenomena modelled via Smoluchowski kinetic equations has also expanded over time, from molecular scales [3][4][5][6] to astronomical [7][8][9]. More detailed information about possible applications of aggregation-based models can be found in extensive reviews [10,11] and references therein.
Whether an original discrete system is used, or a discretization of the continuous, one still has to deal with a rather large systems of nonlinear differential equations, especially if particle masses differ by several orders of magnitude. Because of that, accurate numerical simulation of these systems is quite challenging. While there has been some recent progress on this front, bringing complexity for some classes of coagulation kernels down to almost linear [12,13] it is still insufficient for some of the larger systems arising in practice.
In this paper, we will attack this problem using the ideas of model reduction via Proper Orthogonal Decomposition, as outlined in [14]. Specifically, we are interested in the method of snapshots, introduced in [15]. The main idea of the method is to construct a low-dimensional vector space containing the solution or its approximation by examining its snapshots at different time moments. The end goal here is to create an opportunity to describe and approximate the solution using significantly fewer parameters than the full dimensionality of the system.
In this paper we demonstrate that • a low-dimensional space in which the solution can be approximated with reasonable accuracy exists; • once such a space is found, it is possible to model the system within the complexity depending only its dimension but not on the original problem size; • at least for some cases, it is possible to find the necessary space without constructing the solution of the full original problem.
Even though the results presented here do not seem to be immediately applicable for complex industrial applications, we believe that we suggest a novel concept for solving the aggregation-fragmentation equations leading to a fruitful and challenging avenue of further research. In some sense, our approach gives an alternative view at developing deep learning-based methods [16] for non-linear time-dependent problems with attractors and cycles. In contrast to [16] we deal with much larger systems of ODEs (tens of thousands in our work instead of dozens or hundreds).
The rest of the paper is organized as following: in Section 2 we discuss the target set of kinetic equations and recall the necessary facts about properties of the solution and the coefficients. In Section 3 we introduce a numerical method allowing to solve the target equations in approximate form using the reduction basis. The next Section 4 is devoted to algorithm allowing to construct such a basis via Proper Orthogonal Decomposition (namely, the method of snapshots). In Section 5 we demonstrate the results of numerical experiments and validation of the proposed methodology. In our experiments, we demonstrate the existence of the required low-dimensional reduced basis allowing one to accelerate the computations of numerical solutions of aggregation equations. In this Section, we also discuss the drawbacks of the proposed approach and further accumulate our findings in the conclusions of Section 6.

Problem setting
In our work, we consider the model similar to one originally posed in [1], with the addition of a constant source of particles [17,18]: In this system, n k stands for the concentration of particles of mass k; C i j is a coagulation kernel, characterising the frequencies of collisions between particles of size i and j; J k is a uniform source of particles of size k.
We additionally put some physically relevant constraints on these variables: n k ≥ 0 (there cannot be a negative concentrations of any kind of particles in the system); C i j = C j i ≥ 0 (the coagulation kernel is symmetric and non-negative); J k ≥ 0 (this term corresponds to the source of new particles).
For modelling purposes, we truncate (1) to get a finite system; this is equivalent to postulating an immediate removal of large particles from the system (see e.g. [18]): Given a sufficiently large N , system (2) approximates (1) with reasonable accuracy either in steady-state form [19] or quasi-steady-state [20]. For some cases of kernel coefficients with finite N a steady collective oscillatory solutions of aggregation equations [18] exists, which cannot be expected for the pure infinite aggregation system with source but no sink. However, the required value of N in practice can still be fairly large, so our aim for the rest of the paper is to reduce the number of parameters in (2).

Model reduction
In order to reduce the number of variables in the system (2), we employ the model reduction concept via Proper Orthogonal Decomposition (POD) [14]. The output of the method is an orthonormal basis of a low-dimensional subspace containing the solution allowing at least to construct its approximation. For now, let us assume that we have already found the basis, and see how it can help to work with the aggregation equations (2).
Let us start by rewriting (2) in a more general form. Namely, we start by introducing a tensor S ∈ R N ×N ×N : where δ i,j is the Kronecker symbol. Armed with this tensor, we rewrite (2) as Further, we assume the existence of an orthonormal basis, gathered as columns of a matrix V ∈ R N ×R , such that where n(t) is the solution to (4). We will hereafter abbreviate inequalities of this sort to n(t) ≈ V V T n(t).
Then we can introduce so that the equation (5) Multiplying this last system by V T and rearranging the sums a bit we arrive at a reduced form of the original system (note that doing so does not increase the second-norm absolute error, although it may well increase the relative one): or introducing some extra notatioñ we rewrite it as d dt Finally, instead of defining x via n, we can recast (11) as a system of ordinary differential equations for a new variablex, that approximates x: The important thing to notice here is that evaluation of the right-hand part of (12) only takes O(R 3 ) operations. Hence, we reach our initial goal of completely decoupling the dimensionality of the reduced system from N . If R N it may lead to a significant speedup of computations.
The reduced solutionx can then be used to reconstruct an approximation to the full solution by further approximating the original equation (5):

Constructing a basis
In this section we describe a method used to construct the basis V which we have been using in the previous section. To fullfill this aim we use the snapshot method from [14].
In this method, the basis V is constructed via the snapshots of the original solution at some fixed moments in time t k , for k = 1, . . . , m. The exact method of basis construction may vary in technical details in different publications about its applications, but the specific method from [14] ends up being equivalent to taking leading left singular vectors of the N × m matrix of snapshots. Specifically, in our case, V is taken to be a matrix of senior left singular vectors of a matrix composed of 'snapshots' n(t k ), with time moments t k uniformly spaced across the interval of interest. The number of singular vectors depends on the specified approximation requirements; in practice, we use the same criteria as when combining bases (see below).
Unfortunately, we essentially need to know the solution for construction a reduced basis which we are going to to use to find of the approximation of the solution. To resolve this circularity, we split the initial time-interval into a number of 'windows', and use the snapshot method to construct a basis for each of them in turn instead of finding just one basis for the entire time segment of our interest.
Specifically, let τ be some fixed time-window width, and assume we haveV k such that Each of these can be constructed via the method of snapshots by numerically solving of the full system (2) at each 'window' in turn by use of any standard numerical method for ODE systems.
To combine them into a final, common basis, we introduce an auxiliary operation ⊕ δ for any given δ > 0: given two matrices A ∈ R N ×r 1 and B ∈ R N ×r 2 , A ⊕ δ B is a N × r 3 matrix composed of the senior r 3 left singular vectors of an N × (r 1 + r 2 ) matrix C = (A | B) (that is, a matrix composed of columns of A and B -in principle, in any order), where r 3 is chosen so that σ r 3 ≥ δ > σ r 3 +1 , where σ k are singular values of C.
As a measure of the quality of our basis, we measure an error of approximation of the next window's basis by the "current" one, with some small positive tolerance ε 1. In other words, our algorithm for the basis construction can be formulated as following: Step 1 Set k ← 1, V 0 = 0 ∈ R N ×0 .
Step 2 CalculateV k via the method of snapshots as an approximate reduction basis for the time span [(k − 1)τ, kτ ].
Step 5 Set k ← k + 1 and repeat from step 2.
This algorithm is 'greedy' in some sense: it tries to approximate the entire solution by aggressively approximating each subsequent time 'window'. Hence, it probably may lead to overestimation of the eventual basis dimensionality. In our concrete implementation we execute Step 4 only if the approximation error at Step 3 is larger than an additional auxiliary parameter ε > ε > 0.

Numerical experiments
In this section we present the tests of the implementation of our method from previous sections. In our simulations we use a classical Browninan-type kernel Even though problems with such kernel and its closest generalizations are rather well-studied by nowadays [17,18,[21][22][23] the exact analytical solutions for time-dependent cases are still unknown especially for the cases with steady oscillations [18,23]. Moreover, researchers are still interested in the exploitation of such kernels for practical modelling [24] and theoretical analysis [25] as well.
In [26], a fast numerical method is given for a Cauchy problem with this kernel, evaluating the right-hand side of the equation (2) in just O(N log N ) operations and we want to out-perform this approach. In [18,23], steady collective oscillations in time of n k were detected for systems with this family of kernels with a > 0.5.
We have chosen this system specifically because the presence of the cycles in the solution gives us a strong a priori reason to expect that our method works. Namely, at least after the first iteration of the cycle, any basis which adequately approximates the solution should also approximate the further solution and can be used at least to verify the cyclical behaviour. However, we also note that due to the use of two thresholds ε > ε > 0 we cannot guarantee that the algorithm terminates, but, as we soon see in our experiments, it frequently does.
At first, we demonstrate the principle feasibility and inner workings of the algorithm. For this purpose, we present a sequence of experiments with the following set of model parameters: where δ is for our 'basis addition operator' ⊕ δ from Section 4, with m = 65 snapshots in each window for the snapshot method. As an ODE solver, we utilize a classical explicit midpoint time-integration method with a time-step of 2 −12 ≈ 2.4 × 10 −4 . In the full system, we evaluate right-hand side the via a fast method from [26]. Figure 1 demonstrates the inner working process of the algorithm during the basis construction for kernel parameter a equal to 0.6, 0.7 or 0.8. We can see that the projection error decreases rapidly at the onset, then oscillates a bit around ε , and eventually crosses the ε boundary. This, incidentally, highlights another role of the ε parameter -it effectively prevents the algorithm from over-approximating an initial segment of the solution; we show below the reason why this is important. Figure 2 demonstrates an error for a recomputed solution of a reduced system (12) for the time segment t ∈ [0, 256], as compared to a solution of a full system (2). As can be seen from the figure, the error remains small on an interval where the basis was originally constructed, and sharply increases at its end (where the algorithm effectively switches from interpolation to extrapolation)but, crucially, it still remains bounded around acceptable level 10 −3 .
For the same set of simulations, Table 1 provides the CPU time required to solve full and reduced systems, as well as the eventual size of the basis. As the table clearly demonstrates, the use of the reduction is not always beneficial, especially if fast algorithms for evaluation of the full operator are available; specifically, in the case a = 0.6, with the basis size of 216, the reduced system is actually more expensive to solve than the full problem.
Since this observation, together with Figure 2, strongly hint at a trade-off between performance and precision, one might be tempted to tweak the ε parameter to manage it. Unfortunately, as our second set of experiments demonstrates, this   is not always straight-forward in practice, which we find rather surprising. These experiments are performed with N = 65536 (we find that the effect is more visible at this dimensionality), with δ = ε and ε = 10 3 · ε. The results are available in Table 2; all the other parameters are the same as in the case above. As can be readily observed, the resulting error does not depend monotonically on ε, and seems to depend more on the actual time span which was used to construct the basis. Note that in the very last row, corresponding to ε = 10 −14 , the algorithm has simply used up the entire time span under evaluation for basis construction, and therefore the error reflects the 'interpolation' mode, as seen in Figure 2.
Finally, to test the scalability of our approach, we have tested the algorithm with a larger system with N = 131072 and a ∈ {0.6, 0.7}. Figure 3 demonstrates the relative error of the reconstructed solution for ε = 10 −10 , ε = 10 −12 and ε = 10 −14 , on a longer time segment of t ∈ [0, 512]. The same sharp 'interpolation-extrapolation' transition is visible here; and results from Table 3 confirm that, in both cases, the 'interpolation' region is in the neighbourhood of the transition visible on the graph, except for a = 0.6, ε = 10 −14 , where almost the entire time span was used for basis construction, and thus there is no visible transition at all.
Finally, Figure 4 demonstrates full and reduced solutions at the far end of the simulation time-interval. The solution with a = 0.6, ε = 10 −14 is indistinguishable  Table 3: Solution time and basis information for N = 131072; the test with a = 0.6 and ε = 10 −14 was run on a different hardware than the rest.
from the precise one, as it is effectively computed in 'interpolation' mode, but even the solution with a = 0.7, ε = 10 −14 is fairly quantitatively close to the full solution, diverging only for the smallest values concentration of n k -which, surprisingly, do not affect significantly the solution to the either side of their mass range. Less precise solutions with ε = 10 −10 and ε = 10 −12 do not deliver a good quantitative fit (as can already be seen in Figure 3) but, nevertheless, reproduce the qualitative shape of the solution well, despite being significantly cheaper to compute in terms of CPU-time.

Conclusion
We have suggested an application of the popular and well-established method of model reduction to the problem of a system of Smoluchowski ODEs, including a candidate method for construction of a reduced basis from an automatically selected prefix of the modelled time span. In our numerical experiments, we demonstrate the existence of such a low-dimensional basis, noticeable speed-ups of computations, and reasonable approximation to the full solution by the reduced model.
At the same time, we also demonstrate problematic sides of the chosen approach -the precision control of the reduced solution seems to be not straightforward at all due to the nonlinearity of both method and model. Hence, control of accuracy requires more theoretical analysis. In light of these shortcomings, we find our concept very promising for future development in more complicated applied cases and also consider it as fruitful directions of further research.

Acknowledgements
We are grateful to Nikolai Zamarashkin for comprehensive discussions during preparation of this work. The work was supported by the Russian Science Foundation, grant 19-11-00338.