M-GAPS: Scalable GAP Fitting

• M-GAPS is a data-parallel software package that scales Gaussian Approximation Potentials fitting by distributing descriptor evaluations and QR solves across many nodes.
• It separates the workload into a fully local descriptor and kernel calculation phase and a distributed, ScaLAPACK-based least-squares QR factorization phase.
• The approach dramatically reduces training time and memory constraints, enabling efficient hyperparameter optimization and active learning in complex materials modelling.

Searching arXiv for the M-GAPS paper and closely related Gaussian Approximation Potentials/SOAP references. M-GAPS, short for “Massively Parallel Fitting of Gaussian Approximation Potentials,” is a data-parallel software package for fitting Gaussian Approximation Potentials (GAPs) on multiple nodes using the ScaLAPACK library with MPI and OpenMP. It addresses a central bottleneck in large-scale GAP training: until this implementation, the maximum training set size for GAP models had been limited by the available memory on a single compute node. In M-GAPS, descriptor evaluation is carried out in parallel with no communication requirement, while the subsequent linear solve required to determine the model coefficients is parallelised with ScaLAPACK. The approach scales to thousands of cores, lifts the single-node memory limitation, and expands the applicability of the GAP approach to more complex systems as well as to higher-level workflows such as committee models or hyperparameter optimisation (Klawohn et al., 2022).

1. Problem setting and model class

GAP approximates the total energy of a configuration $A$ by a sum of local contributions,

$E_A \;=\;\sum_{i=1}^{N_A} \varepsilon(x_i)\,,$

where each local energy is represented by a sparse Gaussian-process expansion over $M$ representative descriptors $\{x_j\}$,

$$\varepsilon(x_i)\;=\;\sum_{j=1}^M c_j\,k(x_i,x_j)\,.$$

Here $k(x,x')$ measures the similarity of two local environments. In the most common SOAP variant, one first expands each atomic neighbor density in spherical harmonics and radial basis functions, collects power spectra, and then defines

$$k_\mathrm{SOAP}(x,x') \;=\;\bigl(\hat{p}(x)\cdot \hat{p}(x')\bigr)^\zeta,$$

with $\hat p$ the normalized power spectrum and $\zeta$ a tunable exponent (Klawohn et al., 2022).

The computational challenge arises from fitting on large collections of ab initio observations comprising energies, forces, and stresses. In that regime, the descriptor stage and the dense linear-algebra stage both become substantial. M-GAPS is specifically organized to separate a fully local descriptor/kernel phase from a distributed QR-based least-squares phase. This division is fundamental to the scheme’s scalability because the first phase is embarrassingly parallel, whereas the second is delegated to ScaLAPACK.

2. Regression formulation and linear system

Training uses a Tikhonov-regularized loss for the coefficient vector $c$,

$E_A \;=\;\sum_{i=1}^{N_A} \varepsilon(x_i)\,,$0

with $E_A \;=\;\sum_{i=1}^{N_A} \varepsilon(x_i)\,,$1 and a diagonal noise-precision matrix $E_A \;=\;\sum_{i=1}^{N_A} \varepsilon(x_i)\,,$2. This leads to the sparse-GP normal equations

$E_A \;=\;\sum_{i=1}^{N_A} \varepsilon(x_i)\,,$3

M-GAPS also uses an equivalent least-squares formulation. Writing

$E_A \;=\;\sum_{i=1}^{N_A} \varepsilon(x_i)\,,$4

where $E_A \;=\;\sum_{i=1}^{N_A} \varepsilon(x_i)\,,$5, the weights are obtained from

$E_A \;=\;\sum_{i=1}^{N_A} \varepsilon(x_i)\,,$6

after the numerically stable QR factorization $E_A \;=\;\sum_{i=1}^{N_A} \varepsilon(x_i)\,,$7 (Klawohn et al., 2022).

This reformulation is important because it converts sparse-GP fitting into a dense distributed least-squares problem with a well-defined matrix shape, $E_A \;=\;\sum_{i=1}^{N_A} \varepsilon(x_i)\,,$8, and makes the solver stage directly compatible with ScaLAPACK’s QR routines.

3. Data-parallel descriptor evaluation

The MPI stage splits the $E_A \;=\;\sum_{i=1}^{N_A} \varepsilon(x_i)\,,$9 target properties—energies, forces, and virials—by distributing whole atomic configurations across $M$0 MPI ranks. A small greedy heuristic first sorts configurations by their number of targets and then assigns them to ranks to balance the total $M$1 per rank. Each MPI rank reads only its subset of structures, builds the local descriptor vectors $M$2, and, if forces are included, their gradients $M$3 (Klawohn et al., 2022).

