Over the past decades, the study of systems of particles has become an important part of many research areas, from theoretical physics to applied biology and computational mathematics. Some current research trends include the following.
- Since the 80's, there is a joint effort from biologists, physicists, computer graphics scientists and more recently mathematicians to propose accurate models for large self-organized animal societies. While the motivations of the different communities are largely independent, a commmon goal is to be able to explain and reproduce the emergence of complex patterns from simple interaction rules. Typical examples of complex systems include flocks of birds, fish schools, herds of mammals or ant colonies. It seems reasonable to consider that none of these animals has an accute consciousness of the whole system organization but rather follows simple behavioral patterns: stay close to the group, do not collide with another individual, copy the attitude of the close neighbours etc. These models are often called swarming models or collective dynamics models
vicsek_collective_2012, degond_mathematical_2018, albi_vehicular_2019, naldi_particle_2010. - The microscopic world if full of complex systems made of many simple entities. Colonies of bacteria
wensink_meso-scale_2012are a natural example reminiscent of the macroscopic animal societies described above. But complex patterns can also emerge from systems of non-living particles. Our body is actually a very complex self-organized assembly of cells which dictate every aspect of our life: our brain works thanks to the communications between billions of neuronsfournier_toy_2016, andreis_mckeanvlasov_2018, the reproduction is based on the competition between millions of spermatozoacreppy_turbulence_2015and sometimes, death is sadly caused by the relentless division of cancer cells. Unlike the previous models, the communication between the particles cannot be based on their visual perception, but rather on chemical agents, on physical (geometrical) constraints etc. - On a completely different side, there is an ever growing number of methods in computational mathematics which are based on the simulation of systems of particles. The pioneering Particle Swarm Optimization method
kennedy_particle_1995has shown how biologically-inspired artificial systems of particles can be used to solve tough optimization problems. Following these ideas, other Swarm Intelligence optimization algorithms have been proposed up to very recentlypinnau_consensus-based_2017, totzeck_numerical_2018, grassi_particle_2020, totzeck_trends_2021. Since the beginning of the 2000's, particle systems are also at the core of filteringcrisan_survey_2002, doucet_sequential_2001, kantas_overview_2009, del_moral_feynman-kac_2004, del_moral_mean_2013and sampling methodscappe_population_2004, clarte_collective_2021. Recently, a particle-based interpretationmei_mean_2018, chizat_global_2018, rotskoff_trainability_2019, de_bortoli_quantitative_2020, sirignano_mean_2020of the training task of neural networks has lead to new theoretical convergence results.
Beyond the self-explanatory applications for particle-based algorithms in computational mathematics, the simulation of systems of particles is also a crucial modelling tool in Physics and Biology.
- On the one hand, field experiments are useful to collect data and trigger new modelling ideas. On the other hand, numerical simulations become necessary to test these ideas and to calibrate the models. Numerical experiments can be conducted to identify which mechanisms are able to produce a specific phenomena in a controlled environment
chate_collective_2008. - On a more theoretical side, the direct mathematical study of particle systems can rapidly become incredibly difficult. Inspired by the kinetic theory of gases, many mesoscopic and macroscopic models have been proposed to model the average statistical behavior of particle systems rather than the individual motion of each particle
toner_flocks_1998, degond_continuum_2008, naldi_particle_2010. Within this framework, a particle system is rather seen as a fluid and it is decribed by Partial Differential Equations (PDE) reminiscent from the theory of fluid dynamics. PDE models are more easily theoretically and numerically tractable. However, cheking the validity of these models is not always easy and is sometimes only postulated based on phenomenological considerations. In order to design good models, it is often necessary to go back-and-forth between the PDE models and the numerical simulation of the underlying particle systems.
The development of the SiSyPHE library was initially motivated by the study of body-oriented particles <sisyphe.particles.BOParticles> giacomin_alignment_2019. The (formal) derivation of a macroscopic PDE model from the particle system has lead to a novel conjecture which postulates the existence of a class of so-called bulk topological states degond_bulk_2021. The quantitative comparison between this theoretical prediction and the numerical simulation of the particle system in a suitable regime (with more than 106 particles) has confirmed the existence of these new states of matter. The study of their physical properties which are observed in the numerical experiments but not readily explained by the PDE model is an ongoing work.
Currently, the models implemented in the SiSyPHE library belong to the family of mean-field models. It means that the motion of each particle is influenced by the average behavior of the whole system. The Vicsek model vicsek_novel_1995, degond_continuum_2008 and the Cucker-Smale model cucker_mathematics_2007, ha_emergence_2009, naldi_particle_2010 are two popular examples of mean-field models where each particle tries to move in the average direction of motion of its neighbours (it produces a so-called flocking behavior).
From a mathematical point of view, a mean-field particle system with N particles is defined as a Markov process in (ℝd)N with a generator ℒN whose action on a test function φN is of the form:
where xN = (x1, …, xN) ∈ (ℝd)N and
is the so-called empirical measure. The operator LμxN depends on the empirical measure and acts on one-variable test functions on ℝd. Given an operator L and a test function φN, the notation L⋄iφN denotes the function
L⋄iφN : (x1, …, xN) ∈ (ℝd)N ↦ L[x ↦ φN(x1, …, xi − 1, x, xi + 1, …, xN)](xi) ∈ ℝ.
The dynamics of the mean-field system depends on the operator LμxN which is typically either a diffusion operator degond_continuum_2008, ha_emergence_2009 or a jump operator dimarco_self-alignment_2016, andreis_mckeanvlasov_2018. In the first case, the mean-field particle systems can be alternatively defined by a system of N coupled Stochastic Differential Equations of the form:
dXti = b(Xti, μ𝒳tN)dt + σ(Xti, μ𝒳tN)dBti, i ∈ {1, …, N},
where (Bti)t are N independent Brownian motions and the coefficients b and σ are called the drift and diffusion coefficients. In most cases, these coefficients are linear or quasi-linear which means that they are of the form
for two given functions b̃ : ℝd × ℝn → ℝd and K : ℝd × ℝd → ℝn.
When the particles are initially statistically independent, it can be shown that when N → + ∞, each particle Xti converges towards an independent copy of the solution of the so-called McKean-Vlasov diffusion process mckean_propagation_1969, sznitman_topics_1991, meleard_asymptotic_1996 defined by the Stochastic Differential Equation
where (Bt)t is a Brownian motion and ft is the law of the process
with x = (x1, …, xd) ∈ ℝd and a = σσT. This phenomenon is called propagation of chaos, following the terminology introduced by Kac in the 50's kac_foundations_1956, mckean_propagation_1969, sznitman_topics_1991, meleard_asymptotic_1996, hauray_kacs_2014.
Note
A popular example of mean-field particle system is the (stochastic) Cucker-Smale model cucker_mathematics_2007, ha_emergence_2009, naldi_particle_2010. Each particle is defined by its position Xti ∈ ℝd and its velocity Vti ∈ ℝd which evolve according to the system of 2N Stochastic Differential Equations:
where K : [0, + ∞) → [0, + ∞) is an observation kernel which models the visual perception of the particles. In this example, the function K is a smooth function vanishing at infinity and the communication between the particles is based on the distance between them. The motion of the particles follows the Newton's laws of motion with an additional stochastic term. The term Vtj − Vti is a relaxation force (for a quadratic potential) which tends to align the velocities of particle i and particle j.
On a computer, the above examples which are time-continuous needs to be discretized, using for instance one of the classical numerical schemes for Stochastic Differential Equations. Then the difficulty lies in the evaluation of the empirical measure at each time-step of the numerical scheme.
In the example of the Cucker-Smale model, the force exerted on a particle is the sum of N small relaxation forces of order 1/N. The total number of operations required is thus of order 𝒪(N) for each particle. Since there are N particles, the total time complexity of the algorithm is thus 𝒪(N2T) where T is the total number of iterations of the numerical scheme.
This quadratic cost is the main bottleneck in the simulation of mean-field particle systems. As explained in the documentation of the KeOps library, the evaluation of the N forces at time t
is called a kernel operation and can be understood as a discrete convolution operation (or matrix-vector product) between the matrix of distances (Kij)i, j ∈ {1 …, N} where Kij = K(|Xtj − Xti|) and the vector of velocities (Vtj − Vti)j ∈ {1, …, N}. When N is large (say N > 104), such operation is too costly even for array-based programming languages such as Matlab: the N × N kernel matrix (Kij)i, j would simply not fit into the memory. With lower level languages (Fortran, C), this operation can be implemented more efficiently with two nested loops but with a significantly higher global coding effort and with less versatility.
Over the past decades, several workarounds have been proposed. Popular methods include
- the low-rank decomposition of the kernel matrix,
- the fast-multipole methods
greengard_fast_1987which to treat differently short- and long-range interactions, - the Verlet list method which is based on a grid decompostion of the spatial domain to reduce the problem to only short-range interactions between subsets of the particle system,
- the Random Batch Method
jin_random_2019which is based on a stochastic approximation where only interactions between randomly sampled subsets (batches) of the particle system are computed.
All these methods require an significant amount of work, either in terms of code or to justify the approximation procedures.
The present implementation is based on recent libraries originally developed for machine learning purposes to significantly accelerate such tensor (array) computations, namely the PyTorch package and the KeOps library charlier_kernel_2021. Using the KeOps framework, the kernel matrix is a symbolic matrix defined by a mathematical formula and no approximation is required in the computation of the interactions (up to the time discretization). The SiSyPHE library speeds up both traditional Python and low-level CPU implementations by one to three orders of magnitude for systems with up to several millions of particles. Although the library is mainly intended to be used on a GPU, the implementation is fully functional on the CPU with a significant gain in efficiency.
Moreover, the versatile object-oriented Python interface is well suited to the comparison and study of new and classical many-particle models. This aspect is fundamental in applications in order to conduct ambitious numerical experiments in a systematic framework, even for particles with a complex structure and with a significantly reduced computational cost