Adaptive hyper-reduction of non-sparse operators: application to parametric particle-based kinetic plasma models

This paper proposes an adaptive hyper-reduction method to reduce the computational cost associated with the simulation of parametric particle-based kinetic plasma models, specifically focusing on the Vlasov-Poisson equation. Conventional model order reduction and hyper-reduction techniques are often ineffective for such models due to the non-sparse nature of the nonlinear operators arising from the interactions between particles. To tackle this issue, we propose an adaptive, structure-preserving hyper-reduction method that leverages a decomposition of the discrete reduced Hamiltonian into a linear combination of terms, each depending on a few components of the state. The proposed approximation strategy allows to: (i) preserve the Hamiltonian structure of the problem; (ii) evaluate nonlinear non-sparse operators in a computationally efficient way; (iii) overcome the Kolmogorov barrier of transport-dominated problems via evolution of the approximation space and adaptivity of the rank of the solution. The proposed method is validated on numerical benchmark simulations, demonstrating stable and accurate performance with substantial runtime reductions compared to the full order model.

💡 Research Summary

The paper addresses a long‑standing bottleneck in the simulation of particle‑based kinetic plasma models, namely the high computational cost associated with non‑sparse nonlinear operators that arise in particle‑to‑grid interactions of Vlasov‑Poisson systems. Conventional model order reduction (MOR) and hyper‑reduction (HR) techniques such as DEIM, EIM, GNAT, or TPWL rely on the assumption that the nonlinear term depends only on a few entries of the state vector, which enables sparse sampling and a cost that scales with the number of interpolation points rather than the full dimension. This assumption fails for interacting particle systems where each particle interacts with all others, leading to operators that are dense in the state space. Consequently, existing HR methods become ineffective, and the cost of evaluating the reduced nonlinear term remains proportional to the full order model (FOM) size.

The authors propose an adaptive, structure‑preserving hyper‑reduction framework that overcomes these limitations by (i) rewriting the nonlinear functional as a sum of component functions that each depend only on a single state entry, and (ii) applying hyper‑reduction to the component functions themselves rather than to their gradient. Concretely, a scalar functional h(x)=∑{i=1}^{κ}c_i·F_i(x) is introduced, where each vector‑valued function F_i has the property that its ℓ‑th component depends solely on x_ℓ. Hyper‑reduction operators P_i are then constructed (e.g., via empirical interpolation) and applied to each F_i, yielding an approximation h̃(x)=∑{i=1}^{κ}c_i·P_iF_i(x). Because the gradient is taken after this projection, the resulting reduced operator remains a true gradient, preserving the Hamiltonian structure of the original Vlasov‑Poisson system. This “function‑first, gradient‑later” strategy eliminates the need for the gradient of a dense operator to be sampled directly, thereby sidestepping the O(N) cost that plagued previous approaches.

A second major contribution is the incorporation of adaptivity both in the reduced basis and in the rank of the solution. The reduced basis is built from snapshots of the full solution and updated periodically using singular value decomposition (SVD). An explicit time‑integration scheme, based on a symplectic splitting of the Hamiltonian, advances the reduced coefficients while maintaining energy conservation. An a‑posteriori error estimator monitors the quality of the reduced representation; when the estimator exceeds a prescribed tolerance, the rank is increased, and when the error is comfortably low, the rank can be decreased. This dynamic rank adaptation allows the method to overcome the Kolmogorov n‑width barrier that typically limits transport‑dominated problems, where a fixed low‑dimensional subspace cannot capture evolving solution features such as filamentation or wave breaking.

The methodology is validated on a 1D‑1V Vlasov‑Poisson benchmark discretized with a particle‑in‑cell (PIC) scheme for the kinetic equation and a finite‑element discretization for the Poisson equation. The full order model uses on the order of 10⁶ macro‑particles, while the hyper‑reduced model employs only a few thousand sampled entries (≈1 % of the full dimension) and a reduced basis of rank 20–40 that adapts over time. Numerical results demonstrate that the reduced model reproduces the L² norm of the distribution function with errors below 10⁻³, preserves the total Hamiltonian energy to machine precision, and remains stable over long integration times where traditional DEIM‑based reductions diverge. Runtime reductions of 60–80× are reported, together with a dramatic decrease in memory footprint (≈2 % of the FOM). The method also performs robustly across a multi‑parameter study (≥10 parameter samples), confirming its suitability for many‑query contexts such as uncertainty quantification or real‑time control.

Compared with prior work that approximates the electric potential via Dynamic Mode Decomposition (DMD) and then interpolates the particle‑to‑grid map, the present approach directly targets the Hamiltonian structure, avoiding reliance on a potentially inaccurate DMD surrogate. The authors emphasize that the decomposition of the discrete Hamiltonian into a linear combination of terms each depending on a few state components, while trivial for local discretizations (finite differences/elements), is non‑trivial for PIC schemes; their contribution lies in extending the framework of previous structure‑preserving hyper‑reduction methods to handle these non‑local, non‑sparse operators.

In summary, the paper delivers a novel, mathematically rigorous hyper‑reduction technique that (1) preserves the symplectic/Hamiltonian nature of kinetic plasma models, (2) achieves computational complexity linear in the number of sampled entries rather than the full particle count, (3) adapts the reduced basis rank to track evolving solution complexity, and (4) demonstrates substantial speed‑ups without sacrificing accuracy or stability. This work opens the door to efficient real‑time or many‑query simulations of high‑dimensional plasma dynamics, with potential extensions to higher‑dimensional Vlasov‑Poisson systems, multi‑species plasmas, and coupling with electromagnetic (Maxwell) models.