Because each rank holds the same set of $M$4 sparse descriptors $M$5, evaluation of the local kernel blocks

$M$6

becomes completely embarrassingly parallel: no inter-rank communication is needed until the global matrices must be assembled. Memory is saved by computing and immediately contracting each descriptor block into the global matrix $M$7, while discarding the large per-structure descriptor arrays as soon as they have been used (Klawohn et al., 2022).

A direct consequence is that the descriptor/covariance stage is all-local. The paper reports that this part is virtually perfectly scalable, a property that contrasts with the later QR stage, where collective communication becomes unavoidable.

4. Distributed QR solve and process layout

The full matrix $M$8 is formed by concatenating the scaled $M$9 rows with the $\{x_j\}$0 lower-triangular block $\{x_j\}$1. In ScaLAPACK terms, M-GAPS uses a $\{x_j\}$2 process grid and block-cyclically distributes rows of $\{x_j\}$3 so that each MPI rank holds approximately $\{x_j\}$4 rows. A small number of zero-row padding is inserted so that each local submatrix size is a multiple of the chosen block size (Klawohn et al., 2022).

The solver sequence is explicit:

1. pdgeqrf: compute the QR factorization $\{x_j\}$5;
2. pdormqr: apply $\{x_j\}$6 to $\{x_j\}$7, yielding $\{x_j\}$8;
3. pdtrtrs: solve the triangular system $\{x_j\}$9 by back-substitution.

Block-size tuning is nontrivial. pdgeqrf and pdormqr prefer a moderate column block size, approximately $$\varepsilon(x_i)\;=\;\sum_{j=1}^M c_j\,k(x_i,x_j)\,.$$0, to trade off local BLAS performance against inter-node communication. pdtrtrs enforces equal row-column block sizes; by switching the block-size setting between the second and third calls, the implementation minimizes total extra memory overhead (Klawohn et al., 2022).

The formal cost of the QR factorization scales as

$$\varepsilon(x_i)\;=\;\sum_{j=1}^M c_j\,k(x_i,x_j)\,.$$1

once the data are distributed, plus a small

$$\varepsilon(x_i)\;=\;\sum_{j=1}^M c_j\,k(x_i,x_j)\,.$$2

Cholesky of $$\varepsilon(x_i)\;=\;\sum_{j=1}^M c_j\,k(x_i,x_j)\,.$$3 on rank $$\varepsilon(x_i)\;=\;\sum_{j=1}^M c_j\,k(x_i,x_j)\,.$$4. Communication is limited to the collective operations inside pdgeqrf and pdormqr; the earlier data-parallel kernel part has zero communication overhead (Klawohn et al., 2022).

5. Strong and weak scaling results

The implementation was tested on the Raven cluster, using 72-core Intel IceLake nodes with 100 Gb/s IB, on two real-world fitting problems: a MoNbTaVW high-entropy alloy dataset and a silicon carbide dataset (Klawohn et al., 2022).

System	Core problem size	Reported best multi-node result
HEA (MoNbTaVW high-entropy alloy)	$$\varepsilon(x_i)\;=\;\sum_{j=1}^M c_j\,k(x_i,x_j)\,.$$5, $$\varepsilon(x_i)\;=\;\sum_{j=1}^M c_j\,k(x_i,x_j)\,.$$6	$$\varepsilon(x_i)\;=\;\sum_{j=1}^M c_j\,k(x_i,x_j)\,.$$7 speedup on 64 nodes, wall-time $$\varepsilon(x_i)\;=\;\sum_{j=1}^M c_j\,k(x_i,x_j)\,.$$8 s
SiC (silicon carbide)	$$\varepsilon(x_i)\;=\;\sum_{j=1}^M c_j\,k(x_i,x_j)\,.$$9, $k(x,x')$0	$k(x,x')$1 speedup on 64 nodes, wall-time $k(x,x')$2 s

For the HEA system, the number of sparse points was $k(x,x')$3-body $k(x,x')$4SOAP $k(x,x')$5. The matrix $k(x,x')$6 had $k(x,x')$7 elements, corresponding to about $k(x,x')$8 GB in double precision. The single-node reference run, using $k(x,x')$9 MPI processes and $$k_\mathrm{SOAP}(x,x') \;=\;\bigl(\hat{p}(x)\cdot \hat{p}(x')\bigr)^\zeta,$$0 OpenMP threads on $$k_\mathrm{SOAP}(x,x') \;=\;\bigl(\hat{p}(x)\cdot \hat{p}(x')\bigr)^\zeta,$$1 cores, took $$k_\mathrm{SOAP}(x,x') \;=\;\bigl(\hat{p}(x)\cdot \hat{p}(x')\bigr)^\zeta,$$2 s. The maximum speedup on $$k_\mathrm{SOAP}(x,x') \;=\;\bigl(\hat{p}(x)\cdot \hat{p}(x')\bigr)^\zeta,$$3 nodes was approximately $$k_\mathrm{SOAP}(x,x') \;=\;\bigl(\hat{p}(x)\cdot \hat{p}(x')\bigr)^\zeta,$$4, reducing the wall-time to about $$k_\mathrm{SOAP}(x,x') \;=\;\bigl(\hat{p}(x)\cdot \hat{p}(x')\bigr)^\zeta,$$5 s. Per-MPI-rank memory was reported as approximately $$k_\mathrm{SOAP}(x,x') \;=\;\bigl(\hat{p}(x)\cdot \hat{p}(x')\bigr)^\zeta,$$6–$$k_\mathrm{SOAP}(x,x') \;=\;\bigl(\hat{p}(x)\cdot \hat{p}(x')\bigr)^\zeta,$$7 GB, growing linearly with the number of ranks, while total memory at $$k_\mathrm{SOAP}(x,x') \;=\;\bigl(\hat{p}(x)\cdot \hat{p}(x')\bigr)^\zeta,$$8 nodes was about $$k_\mathrm{SOAP}(x,x') \;=\;\bigl(\hat{p}(x)\cdot \hat{p}(x')\bigr)^\zeta,$$9 TB (Klawohn et al., 2022).

For the SiC system, $\hat p$0 with the same $\hat p$1. The matrix $\hat p$2 had approximately $\hat p$3 elements, about $\hat p$4 GB. The single-node reference run, using $\hat p$5 MPI processes and $\hat p$6 threads, took $\hat p$7 s. The best speedup was approximately $\hat p$8 on $\hat p$9 nodes, with $\zeta$0 MPI processes and $\zeta$1 threads per node, reducing the run time to about $\zeta$2 s. Peak total memory reached about $\zeta$3 TB for pure-MPI runs with one thread on $\zeta$4 nodes (Klawohn et al., 2022).

In both cases, the covariance-matrix build was virtually perfectly scalable because it was all-local, whereas the QR solve exhibited the usual saturation once inter-node communication began to dominate. This pattern is central to the performance profile of M-GAPS.

6. Workflow role, operating recommendations, and limitations

By distributing both memory and compute across many nodes, M-GAPS completely removes the single-node memory ceiling, splitting a $\zeta$5 TB problem into $\zeta$6–$\zeta$7 GB per rank. The fitting step, which once took days on a fat-node, can instead be done in minutes. The paper identifies three immediate workflow consequences: automated hyperparameter searches involving kernel widths, regularization strength, and SOAP cutoffs; active-learning loops or on-the-fly committee-model retraining; and integration into larger materials-design pipelines (Klawohn et al., 2022).

The implementation also gives specific operating guidance. Recommended practice is to keep at least some OpenMP threads per MPI rank, for example $\zeta$8–$\zeta$9 threads, in order to reduce per-rank memory and control the communication footprint. It is also recommended to distribute only up to approximately $c$0 the number of structures as MPI ranks, so that each rank still holds at least one structure’s worth of data (Klawohn et al., 2022).

Two practical limits are explicit. First, the QR stage remains communication-sensitive and therefore eventually saturates as node counts increase. Second, sparse-point selection, for example CUR clustering, is still carried out in a first, single-node sparsify_only pass; future work could hide this behind MPI as well (Klawohn et al., 2022).

7. Significance within GAP fitting

M-GAPS establishes a concrete HPC realization of sparse-GP GAP fitting in which the descriptor/covariance stage has zero communication overhead and the least-squares stage is handled by ScaLAPACK QR factorization. The resulting architecture is not merely an implementation detail: it determines the method’s scalability envelope, memory behavior, and suitability for repeated fitting inside larger workflows (Klawohn et al., 2022).

A plausible implication is that M-GAPS is most consequential in settings where the fitting stage is no longer an isolated terminal step but a repeatedly invoked component in active-learning, committee-model, and hyperparameter-optimization loops. Within that regime, the paper’s main contribution is the conversion of multi-node GAP fitting from a single-node memory-constrained procedure into a routine distributed workflow element with well-characterized performance tradeoffs between an all-local kernel stage and a communication-bound QR stage.