# Nonlinear reduced basis using mixture Wasserstein Barycenters: application to an eigenvalue problem inspired from quantum chemistry

Maxime Dalery, Geneviève Dusson, Virginie Ehrlacher, Alexei Lozinski

► To cite this version:

Maxime Dalery, Geneviève Dusson, Virginie Ehrlacher, Alexei Lozinski. Nonlinear reduced basis using mixture Wasserstein Barycenters: application to an eigenvalue problem inspired from quantum chemistry. 2025. hal-05401051

HAL Id: hal-05401051

https://hal.science/hal-05401051v1

Preprint submitted on 5 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# NONLINEAR REDUCED BASIS USING MIXTURE WASSERSTEIN BARYCENTERS: APPLICATION TO AN EIGENVALUE PROBLEM INSPIRED FROM QUANTUM CHEMISTRY

MAXIME DALERY*, GENEVIÈVE DUSSON*, VIRGINIE EHRLACHER†, AND ALEXEI LOZINSKI*

**Abstract.** The aim of this article is to propose a new reduced order modelling approach for parametric eigenvalue problems arising in electronic structure calculations. Namely, we develop nonlinear reduced basis techniques for the approximation of parametric eigenvalue problems inspired from quantum chemistry applications. More precisely, we consider here a one-dimensional model which is a toy model for the computation of the electronic ground state wavefunction of a system of electrons within a molecule, solution to the many-body electronic Schrödinger equation, where the varying parameters are the positions of the nuclei in the molecule. We estimate the decay rate of the Kolmogorov *n*-width of the set of solutions for this parametric problem in several settings, including the standard $L^2$-norm as well as with distances based on optimal transport. The fact that the latter decays much faster than in the traditional $L^2$-norm setting motivates us to propose a practical nonlinear reduced basis method, which is based on an offline greedy algorithm, and an efficient stochastic energy minimization in the online phase. We finally provide numerical results illustrating the capabilities of the method and good approximation properties, both in the offline and the online phase.

**Key words.** reduced basis, eigenvalue problem, optimal transport, Wasserstein barycenters

**AMS subject classifications.** 65D05, 65K10, 41A63, 60B05, 47N50, 47A75

**1. Introduction.** In many academic and industrial applications, model order reduction techniques are used to accelerate the computation of the solutions to parametric partial differential equations. Techniques such as the reduced basis method give outstanding results for a large class of problems, see [14, 20]. A key factor determining the success of the method is the approximability of the solution by a linear combination in a fixed vector space, possibly spanned by solutions for specific values of the parameters called snapshots. This ability is characterized by a fast-decreasing of the so-called Kolmogorov *n*-width, as described below. Such approach works very well in numerous cases, such as for linear elasticity equations [16], thermal equations [24], see also [21] and references therein. However, as was pointed out in [10], this approach proves ineffective in several cases, especially when the solutions exhibit some transport of mass over parameter or time variation. This is, for instance, the case of the pure transport equation [18]. As the solution is translated over time, it is ineffective to approximate the solution as a linear combination of previous time steps, which would be the standard approximation using a linear reduced basis method. Another illustrative example is the electronic structure problem – an eigenvalue problem where the solution tends to localize around the nuclei. This will be the problem of interest in this article. Also, again in [10], it is shown that for a simple Burgers equation, the Kolmogorov *n*-width decreases faster if one uses not a linear combination of previous snapshots, but Wasserstein barycenters between solutions, i.e. if there is some nonlinear transformation involved.

Recently, several works have built on this idea to propose nonlinear interpolations between solutions based on optimal transport. This is for example the case in [15], where the authors propose a method based on an affine transformation of the snapshots to construct new approximations, or in [17], where a preprocessing step using optimal transport is added to the offline phase. In [8], a new method based on sparse Wasserstein barycenters is proposed. Other works include some machine-learning techniques to construct the nonlinear map, see e.g. [23], or to reconstruct higher frequency modes of the linear reduced basis operator from the low frequency modes, using trees or random forests [5].

The main limitation to the works based on optimal transport is the computational cost of Wasserstein barycenters, which do not scale well with the space dimension and the number of snapshots involved in the computation of the barycenter (solved through the so-called multi-marginal problem). A recent article [6] has proposed a modified Wasserstein distance between mixtures of Gaussians, which was extended for more general mixtures in [9]. This is particularly interesting for the electronic structure problem, where the solutions are often represented by a small number of functions of the same type, typically Gaussians or Slater functions. Indeed, this modified Wasserstein distance enables the computation of barycenters without suffering from the curse of dimensionality, since the problem dimension depends on the number of components

<footnote>* Université Marie et Louis Pasteur, CNRS, LmB (UMR 6623), F-25000 Besançon, France (genevieve.dusson@math.cnrs.fr)</footnote>

<footnote>† CERMICS, École des Ponts & Inria Paris, Paris, France</footnote>
in the mixtures, and not on a potentially large spatial grid parameter.

In this article, we fully use this mixture distance to propose a nonlinear reduced basis method based on optimal transport. As is standard in reduced basis methods, our algorithm works in two stages. In the offline stage, a collection of snapshots is gathered using a greedy algorithm, selecting, at each step, the worse-approximated snapshot in the set of barycenters between previously selected snapshots for this particular mixture distance. In the online stage, i.e. when one wants to compute the solution for a new parameter, we minimize the energy on the set of barycenters between selected snapshots. It is a nonlinear problem but in a low-dimensional parameter space, so that the online cost stays reasonable.

The outline of this article is as follows. In Section 2, we present the settings, namely we present the one-dimensional eigenvalue toy problem we mainly focus on in this work, and we detail a few preliminaries on the Kolmogorov n-width, as well as on optimal transport and Wasserstein barycenters. In Section 3 we state our main theoretical results about the estimations for the Kolmogorov n-width of the set of solutions of the one-dimensional toy problem, in different settings: with a linear approximation, with a Wasserstein transport metric, as well as a Wasserstein-type metric for mixtures. Our results are similar in spirit to those proved of [10], with the specificity that we focus on solution sets stemming from parametrized electronic structure calculation problems and that we consider the mixture Wasserstein metric [6, 9], in addition to the exact Wasserstein and L² metrics. The proofs of the theoretical results are postponed to Section 6. In Section 4, we present the nonlinear reduced basis method proposed in this work, which consists of an offline and an online stage. Let us emphasize that, despite the fact that we only illustrate its numerical behaviour on the one-dimensional toy problem mentioned above, the proposed numerical strategy can in principle be adapted to solve realistic electronic structure problems in any dimension. This is the focus of ongoing work. The main novelty of the approach, in particular compared to the one proposed in [10] is that, in the online phase, we use an energetic variational principle to determine the solution of the nonlinear reduced order model. Finally, we present numerical results in Section 5.

**2. Preliminaries.** The aim of this section is to introduce some preliminaries. We first present in Section 2.1 the considered parametric eigenvalue problem, which is motivated by quantum chemistry applications. We then provide some definitions about Kolmogorov widths in Section 2.2, and we recall some fundamentals about the Wasserstein and Mixture Wasserstein metrics in Section 2.3.

**2.1. An eigenvalue problem inspired from quantum chemistry.** In this article, we focus on the following one-dimensional eigenvalue partial differential equation parameterized by $\mathbf{r} := (r_1, \dots, r_M) \in \mathbb{R}^M$ and $\mathbf{z} := (z_1, \dots, z_M) \in (\mathbb{R}_+)^M$ for $M \in \mathbb{N}^*$. More precisely, we are looking for the lowest eigenvalue $E_{\mathbf{r},\mathbf{z}} \in \mathbb{R}$ and a corresponding eigenstate $u_{\mathbf{r},\mathbf{z}} \in H^1(\mathbb{R})$ satisfying

$$ (2.1) \qquad -\frac{1}{2}u''_{\mathbf{r},\mathbf{z}} + \left(-\sum_{m=1}^{M} z_m \delta_{r_m}\right) u_{\mathbf{r},\mathbf{z}} = E_{\mathbf{r},\mathbf{z}} u_{\mathbf{r},\mathbf{z}}. $$

This problem can be seen as a toy ground state electronic structure problem, with an Hamiltonian of the form $-\frac{1}{2}\Delta + V$, with a potential $V$ taken as a sum of Dirac masses $V := -\sum_{m=1}^{M} z_m \delta_{r_m}(x)$ localized at atomic positions $\mathbf{r}$ with charges $\mathbf{z}$. In this simple framework, it can be easily checked that any eigenvector solution to this problem actually belongs to $L^1(\mathbb{R})$. More precisely, there exists a unique strictly positive eigenvector solution to this problem such that $\|u_{\mathbf{r},\mathbf{z}}\|_{L^1(\mathbb{R})} = 1$, and it is explicitly given by [19, Section 3.1]

$$ (2.2) \qquad u_{\mathbf{r},\mathbf{z}} = \sum_{m=1}^{M} \pi_{m}^{\mathbf{r},\mathbf{z}} S_{\zeta_{\mathbf{r},\mathbf{z}},r_m}, $$

with $\pi_{m}^{\mathbf{r},\mathbf{z}} = (\pi_{m}^{\mathbf{r},\mathbf{z}})_{m=1}^{M} \in (\mathbb{R}_+)^M$ of total sum equal to 1, $\zeta_{\mathbf{r},\mathbf{z}} > 0$, and where for all $\zeta > 0$ and all $r \in \mathbb{R}$, the Slater function $S_{\zeta,r}$ is defined by

$$ (2.3) \qquad S_{\zeta,r} : x \mapsto \frac{\zeta}{2} e^{-\zeta|x-r|}. $$

Note that the normalization with respect to the $L^1$-norm in problem (2.1) is not standard, but $u_{\mathbf{r},\mathbf{z}}$ can thus be interpreted as the density associated with a probability measure on $\mathbb{R}$, which will be an essential
feature in the following. In the rest of the article, we will make an abuse of notation and identify an absolutely continuous probability measure with its associated probability density. We refer the reader to [19] for an extensive review on the link between this toy one-dimensional model and actual electronic structure calculation problems in molecules.

The eigenvalue problem (2.1) can be equivalently formulated as an energy minimization problem

$$ (2.4) \qquad \min_{\substack{u \in H^1(\mathbb{R}) \\ \|u\|_{L^1(\mathbb{R})} = 1}} \frac{E_{\mathbf{r},\mathbf{z}}(u)}{\|u\|_{L^2(\mathbb{R})}^2}, $$

with

$$ E_{\mathbf{r},\mathbf{z}}(u) := \frac{1}{2} \int_{\mathbb{R}} |u'|^2 - \sum_{m=1}^{M} z_m u(r_m)^2. $$

Let us mention here some particular explicit formulas available for small values of $M$.

**Case 1:** For $M = 1$, $\mathbf{r} = (r)$ and $\mathbf{z} = (z)$ for some $r \in \mathbb{R}$ and $z > 0$, it holds that (see [19, Theorem 3.4])

$$ (2.5) \qquad \zeta_{\mathbf{r},\mathbf{z}} = z. $$

**Case 2:** For $M = 2$, $\mathbf{r} = (-r, r)$ and $\mathbf{z} = (z, z)$ with $r, z > 0$, it holds that (see [19, Corollary 3.3])

$$ (2.6) \qquad \zeta_{\mathbf{r},\mathbf{z}} = z + \frac{W(2zre^{-2zr})}{2r} \quad \text{and} \quad \pi^{\mathbf{r},\mathbf{z}} = \left(\frac{1}{2}, \frac{1}{2}\right), $$

where $W$ is the Lambert function defined as the inverse of the function $x \mapsto xe^x$. In this special case, we also have the following equality on $\zeta_{\mathbf{r},\mathbf{z}}$:

$$ (2.7) \qquad \zeta_{\mathbf{r},\mathbf{z}} = z(1 + e^{-2\zeta_{\mathbf{r},\mathbf{z}}r}). $$

**2.2. Kolmogorov widths in metric spaces.** Let us now introduce some definitions, which will play the role of generalized Kolmogorov $n$-widths in our setting. We first start by recalling the definition of the Kolmogorov $n$-width in a Hilbert space $\mathbb{H}$ endowed with a scalar product $\langle \cdot, \cdot \rangle$ and associated norm $\| \cdot \|$. We denote by $P_V : \mathbb{H} \to V$ the projection operator onto a closed vector subspace $V \subset \mathbb{H}$.

We also denote by $\mathbf{Z} \subset \mathbb{R}^p$ a subset of parameter values, and for all $z \in \mathbf{Z}$, we assume that $u(z)$ is an element of $\mathbb{H}$. Finally, we denote by $\mathcal{E}$ the following set

$$ \mathcal{E} := \{u(z), z \in \mathbf{Z}\} \subset \mathbb{H}. $$

**Definition 2.1.** The $L^\infty$ Kolmogorov $n$-width of $\mathcal{E}$ is defined by

$$ (2.8) \qquad d_{\infty,n}(\mathcal{E}, \mathbb{H}) := \inf_{\substack{V_n \subset \mathbb{H} \\ \dim V_n = n}} \sup_{z \in \mathbf{Z}} \|u(z) - P_{V_n} u(z)\|. $$

The $L^2$ Kolmogorov $n$-width of $\mathcal{E}$ is defined by

$$ (2.9) \qquad d_{2,n}(\mathcal{E}, \mathbb{H}) := \inf_{\substack{V_n \subset \mathbb{H} \\ \dim V_n = n}} \left( \int_{z \in \mathbf{Z}} \|u(z) - P_{V_n} u(z)\|^2 dz \right)^{1/2}. $$

Let us point out that it can easily be checked that

$$ (2.10) \qquad d_{2,n}(\mathcal{E}, \mathbb{H}) \le |\mathbf{Z}|^{1/2} d_{\infty,n}(\mathcal{E}, \mathbb{H}) $$

where $|\mathbf{Z}|$ refers to the Lebesgue measure of the set $\mathbf{Z}$.

In the following, we extend this definition in a meaningful way to the case where $\mathcal{E}$ is not a subset of a Hilbert space $\mathbb{H}$, but of a metric space $\mathcal{M}$ equipped with a distance $\delta$. First a natural generalization consists
in replacing the quantity $||u(z) - P_{V_n}u(z)||$ by the quantity $\inf_{v_n \in V_n} \delta(u(z), v_n)$. Note however that there is no notion of vectorial subspace in a metric space. Instead, we consider the notion of barycenters whenever this is well-defined. In a formal way, this corresponds for a given $\mathbf{u} := (u_1, \dots, u_n) \in \mathbb{M}^n$ to generating all possible barycenters defined as

$$ (2.11) \qquad \mathrm{Bar}_{\delta}^{\lambda}(\mathbf{u}) \in \underset{u \in \mathbb{M}}{\operatorname{argmin}} \sum_{i=1}^{n} \lambda_i \delta(u, u_i)^2, $$

for $\lambda = (\lambda_1, \dots, \lambda_n) \in \Lambda_n$ in a subset $\Lambda_n$ of $\mathbb{R}^n$ guaranteeing that problem (2.11) is well-posed. We assume in the following that the metric $\delta$ is such that there exists at least one solution $\mathrm{Bar}_{\delta}^{\lambda}(\mathbf{u})$ to problem (2.11). In the case when $\mathbb{M}$ is a Hilbert space $\mathbb{H}$ and the metric $\delta$ is defined by $\delta(u, v) = ||u - v||$ for all $u, v \in \mathbb{H}$, the solution to problem (2.11) is unique if and only if $\sum_{i=1}^n \lambda_i > 0$ and is given by $\sum_{i=1}^n \lambda_i u_i$ provided $\sum_{i=1}^n \lambda_i = 1$. In general, we will denote in the sequel $\mathcal{B}_{\delta}^{\lambda}(\mathbf{u})$ the set of minimizers to (2.11).

Assume now that $\mathcal{E} \subset \mathbb{M}$. The most straightforward extension of the notion of Kolmogorov $n$-width in a metric space setting is given in the following definition.

**DEFINITION 2.2.** The metric $L^{\infty}$ Kolmogorov $n$-width of the set $\mathcal{E} \subset \mathbb{M}$ is defined by

$$ d_{\infty,n}(\mathcal{E}, \mathbb{M}) = \inf_{\mathbf{m} \in \mathbb{M}^n} \sup_{z \in \mathcal{Z}} \inf_{\lambda \in \Lambda_n} \inf_{b \in \mathcal{B}_{\delta}^{\lambda}(\mathbf{m})} \delta(u(z), b). $$

Similarly, the metric $L^2$ Kolmogorov $n$-width is defined by

$$ d_{2,n}(\mathcal{E}, \mathbb{M}) := \inf_{\mathbf{m} \in \mathbb{M}^n} \left( \int_{z \in \mathcal{Z}} \inf_{\lambda \in \Lambda_n} \inf_{b \in \mathcal{B}_{\delta}^{\lambda}(\mathbf{m})} \delta(u(z), b)^2 \right)^{1/2}. $$

**2.3. Wasserstein spaces.** The aim of this section is to present some preliminaries about Wasserstein and mixture Wasserstein spaces.

**2.3.1. Wasserstein metric.** Let $\mathcal{P}_2(\mathbb{R})$ denote the set of probability measures on $\mathbb{R}$ with finite second-order moments. The 2-Wasserstein distance over $\mathcal{P}_2(\mathbb{R})$ is defined for $u, v \in \mathcal{P}_2(\mathbb{R})$ as

$$ W_2(u, v) := \inf_{\pi \in \Pi(u, v)} \left( \int_{\mathbb{R}^2} (x - y)^2 d\pi(x, y) \right)^{1/2}, $$

where $\Pi(u, v)$ is the set of probability measures over $\mathbb{R}^2$ with marginals $u$ and $v$, which is called the set of transport plans between $u$ and $v$. Note that $(\mathcal{P}_2(\mathbb{R}), W_2)$ is a geodesic metric space [6].

For any measurable map $T: \mathbb{R} \to \mathbb{R}$ and any probability measure $\rho \in \mathcal{P}_2(\mathbb{R})$, the push-forward measure $T_{\#}\rho$ is defined as the probability measure on $\mathbb{R}$ defined so that for all measurable sets $B \subset \mathbb{R}$, $T_{\#}\rho(B) = \rho(T^{-1}(B))$.

For $n \in \mathbb{N}^*$, let us denote by

$$ \Lambda_n = \left\{ (\lambda_1, \dots, \lambda_n) \in (\mathbb{R}_+)^n, \sum_{i=1}^n \lambda_i = 1 \right\} $$

the set of barycentric weights of cardinality $n$. The Wasserstein barycenter of a collection of $n$ probability measures $\mathbf{u} := (u_1, \dots, u_n) \in \mathcal{P}_2(\mathbb{R})^n$ associated to a set of barycentric weights $\lambda := (\lambda_1, \dots, \lambda_n) \in \Lambda_n$ is then defined (see [1]) as the unique solution to the problem

$$ (2.12) \qquad \inf_{u \in \mathcal{P}_2(\mathbb{R})} \sum_{i=1}^{n} \lambda_i W_2(u, u_i)^2. $$

The unique minimizer of (2.12) is denoted by $\mathrm{Bar}_{W_2}^{\lambda}(\mathbf{u})$. This barycenter is also related to the so-called
multi-marginal optimal transport problem [1, 12], defined, given $n$ elements $\mathbf{u} = (u_1, \dots, u_n)$ in $\mathcal{P}_2(\mathbb{R})^n$, as

$$ (2.13) \qquad mW_2(\mathbf{u}; \boldsymbol{\lambda}) := \inf_{\pi \in \Pi(u_1, \dots, u_n)} \left( \int_{\mathbb{R}^n} \frac{1}{2} \sum_{i,j=1}^n \lambda_i \lambda_j (x_i - x_j)^2 d\pi(x_1, \dots, x_n) \right)^{1/2}, $$

where $\Pi(u_1, \dots, u_n)$ is the set of probability measures over $\mathbb{R}^n$ with marginals $u_1, \dots, u_n$, and there holds

$$ mW_2(\mathbf{u}; \boldsymbol{\lambda}) = \left[ \sum_{i=1}^{n} \lambda_i W_2(\mathrm{Bar}_{W_2}^{\boldsymbol{\lambda}}(\mathbf{u}), u_i) \right]^{1/2}. $$

In the present one-dimensional setting, the Wasserstein distance and barycenter can be expressed in a more direct way using the inverse cumulative distribution function of the considered probability measures. More precisely, we introduce the cumulative distribution function (cdf) of an element $u \in \mathcal{P}_2(\mathbb{R})$ as

$$ \mathrm{cdf}_u : x \in \mathbb{R} \mapsto \int_{-\infty}^{x} du(t) $$

and its inverse cumulative distribution function (icdf) as the generalized inverse of the cdf:

$$ \mathrm{icdf}_u : \begin{cases} [0, 1] & \to \mathbb{R}, \\ s & \mapsto \mathrm{cdf}_u^{-1} := \inf\{x \in \mathbb{R}, \mathrm{cdf}_u(x) > s\}. \end{cases} $$

Then, for any $(u, v) \in \mathcal{P}_2(\mathbb{R})^2$, there holds

$$ (2.14) \qquad W_2(u, v) = \|\mathrm{icdf}_u - \mathrm{icdf}_v\|_{L^2([0,1])}, $$

and for any set of barycentric weights $\boldsymbol{\lambda} := (\lambda_1, \dots, \lambda_n) \in \Lambda_n$ and $\mathbf{u} := (u_1, \dots, u_n) \in \mathcal{P}_2(\mathbb{R})^n$, the icdf of the barycenter $\mathrm{Bar}_{W_2}^{\boldsymbol{\lambda}}(\mathbf{u})$ satisfies

$$ (2.15) \qquad \mathrm{icdf}_{\mathrm{Bar}_{W_2}^{\boldsymbol{\lambda}}}(\mathbf{u}) = \sum_{i=1}^{n} \lambda_i \mathrm{icdf}_{u_i}. $$

We will significantly use this convenient property (2.15), which is specific to the one-dimensional setting, in our analysis.

### 2.3.2. Slater mixture Wasserstein metric.
In higher dimensional spaces, such a characterization does not exist, so that the computation of Wasserstein distances and barycenters is more involved. However, for some specific classes of probability distributions, the Wasserstein distances and barycenters are explicit. This is the case for Gaussian distributions, and more generally, for all location-scatter distributions [2], i.e. all distributions that can be related with an affine transportation map, see also [9, Section 4.1]. In this contribution, we will draw a particular interest to Slater distributions, as defined in (2.3). Noting that the mean of the Slater distribution $S_{\zeta,r}$ is $r$ and the variance is $2/\zeta^2$, we can easily obtain the explicit expression of the Wasserstein distance between two Slater distributions from [2, Theorem 2.3]. More precisely, for $\zeta_1, \zeta_2 > 0$ and $r_1, r_2 \in \mathbb{R}$, there holds

$$ (2.16) \qquad W_2(S_{\zeta_1, r_1}, S_{\zeta_2, r_2})^2 = (r_1 - r_2)^2 + 2 \left( \frac{1}{\zeta_1} - \frac{1}{\zeta_2} \right)^2. $$

Moreover, thanks to [2, Theorem 2.4], it is also possible to obtain an explicit expression for the barycenter between $n \in \mathbb{N}^*$ Slater distributions. Indeed, for $\boldsymbol{\lambda} := (\lambda_1, \dots, \lambda_n) \in \Lambda_n$, $\zeta := (\zeta_1, \dots, \zeta_n) \in (\mathbb{R}_+^*)^n$, $\mathbf{r} := (r_1, \dots, r_n) \in \mathbb{R}^n$, denoting by $S_{\zeta,r} := (S_{\zeta_1, r_1}, \dots, S_{\zeta_n, r_n})$,

$$ (2.17) \qquad \mathrm{Bar}_{W_2}^{\boldsymbol{\lambda}}(S_{\zeta,r}) = S_{\zeta,r}^{\boldsymbol{\lambda}}, $$
with

$$ (2.18) \qquad \zeta^{\lambda} = \left[ \sum_{i=1}^{n} \frac{\lambda_i}{\zeta_i} \right]^{-1} \quad \text{and} \quad r^{\lambda} = \sum_{i=1}^{n} \lambda_i r_i. $$

Since the solutions of problem (2.1) are convex combinations of Slater distributions, as stated in (2.2), we are interested in approaches to efficiently compute Wasserstein-like distances and barycenters for such objects. To this aim, let us start by precisely defining mixtures of Slater distributions. A mixture of Slater distributions $m$ is a finite convex combination of Slater distributions i.e. a probability distribution such that there exists $K \in \mathbb{N}^*$, a $K$-tuple $\zeta = (\zeta_1, \dots, \zeta_K) \in (\mathbb{R}_+^*)^K$, a position tuple $\mathbf{r} = (r_1, \dots, r_K) \in \mathbb{R}^K$ and barycentric weights $\pi = (\pi_1, \dots, \pi_K) \in \Lambda_K$ such that

$$ m = \sum_{k=1}^{K} \pi_k S_{\zeta_k, r_k}. $$

Such a probability distribution will be denoted by $m_{\zeta,r,\pi}$ in the following. Therefore, the solution of (2.1) with parameters $M$, $\mathbf{r}$, and $\mathbf{z}$ is a mixture of Slater distributions with $K = M$, position parameters $\mathbf{r}$, scale parameters $\zeta_1 = \dots = \zeta_K = \zeta_{\mathbf{r},\mathbf{z}}$ and weights $\pi = \pi^{\mathbf{r},\mathbf{z}}$.

We denote by SM($\mathbb{R}$) the subset of $\mathcal{P}_2(\mathbb{R})$ of Slater mixtures, that is

$$ (2.19) \qquad \text{SM}(\mathbb{R}) := \left\{ \begin{array}{l} m \in \mathcal{P}_2(\mathbb{R}), \exists K \in \mathbb{N}^*, \exists \zeta_1, \dots, \zeta_K \in \mathbb{R}_+^*, \exists r_1, \dots, r_K \in \mathbb{R}, \\ \exists \pi = (\pi_1, \dots, \pi_K) \in \Lambda_K, m = \sum_{k=1}^{K} \pi_k S_{\zeta_k, r_k} \end{array} \right\}. $$

In order to compute meaningful and easily computable distances, we in fact endow this space SM($\mathbb{R}$) not with the Wasserstein distance, but with a modified Wasserstein distance, as was proposed in [6, 3] for mixtures of Gaussians, and recently extended to the case of Slater functions in [9]. Therefore, we endow SM($\mathbb{R}$) with this Wasserstein-type distance denoted by MW$_2$. More precisely, for $m^1, m^2 \in \text{SM}(\mathbb{R})$, with parameters $K^1, K^2 \in \mathbb{N}^*, \zeta^1 := (\zeta_1^1, \dots, \zeta_{K_1}^1) \in (\mathbb{R}_+^*)^{K_1}, \zeta^2 := (\zeta_1^2, \dots, \zeta_{K_2}^2) \in (\mathbb{R}_+^*)^{K_2}, \mathbf{r}^1 := (r_1^1, \dots, r_{K_1}^1) \in \mathbb{R}^{K_1}, \mathbf{r}^2 := (r_1^2, \dots, r_{K_2}^2) \in \mathbb{R}^{K_2}, \pi^1 = (\pi_1^1, \dots, \pi_{K_1}^1) \in \Lambda_{K^1}, \pi^2 = (\pi_1^2, \dots, \pi_{K_2}^2) \in \Lambda_{K^2}$, namely

$$ m^1 = m_{\zeta^1, r^1, \pi^1} = \sum_{k^1=1}^{K^1} \pi_{k^1}^1 m_{k^1}^1, \quad m^2 = m_{\zeta^2, r^2, \pi^2} = \sum_{k^2=1}^{K^2} \pi_{k^2}^2 m_{k^2}^2, $$

where

$$ \text{for } i = 1, 2, \quad \text{for } k^i = 1, \dots, K^i, \quad m_{k^i}^i = S_{\zeta_{k^i}^i, r_{k^i}^i}, $$

then the modified Wasserstein distance on the set of Slater mixtures is defined through the following minimization problem

$$ (2.20) \qquad \text{MW}_2(m^1, m^2) := \left[ \min_{w \in \Pi(\pi^1, \pi^2)} \sum_{k^1=1}^{K^1} \sum_{k^2=1}^{K^2} w_{k^1, k^2} W_2^2(m_{k^1}^1, m_{k^2}^2) \right]^{1/2}, $$

where

$$ \Pi(\pi^1, \pi^2) = \left\{ w \in \mathbb{R}_+^{K^1 \times K^2}, \sum_{k^1=1}^{K^1} w_{k^1,k^2} = \pi_{k^2}^2, \sum_{k^2=1}^{K^2} w_{k^1,k^2} = \pi_{k^1}^1 \right\}. $$

The minimization problem (2.20) is well-posed as the function to minimize is continuous on the bounded closed finite dimensional set $\Pi(\pi^1, \pi^2)$. Note that the Wasserstein distances $W_2^2(m_{k^1}^1, m_{k^2}^2)$ can be analytically computed using formula (2.16). In [9], it is proved that (SM($\mathbb{R}$), MW$_2$) is a geodesic metric space.
A multi-marginal optimal transport problem can similarly be defined for the distance MW₂. More precisely, for all $\lambda = (\lambda_1, \dots, \lambda_n) \in \Lambda_n$ and all $\mathbf{m} = (m^1, \dots, m^n) \in \mathrm{SM}(\mathbb{R})^n$, we introduce

$$ (2.21) \qquad m\mathrm{MW}_2(\mathbf{m}; \boldsymbol{\lambda}) := \left[ \min_{w \in \Pi(\boldsymbol{\pi}^1, \dots, \boldsymbol{\pi}^n)} \sum_{\mathbf{k} \in \mathbf{K}} w_{\mathbf{k}} mW_2^2(m_{k1}^1, \dots, m_{kn}^n; \boldsymbol{\lambda}) \right]^{1/2}, $$

where $\mathbf{k} = (k^1, \dots, k^n)$ is in $\mathbf{K} = \{1, \dots, K^1\} \times \dots \times \{1, \dots, K^n\}$ with

$$ \forall i = 1, \dots, n, \quad m^i = \sum_{k^i=1}^{K^i} \pi_{k^i}^i m_{k^i}^i, $$

so that $(w_{\mathbf{k}})_{\mathbf{k} \in \mathbf{K}}$ is an $n$-order tensor with respective dimensions $(K^1, \dots, K^n)$. Recall that the quantities $mW_2^2(m_{k1}^1, \dots, m_{kn}^n; \boldsymbol{\lambda})$ are defined using formula (2.13). Moreover, $\Pi(\boldsymbol{\pi}^1, \dots, \boldsymbol{\pi}^n)$ is the set of such tensors with non-negative coefficients satisfying the constraints

$$ \forall i = 1, \dots, n, \quad \sum_{k^1, \dots, k^{i-1}, k^{i+1}, \dots, k^n} w_{k^1, \dots, k^n} = \pi_{k^i}^i, $$

where $\boldsymbol{\pi}^i = (\pi_{k^i}^i)_{1 \le k^i \le K^i}$.

A barycenter between $n$ Slater mixtures $\mathbf{m} = (m^1, \dots, m^n) \in \mathrm{SM}(\mathbb{R})^n$ with barycentric weights $\boldsymbol{\lambda} = (\lambda_1, \dots, \lambda_n) \in \Lambda_n$ for the mixture distance MW₂ (see (2.11)) is defined as an element $\mathrm{Bar}_{\mathrm{MW}_2}^{\boldsymbol{\lambda}}(\mathbf{m}) \in \mathrm{SM}(\mathbb{R})$ solution to the following minimization problem:

$$ (2.22) \qquad \mathrm{Bar}_{\mathrm{MW}_2}^{\boldsymbol{\lambda}}(\mathbf{m}) \in \underset{b \in \mathrm{SM}(\mathbb{R})}{\operatorname{argmin}} \sum_{i=1}^{n} \lambda_i \mathrm{MW}_2(b, m^i)^2. $$

It is shown in [9] that there always exists at least one solution to (2.22) and that any solution can be expressed as

$$ (2.23) \qquad \mathrm{Bar}_{\mathrm{MW}_2}^{\boldsymbol{\lambda}}(\mathbf{m}) = \sum_{\mathbf{k} \in \mathbf{K}} (w^*(\mathbf{m}; \boldsymbol{\lambda}))_{\mathbf{k}} \mathrm{Bar}_{W_2}^{\boldsymbol{\lambda}}(m_{k1}^1, \dots, m_{kn}^n), $$

where $w^*(\mathbf{m}; \boldsymbol{\lambda})$ is a solution to the minimization problem (2.21).

REMARK 2.3. Note that all minimizers $w^*(\mathbf{m}; \boldsymbol{\lambda})$ of problem (2.21) have at most $K_1+K_2+\dots+K_n-n+1$ nonzero components, see e.g. [6].

### 3. Theoretical estimations of the decay of Kolmogorov n-widths.
In this section, we provide estimates for metric Kolmogorov n-width of sets of solutions of the one-dimensional toy model presented in Section 2.1 using three different metrics. We state our main results here and postpone their proofs to Section 6. First, we study the Kolmogorov n-width with an underlying $L^2$-norm, which is the standard choice in linear reduced order modelling. We then turn to the metric Kolmogorov n-width with respect to the Wasserstein metric, which in the one-dimensional case, can be recast as a Kolmogorov n-width for the $L^2$-norm on the inverse cumulative distribution functions of the solutions, and we show that the decay rate is faster than for the traditional $L^2$-norm. Finally, we show that this decay can be improved by considering the modified Wasserstein metric defined in (2.20).

#### 3.1. On linear approximations: case $\mathbb{H} = L^2(\mathbb{R})$.
In this section, we state estimates on the $L^2$ Kolmogorov n-width decay of a set of solutions of (2.1) where the parameters vary in a compact set. We consider the charges $\mathbf{z}$ to be fixed and define the set of solutions

$$ \mathcal{M}_{\mathbf{z}}^{R} := \{u_{\mathbf{r},\mathbf{z}}, \mathbf{r} \in [-R, R]^M\}, $$

where $R \in \mathbb{R}_+$ is a given parameter and $u_{\mathbf{r},\mathbf{z}}$ is the solution to (2.1). We first consider the case $M=1$, i.e. the case where the potential is a single Dirac delta. We recall in (2.5) the particular form of $\zeta_{\mathbf{r},\mathbf{z}}$ and denote by $u_r = u_{\mathbf{r},\mathbf{z}}$ the solution of (2.1). We also denote $\mathcal{M}_{\mathbf{z}}^R$ the above set in this case.
**THEOREM 3.1.** There exist positive constants $c_R$, $C_R$, $\tilde{c}_R$ and $\tilde{C}_R$ depending on $R$ such that for all $n \in \mathbb{N}^*$,

$$ (3.1) \qquad c_R n^{-\frac{3}{2}} \le d_{\infty,n}(\mathcal{M}_z^R, L^2(\mathbb{R})) \le C_R n^{-\frac{3}{2}}, $$
and

$$ (3.2) \qquad \tilde{c}_R n^{-\frac{3}{2}} \le d_{2,n}(\mathcal{M}_z^R, L^2(\mathbb{R})) \le \tilde{C}_R n^{-\frac{3}{2}}. $$

The proof of Theorem 3.1 is presented in Section 6.2. Next, we claim that a similar result holds true for any system of several fixed charges $\mathbf{z}$.

**COROLLARY 3.2.** The following lower bound

$$ (3.3) \qquad c_R n^{-\frac{3}{2}} \le d_{\infty,n}(\mathcal{M}_z^R, L^2(\mathbb{R})) $$

holds with the positive constant $c_R$ from (3.1).

*Proof.* Remark that for $\mathbf{r} = (r, \dots, r)$, the problem (2.1) corresponds to $M = 1$, with a unique charge $z = \sum_{m=1}^{M} z_m$. Thus, $d_{\infty,n}(\mathcal{M}_z^R, L^2(\mathbb{R})) \le d_{\infty,n}(\mathcal{M}_z^R, L^2(\mathbb{R}))$ and (3.3) follows from (3.1). $\square$

**REMARK 3.3.** It does not seem trivial to obtain a similar upper bound as the lower bound in (3.3), nor the bounds for $d_{2,n}(\mathcal{M}_z^R, L^2(\mathbb{R}))$. It is not clear either whether the bound (3.3) is optimal. However, since our point is to show that the Kolmogorov $n$-width of the set of solutions in $L^2(\mathbb{R})$ is larger than the generalized Kolmogorov $n$-width using Wasserstein metric, the lower bound (even suboptimal) is sufficient here.

**3.2. On the Wasserstein transport metric: case $M = \mathcal{P}_2(\mathbb{R})$ with $\delta = W_2$.** Let us now give some insight on what happens when we consider the Wasserstein distance $\delta = W_2$ with $M = \mathcal{P}_2(\mathbb{R})$ and comment on the theoretical results.

**3.2.1. Preliminary discussion.** It holds from (2.14) that barycenters can be expressed using the icdf of $m_1, \dots, m_n$. there holds for any $m \in \mathcal{P}_2(\mathbb{R})$ and $\mathbf{m} = (m_1, \dots, m_n) \in \mathcal{P}_2(\mathbb{R})^n$

$$ W_2(m, \text{Bar}_{W_2}^{\lambda}(\mathbf{m})) = \left\| \text{icdf}_m - \sum_{i=1}^{n} \lambda_i \text{icdf}_{m_i} \right\|_{L^2([0,1])}. $$

It is then natural to relate the generalized Kolmogorov $n$-width of a set $\mathcal{E} \subset \mathcal{P}_2(\mathbb{R})$ associated with the Wasserstein metric with the Kolmogorov $n$-width of the set $\mathcal{T} = \{\text{icdf}_{u(z)}, z \in \mathbb{Z}\}$ in the Hilbert space $L^2(0,1)$. It can then be easily checked that

$$ d_{\infty,n}(\mathcal{E}, \mathcal{P}_2(\mathbb{R})) \ge d_{\infty,n}(\mathcal{T}, L^2([0,1))), $$

and similarly that

$$ d_{2,n}(\mathcal{E}, \mathcal{P}_2(\mathbb{R})) \ge d_{2,n}(\mathcal{T}, L^2([0,1))). $$

Proving decay estimates on the quantities $d_{\infty,n}(\mathcal{E}, \mathcal{P}_2(\mathbb{R}))$ and $d_{2,n}(\mathcal{E}, \mathcal{P}_2(\mathbb{R}))$ appears to be a difficult task. In the present work, we manage to prove decay estimates with respect to $n$ of $d_{\infty,n}(\mathcal{T}, L^2([0,1]))$ and $d_{2,n}(\mathcal{T}, L^2([0,1]))$. We first would like to point out that the latter quantities are also the ones for which decay estimates have been proven for some conservative transport equations in [10].

**3.2.2. Main results.** In this section, we give estimates of $L^2(0,1)$ Kolmogorov $n$-widths decays of the set of icdfs of solutions, and show that they converge faster to zero as a function of $n$ than the $L^2(\mathbb{R})$ Kolmogorov $n$-widths of the original solution set. First, in the case $M = 1$ with a fixed charge $z$, we have the following proposition.

**PROPOSITION 3.4.** The $L^2(0,1)$ Kolmogorov $n$-width of the set of icdf of solutions $\{\text{icdf}_{u_r}, r \in \mathbb{R}\}$ with an unbounded position parameter $r \in \mathbb{R}$ is equal to zero for $n > 1$.
Proof. The solutions $u_r$ are translations one to another, hence the set of icdf of solutions is a subset of the 2-dimensional vector space span${\mathrm{icdf}_{u_0}, 1_{(0,1)}}$.

In the rest of Section 3.2, we consider symmetric systems with $M = 2$, i.e. $\mathbf{r} = (-r, r)$ and $\mathbf{z} = (z, z)$. For simplicity of notation, we denote by $u_r = u_{\mathbf{r},z}$ the ground state of the above system and $\zeta_r = \zeta_{\mathbf{r},z}$ which is given by (2.6).

**Transport approximation for** $r \le R$. We first consider the set of icdf of solutions where the position parameter $r$ is bounded by some positive $R > 0$. Let us introduce the set

$$ \mathcal{T}_R^- := \{\mathrm{icdf}_{u_r}, r \in [0, R]\} \subset L^2(0, 1). $$

We show a better Kolmogorov $n$-width decay than its counterpart for the original solution set.

**THEOREM 3.5.** There exists a constant $C > 0$ independent of $R$ such that for all $n \in \mathbb{N}^*$,

$$ d_{\infty,n}(\mathcal{T}_R^-, L^2(0,1)) \le C e^{5\zeta_R R} n^{-\frac{5}{2}}. $$

The proof of Theorem 3.5 is postponed to Section 6.3.

**Transport approximation for** $r \ge R$. We now state an asymptotic result that holds for position parameters $r \ge R$. Introducing the set

$$ \mathcal{T}_R^+ := \{\mathrm{icdf}_{u_r}, r \ge R\} \subset L^2(0, 1), $$

we have the following result.

**THEOREM 3.6.** There exists a constant $C > 0$ independent of $R$ such that for all $n \in \mathbb{N}^*$,

$$ d_{\infty,4}(\mathcal{T}_R^+, L^2(0,1)) \le C R e^{-\frac{1}{2}\zeta_R R}. $$

The proof of Theorem 3.6 is postponed to Section 6.4.

**Transport approximation for** $r \in \mathbb{R}^+$. Now consider the unbounded set of inverse cumulative distribution function of solutions

$$ \mathcal{T} := \{\mathrm{icdf}_{u_r}, r \in \mathbb{R}^+\}. $$

We finally state our final result on the Kolmogorov $n$-width decay of $\mathcal{T}$.

**THEOREM 3.7.** For all $\varepsilon > 0$, there exists a constant $C_\varepsilon > 0$ such that for all $n \in \mathbb{N}^*$,

$$ d_{2,n}(\mathcal{T}, L^2(0,1)) \le C_\varepsilon n^{\frac{25}{21}-\frac{\varepsilon}{2}-\frac{5}{2}}. $$

The proof of Theorem 3.7 is postponed to Section 6.5.

**3.3. On the mixture Wasserstein transport metric: case $M = SM(\mathbb{R})$ with $\delta = MW_2$.** We now provide theoretical results related to this mixture distance and comment on the theoretical results.

**3.3.1. Preliminary discussion.** We first specify the definitions of metric Kolmogorov widths considered in the mixture Wasserstein framework of our theoretical analysis.

First the computation of an exact mixture Wasserstein barycenter using formula (2.23) requires, for any collection of barycentric weights $\lambda \in \Lambda_n$ and of Slater mixtures $m \in SM(\mathbb{R})^n$, the computation of a solution $w^*(m; \lambda)$ to the minimization problem (2.21). In our model order reduction context, this appears to be too costly from a computational point of view. For this reason we rather rely on the use of approximate mixture barycenters, the computation and definition of which do not require the resolution of multiple problems of
the form (2.21). More precisely, denoting by $\bar{\lambda} := (1/n, 1/n, \dots, 1/n) \in \Lambda_n$, for any (possibly nonunique) solution $w^*(m; \bar{\lambda})$ of (2.21), we define (recall formula (2.23) for the exact mixture Wasserstein barycenter)

$$ (3.4) \qquad \overline{\mathrm{Bar}}_{\mathrm{MW}_2}^{\lambda}(m) := \sum_{\mathbf{k} \in \mathcal{K}} (w^*(m; \bar{\lambda}))_{\mathbf{k}} \mathrm{Bar}_{\mathrm{W}_2}^{\lambda}(m_{k^1}^1, \dots, m_{k^n}^n). $$

Note that it is easy to check that the approximate mixture barycenter remains an interpolation between the mixtures $m^1, \dots, m^n$ in the following sense: if $\lambda = (\delta_{n_0,i})_{1 \le i \le n}$ for some $1 \le n_0 \le n$, then $\overline{\mathrm{Bar}}_{\mathrm{MW}_2}^{\lambda}(m) = m^{n_0}$.

In addition, we observed in our numerical tests that restricting the set of weights $\lambda$ to the set $\Lambda_n$ of barycentric weights is quite restrictive. Additionally, the definition of approximate mixture barycenters can easily be extended to the case where the collection of weights $\lambda$ belongs to the larger set

$$ (3.5) \qquad \Omega_n(m) = \left\{ \lambda \in \mathbb{R}^n, \forall \mathbf{k} \in \{1, \dots, K^1\} \times \dots \times \{1, \dots, K^n\}, \sum_{i=1}^n \frac{\lambda_i}{\zeta_{k^i}^i} > 0 \right\}. $$

Indeed, this requires to extend the definition of Wasserstein barycenters between a collection of Slater distributions, which we do via (2.17) using (2.18), since the quantities defined in (2.18) are still well-defined for $\lambda \in \Omega_n(m)$. In the following, for all $\lambda \in \Omega_n(m)$, we denote by $\overline{\mathcal{B}}_{\mathrm{MW}_2}^{\lambda}(m)$ the set of approximate mixture Wasserstein barycenters $\overline{\mathrm{Bar}}_{\mathrm{MW}_2}^{\lambda}(m)$ which can be written under the form (3.4).

We are now in position to consider two following quantities, which are related to (but different from) the metric Kolmogorov widths in the mixture Wasserstein space. We thus define the extended approximate metric $L^\infty$ Kolmogorov width

$$ (3.6) \qquad \tilde{d}_{\infty,n}(\mathcal{E}, SM(\mathbb{R})) = \inf_{m \in SM(\mathbb{R})^n} \sup_{z \in \mathcal{Z}} \inf_{\lambda \in \Omega_n(m)} \inf_{b \in \overline{\mathcal{B}}_{\mathrm{MW}_2}^{\lambda}(m)} \mathrm{MW}_2(u(z), b). $$

and the extended approximate metric $L^2$ Kolmogorov width

$$ \tilde{d}_{2,n}(\mathcal{E}, SM(\mathbb{R})) := \inf_{m \in SM(\mathbb{R})^n} \left( \int_{z \in \mathcal{Z}} \inf_{\lambda \in \Omega_n(m)} \inf_{b \in \overline{\mathcal{B}}_{\mathrm{MW}_2}^{\lambda}(m)} \mathrm{MW}_2(u(z), b)^2 \right)^{1/2}. $$

**3.3.2. Main result.** In this part, we also focus on the $M = 2$ symmetric case, i.e. where $\mathbf{r} = (-r, r)$ and with fixed equal charges $\mathbf{z} = (z, z)$. Once again, we denote by $\mathcal{M}_z$ the set of solutions. In that case, the extended approximate metric Kolmogorov $n$-width relative to the mixture distance $\mathrm{MW}_2$ (2.20) defined in (3.6) is zero for $n \ge 2$. This means that in that case, any solution can be exactly obtained from a collection only two Slater mixtures.

**THEOREM 3.8.** Let $SM(\mathbb{R}) \subset \mathcal{P}_2(\mathbb{R})$ be the metric space of mixtures of Slater distributions defined in (2.19) endowed with its metric $\mathrm{MW}_2$ defined in (2.20). For all $n > 1$, there holds

$$ \tilde{d}_{\infty,n}(\mathcal{M}_z, SM(\mathbb{R})) = \tilde{d}_{2,n}(\mathcal{M}_z, SM(\mathbb{R})) = 0. $$

The proof of Theorem 3.8 is postponed to Section 6.6.

**3.4. Link with other types of Kolmogorov widths.** We now compare the proposed extended metric Kolmogorov widths with two other notions of generalized Kolmogorov widths.

**3.4.1. Link with nonlinear Kolmogorov widths.** Let us point out that the metric Kolmogorov width $d_{\infty,n}(\mathcal{E}, \mathcal{P}_2(\mathbb{R}))$ can be related to nonlinear Kolmogorov widths [25, 7]. Indeed this quantity can be rewritten as

$$ (3.7) \qquad \inf_{m \in \mathcal{M}^n} \sup_{z \in \mathcal{Z}} \delta(u(z), D_m \circ C_m(u(z))), $$
with $M = \mathcal{P}_2(\mathbb{R})$, $\delta = W_2$, and where for all $\mathbf{m} \in \mathbb{M}^n$, $D_m$ refers to a particular decoder map and $C_m$ a particular encoder map defined as

$$D_m : \Lambda_n \ni \lambda \mapsto \mathrm{Bar}_{W^2}^{\lambda}(\mathbf{m}) \in \mathcal{P}_2(\mathbb{R}) \quad \text{and} \quad C_m : \mathcal{P}_2(\mathbb{R}) \ni u \mapsto \underset{\lambda \in \Lambda_n}{\operatorname{argmin}} W_2(u, \mathrm{Bar}_{W^2}^{\lambda}(\mathbf{m})) \in \Lambda_n.$$

However, in the case of $\tilde{d}_{\infty,n}(\mathcal{E}, \mathrm{SM}(\mathbb{R}))$, the situation is more subtle since there is in general no uniqueness of approximate mixture Wasserstein barycenters. It is therefore not obvious to rewrite $\tilde{d}_{\infty,n}(\mathcal{E}, \mathrm{SM}(\mathbb{R}))$ under the form (3.7) with $M = \mathrm{SM}(\mathbb{R})$, $\delta = MW_2$ for some decoder maps $D_m$ and encoder maps $C_m$ parametrized by $\mathbf{m} \in \mathbb{M}^n$.

**3.4.2. Link with Kolmogorov (n, m)- width.** The extended metric Kolmogorov width can also be linked to the Kolmogorov (n, m)-width introduced in [22] as follows. Denoting by $\mathbb{T}$ the set of bijective maps from $\mathbb{R}$ to $\mathbb{R}$, the Kolmogorov (n, m)-width then reads as:

$$d_{n,m}(\mathcal{E}, \mathbb{H}) := \inf_{\substack{V_n \subset \mathbb{H} \\ \dim W_m = n \\ \dim \mathrm{Span}\{W_m\} = m}} \sup_{z \in \mathbb{Z}} \inf_{T_m \in \mathrm{Span}\{W_m\} \cap \mathbb{T}} \|u(z) - v_n \circ T_m^{-1}\|.$$

In the present setting, it is more convenient to measure errors with respect to the Wasserstein metric. This is our motivation for introducing the quantity

$$d_{1,m,W_2}(\mathcal{E}, \mathcal{P}_2(\mathbb{R})) := \inf_{\substack{v \in \mathcal{P}_2(\mathbb{R}) \\ W_m \subset \mathbb{T} \\ \dim \mathrm{Span}\{W_m\} = m}} \sup_{z \in \mathbb{Z}} \inf_{T_m \in \mathrm{Span}\{W_m\} \cap \mathbb{T}W_2(u(z), v \circ T_m^{-1})} \{W_m\}.$$

In the case when $\mathbf{m} = (m_1, \dots, m_n) \in \mathcal{P}_2(\mathbb{R})$ is a family of $n$ probability measures which are absolutely continuous with respect to the Lebesgue measure and the support of which is equal to the whole space $\mathbb{R}$, then it holds that, for all $\lambda = (\lambda_1, \dots, \lambda_n) \in \Lambda_n$,

$$\mathrm{Bar}_{W_2}^{\lambda}(\mathbf{m}) = \left( \sum_{i=1}^{n} \lambda_i T_i \right) \# m_1,$$

where for all $1 \le i \le n$, $T_i$ denotes the optimal transport map between $m_1$ and $m_i$. Thus, in the case of the exact Wasserstein metric, if the probability measures are regular enough and have full support in $\mathbb{R}$, $d_{\infty,n}(\mathcal{E}, \mathcal{P}_2(\mathbb{R}))$ is equal to $d_{1,n,W_2}(\mathcal{E}, \mathcal{P}_2(\mathbb{R}))$. However, when the probability measures are not regular enough, or in the case of the mixture Wasserstein metric, the existence of optimal transport maps is no longer guaranteed, which makes difficult to draw a link between the Kolmogorov width considered in the present work, and the (n, m)-widths of [22].

**4. Nonlinear reduced basis method.** Motivated by the fast decay of the generalized Kolmogorov n-width for the Wasserstein mixture distance, we now propose a nonlinear reduced basis method for this problem. The method, as is common in reduced order modeling, is based on an offline phase followed by an online phase. In the offline phase a few representative snapshots are selected thanks to a greedy algorithm. In the online phase, for any new set of parameters, the energy of the system is minimized over the set of barycenters of selected snapshots, using a quasi-Newton minimization algorithm started at several initial points. Let us emphasize here that, despite the fact that the theoretical results stated in Section 3 only hold for the one-dimensional toy problem introduced in Section 2.1, the numerical strategy presented in this section can be used to build nonlinear reduced order models for more general problems.

**4.1. Greedy algorithm.** We first present the greedy algorithm used in the offline phase. Let $\mathbf{z} \in (\mathbb{R}_+^*)^M$ be fixed positive charges, $\mathcal{M}_{\mathbf{z}}^I = \{u_{\mathbf{r},\mathbf{z}}, \mathbf{r} \in I^M\}$ with $I \subset \mathbb{R}$ an interval be a set of solutions of (2.1) and $\mathcal{M}_{tr} \subset \mathcal{M}_{\mathbf{z}}^I$ a finite training set of already computed solutions called snapshots. The aim here is to select the most representative snapshots in $\mathcal{M}_{tr}$, so that any solution $u \in \mathcal{M}_{\mathbf{z}}^I$ can be efficiently approximated with only a few snapshots. The main idea in this greedy algorithm is to select at each iteration the snapshot in $\mathcal{M}_{tr}$ the approximation of which as a mixture barycenter of previously selected snapshots leads to the highest
error. Since the proposed algorithm is generic to any training set where the elements can be represented by mixtures equipped with a mixture distance, we use the notation m for the elements in $M_{tr}$ instead of u.

Recall from Section 2.3 that we write a mixture as $m = \sum_{k=1}^{K} \pi_k m_k$, so for the rest of this section, we denote an element of $M_{tr}$ as $m$. In our case, these mixtures are solutions to problem (2.1), which means that their elements, denoted by $m_k$ are Slater functions with a scale parameter $\zeta$ independent of $k$ and position parameter $r_k$. In the following, mixtures written as $m^i$ where $i$ is an integer follow the same rules of notation, their parameters being $(\pi_{k^i})_{k^i=1}^{K^i}$ for the weights, $\zeta^i$ for the common scale parameters and $(r_{k^i})_{k^i=1}^{K^i}$ for the position parameters. The proposed greedy algorithm is presented below (Algorithm 4.1).

**Algorithm 4.1 GREEDY ALGORITHM**

**Input:** $M_{tr}$, training set; $N$, number of elements to select
Select $m^1$ and $m^2$ solutions to $\argmax_{(m^1,m^2) \in M_{tr}}$ MW$_2(m^1,m^2)^2$.

$$ \mathcal{B} := \{m^1, m^2\} $$

**for** $n = 2, \dots, N-1$ **do**

Select

$$ (4.1) \qquad m^{n+1} \in \argmax_{m \in \mathcal{M}_{tr}} \inf_{\lambda \in \Omega_n(m^1, \dots, m^n)} \mathrm{MW}_2 \left( m, \overline{\mathrm{Bar}}_{\mathrm{MW}_2}^{\lambda} (m^1, \dots, m^n) \right)^2, $$

where

$$ \Omega_n(m^1, \dots, m^n) = \left\{ \lambda \in \mathbb{R}^n, \sum_{i=1}^n \frac{\lambda_i}{\zeta^i} > 0 \right\}. $$

$$ \mathcal{B} = \mathcal{B} \cup \{m^{n+1}\} $$

**end for**

**Output:** Reduced basis $\mathcal{B} \subset M_{tr}$

The keystone of Algorithm 4.1 is the resolution of problem (4.1), and more precisely for $m^1, \dots, m^n \in M_{tr}$ and $m \in M_{tr}$, the resolution of the following minimization problem

$$ (4.2) \qquad \inf_{\lambda \in \Omega_n(m^1, \dots, m^n)} \mathrm{MW}_2 \left( m, \overline{\mathrm{Bar}}_{\mathrm{MW}_2}^{\lambda} (m^1, \dots, m^n) \right)^2. $$

In practice, we start by solving problem (2.21) with $\lambda = \bar{\lambda} = (1/n, \dots, 1/n)$ to obtain the barycenters weights $w^*(m, \bar{\lambda})$ appearing in (3.4) for the calculation of $\overline{\mathrm{Bar}}_{\mathrm{MW}_2}^{\lambda}(m^1, \dots, m^n)$. It can be computed using any linear programming solver as problem (2.21) is a linear problem with linear constraints. In the sequel, for the sake of simplicity, we denote by $w^* := w^*(m, \bar{\lambda})$. Note that in the representation of the barycenters (2.23), one can in fact only consider the indices $k = (k^1, k^2, \dots, k^n) \in \mathbf{K} \subset \{1, \dots, K\} \times \dots \times \{1, \dots, K^n\}$ for which $w_k^*$ is non zero (see Remark 2.3) to reduce the dimensionality of the problem.

Now, by the definition of the distance MW$_2$, we can say that

$$ \mathrm{MW}_2 \left( m, \overline{\mathrm{Bar}}_{\mathrm{MW}_2}^{\lambda} (m^1, \dots, m^n) \right)^2 = \min_{w \in \Pi(\pi, w^*)} \sum_{\mathbf{k} \in \mathbf{K}} \sum_{k=1}^{K} w_{\mathbf{k},k} \mathrm{W}_2 \left( m_k, \overline{\mathrm{Bar}}_{\mathrm{W}_2}^{\lambda} (m_{k^1}^1, \dots, m_{k^n}^n) \right)^2, $$

where the weights $w$ are matrices with non-negative terms indexed by $\mathbf{k} \in \mathbf{K}$ and $k \in \{1, \dots, K\}$ in the set

$$ (4.3) \qquad \Pi(\pi, w^*) := \left\{ w \in (\mathbb{R}_+)^{|\mathbf{K}| \times K}, \forall k \in \{1, \dots, K\}, \sum_{\mathbf{k} \in \mathbf{K}} w_{\mathbf{k},k} = \pi_k, \forall \mathbf{k} \in \mathbf{K}, \sum_{k=1}^{K} w_{\mathbf{k},k} = w_{\mathbf{k}}^* \right\}. $$
Note that the minimization set $\Pi(\pi, w^*)$ does not depend on the parameter $\lambda$ as the previously computed weights of barycenter $w^*$ do not depend on $\lambda$ either. We also have, recalling that the parameters for the mixture $m$ are $r_k$ and $\zeta$

$$W_2(m_k, \text{Bar}_{W_2}^\lambda(m_{k1}^1, \dots, m_{kn}^n)) = \left(r_k - \sum_{i=1}^n \lambda_i r_{ki}^i\right)^2 + 2 \left(\frac{1}{\zeta} - \sum_{i=1}^n \frac{\lambda_i}{\zeta^i}\right)^2 = \lambda^\top A_k \lambda + b_{k,k}^\top \lambda + c_k,$$

where

$$A_k = \mathbf{r}_k^\top \mathbf{r}_k + 2\zeta^\top \zeta, \quad b_{k,k} = -2 \left( r_k \mathbf{r}_k + \frac{2}{\zeta} \zeta \right), \quad \text{and} \quad c_k = r_k^2 + \frac{2}{\zeta^2},$$

with $\mathbf{r}_k = (r_{k1}^1, \dots, r_{kn}^n)$ and $\zeta = (\frac{1}{\zeta^1}, \dots, \frac{1}{\zeta^n})$. In particular the matrices $A_k$ are non-negative for any $k \in \mathbf{K}$. Hence, problem (4.2) reduces to

$$(4.4) \quad \inf_{\lambda \in \Omega_n(m^1, \dots, m^n)} \min_{w \in \Pi(\pi, w^*)} \lambda^\top A \lambda + b_w^\top \lambda + c = \min_{w \in \Pi(\pi, w^*)} \inf_{\lambda \in \Omega_n(m^1, \dots, m^n)} \lambda^\top A \lambda + b_w^\top \lambda + c,$$

where

$$(4.5) \quad A = \sum_{k \in \mathbf{K}} w_k^* A_k, \quad b_w = \sum_{k \in \mathbf{K}} \sum_{k=1}^K w_{k,k} b_{k,k}, \quad \text{and} \quad c = \sum_{k=1}^K \pi_k c_k.$$

The matrix $A$ is also non-negative as a sum of non-negative matrices. Note that since the matrix $A$ does not depend neither on $\lambda$ and on the mixture $m$, we can indeed compute $A$ at each update of $\mathcal{B}$, just like the weights $w^*$. In all our numerical tests, we have numerically checked that the matrix $A$ is in fact positive definite. Then, the solution to the minimization problem

$$\min_{\lambda \in \mathbb{R}^n} \lambda^\top A \lambda + b_w^\top \lambda + c$$

is $\lambda_w = -\frac{1}{2}A^{-1}b_w$, and we can also check here a posteriori that the solution $\lambda_w \in \Omega_n(m)$, so that the solution of $\inf_{\lambda \in \Omega_n(m)} \lambda^\top A \lambda + b_w^\top \lambda + c$ is also $\lambda_w$, which is always the case in the tested examples. If however it turned out not to be the case, since $\Omega_n(m)$ is a convex set, it is possible to directly solve the minimization problem on $\Omega_n$ using quadratic programming. In the case when $A$ might not be invertible (which never occured in the numerical tests performed in this work), the inverse of $A$ would naturally be replaced by its pseudo-inverse in the expression of $\lambda_w$ above, and in the expressions below.

Hence, by putting $\lambda_w$ back in problem (4.4), we have that problem (4.2) is equivalent to

$$\min_{w \in \Pi(\pi, w^*)} -\frac{1}{4} b_w^\top A^{-1} b_w + c = \min_{w \in \Pi(\pi, w^*)} w^\top \left( -\frac{1}{4} B^\top A^{-1} B \right) w + c,$$

by choosing a vectorization for the weights $w$ and where $B$ is such that $Bw = b_w$, which is a concave quadratic minimization problem because the matrix $-\frac{1}{4}B^\top A^{-1}B$ is negative since $A$ is positive. The matrix $-\frac{1}{4}B^\top A^{-1}B$ of this problem has a size $K|\mathbf{K}|$. For a hint on the size $|\mathbf{K}|$, see Remark 2.3. The solution of the problem is in fact a vertex of the polytope $\Pi(\pi, w^*)$, thanks to the convexity of the polytope and the concavity of the problem [11]. We summarize the whole procedure to solve problem (4.1) in Algorithm 4.2 below.

**REMARK 4.1.** In our implementation, we took advantage of the concave setting of the problem and searched the solution directly among the vertices of the polytope $\Pi(\pi, w^*)$ to ensure global optimality. In practice, global optimization packages such as Gurobi could be used.

**REMARK 4.2.** Algorithm 4.1 is in fact more general and can be used wherever the elements of the training set can be represented by mixtures. For example, we can consider a setting where the solutions are Slater
Algorithm 4.2 OFFLINE PROJECTION MINIMIZATION

Input: $\mathcal{B} = \{m^1, \dots, m^n\}$, selected elements

Compute $w^*$ as in problem (2.21).

Compute the matrix $A$ as in (4.5), check that it is positive definite, and compute $A^{-1}$.

Select

$$m^{n+1} \in \underset{m \in \mathcal{M}_{tr}}{\operatorname{argmin}} \min_{w \in \Pi(\pi, w^*)} -\frac{1}{4} b_w^T A^{-1} b_w + c,$$

where $b_w$ and $c$ are given in (4.5).

Output: $m^{n+1}$

mixtures with different scale parameters $\zeta$. In this scenario, the set of admissible weights for the barycenters become a bit more complex and reads

$$\Omega_n(\mathbf{m}) = \left\{ \lambda \in \mathbb{R}^n, \forall \mathbf{k} \in \{1, \dots, K^1\} \times \dots \times \{1, \dots, K^n\}, \sum_{i=1}^{n} \frac{\lambda_i}{\zeta_{k^i}^i} > 0 \right\}.$$

More generally, one can consider any mixtures for which a Wasserstein mixture distance is well-defined, as presented in [9], which includes e.g. Gaussian mixtures, upon modifying the set $\Omega_n(\mathbf{m})$ with the correct parameters of the distributions to ensure admissibility of barycenters.

4.2. Online algorithm. Once the reduced basis is computed, we want to efficiently compute approximations of solutions, given a new position for the nuclei. Since the projection minimization algorithm used in the offline phase requires a high-fidelity estimate of the exact solution, it is not a viable option for the online phase. Here we instead take advantage of the structure of our problem, which is an energy minimization problem (2.4), and we minimize the energy of the new system over the set of barycenters of the elements in the reduced basis. More precisely, assume that we selected $N$ mixtures solutions $m^1, \dots, m^N$ in the offline phase and we want to obtain an approximation to the solution with molecular parameters $\mathbf{r}$. We consider the following optimization problem

$$(4.6) \qquad \inf_{\lambda \in \Omega_N(m^1, \dots, m^N)} \frac{E_{\mathbf{r},\mathbf{z}}\left(\overline{\text{Bar}}_{\text{MW}_2}^{\lambda}(m^1, \dots, m^N)\right)}{\left\|\overline{\text{Bar}}_{\text{MW}_2}^{\lambda}(m^1, \dots, m^N)\right\|_{L^2(\mathbb{R})}^2},$$

where $\Omega_N(m^1, \dots, m^N)$ is the extended set of admissible barycenters, as in the greedy algorithm. For clarity, we write from now on $E_{\mathbf{r},\mathbf{z}}(\lambda)$ instead of $\frac{E_{\mathbf{r},\mathbf{z}}\left(\overline{\text{Bar}}_{\text{MW}_2}^{\lambda}(m^1, \dots, m^N)\right)}{\left\|\overline{\text{Bar}}_{\text{MW}_2}^{\lambda}(m^1, \dots, m^N)\right\|_{L^2(\mathbb{R})}^2}$. Note that the energy functional can be easily computed with the following formula

$$(4.7) \qquad E_{\mathbf{r},\mathbf{z}}(\lambda) = \left[ \sum_{\mathbf{k} \in \mathcal{K}} \sum_{\mathbf{l} \in \mathcal{K}} w_{\mathbf{k}}^* w_{\mathbf{l}}^* (1 + \zeta^{\lambda} |r_{\mathbf{k}}^{\lambda} - r_{\mathbf{l}}^{\lambda}|) e^{-\zeta^{\lambda} |r_{\mathbf{k}}^{\lambda} - r_{\mathbf{l}}^{\lambda}|} \right]^{-1} \cdot \left( \frac{\zeta^{\lambda}^2}{2} \sum_{\mathbf{k} \in \mathcal{K}} \sum_{\mathbf{l} \in \mathcal{K}} w_{\mathbf{k}}^* w_{\mathbf{l}}^* (1 - \zeta^{\lambda} |r_{\mathbf{k}}^{\lambda} - r_{\mathbf{l}}^{\lambda}|) e^{-\zeta^{\lambda} |r_{\mathbf{k}}^{\lambda} - r_{\mathbf{l}}^{\lambda}|} - \zeta^{\lambda} \sum_{m=1}^{M} z_m \sum_{\mathbf{k} \in \mathcal{K}} w_{\mathbf{k}}^* e^{-\zeta^{\lambda} |r_{\mathbf{k}}^{\lambda} - r_m|} \right)$$

where $\zeta^{\lambda} = \left[ \sum_{i=1}^{N} \frac{\lambda_i}{\zeta_i^i} \right]^{-1}$ is the scale parameter and $r_{\mathbf{k}}^{\lambda} = \sum_{i=1}^{N} \lambda_i r_{k^i}^i$ is the position parameter of a Slater component of a barycenter with weights $\lambda$, and $w^*$ is the solution to problem (2.21) for the reduced basis, and can be computed offline.

REMARK 4.3. Let us point out that the minimization problem (4.6) is in general non-convex. Numerical solutions obtained by the minimization algorithms used here then depends on the initialization.
<figure><img src="image_2.png" /><figcaption>FIG. 5.1. Three example solutions in $M_{tr}$.</figcaption></figure>

Since the energy functional as a function of $\lambda$ is nonconvex, and in practice exhibits many local minima, solving problem (4.6) requires to use a global optimization algorithm, preferably very robust to ensure that the global minimizer is found, to guarantee repeatability of the results. The natural optimization procedure is detailed here, and detailed in Algorithm 4.3. We use a quasi-Newton minimization algorithm (LBFGS) with evenly distributed starting points using a Sobol sequence on a representative set $B_N = [-B, B]^N \cap \Omega_N$ of values of $\lambda$.

**Algorithm 4.3 ONLINE OPTIMIZATION**

**Input:** Reduced basis $m^1, \dots, m^N$, $B_N = [-B, B]^N \cap \Omega_N(m^1, \dots, m^N)$, starting points $\lambda_1, \dots, \lambda_L$ in $B_N$

for $l = 1, \dots, L$ do

Compute $\lambda_l^*$, $E_l^*$ minimizer and energy solution found by optimizing $E_{r,z}(\lambda)$ for $\lambda \in \Omega_N(m^1, \dots, m^N)$ with starting point $\lambda_l$ with a LBFGS algorithm

end for

**Output:** Minimizer $\lambda^* = \operatorname*{argmin}_{l=1,\dots,L} E_l^*$

Note that the LBFGS algorithm requires the explicit computation of the gradient of the energy, which can easily be computed from formula (4.7). Also, to ensure that a solution in the constraint set $\Omega_N(m^1, \dots, m^N)$ is found, points outside of $\Omega_N(m^1, \dots, m^N)$ are penalized by the function $\lambda \mapsto C + [\zeta^\lambda]^{-10}$, where $C$ is a large positive constant to ensure a return in the domain $\Omega_N \Omega_N(m^1, \dots, m^N)$ if a point outside of the constraint set is reached. Note that as the energy function explodes to $+\infty$ at the border of the domain, this is unlikely to happen.

**REMARK 4.4.** To avoid failures in line searches in the LBFGS algorithm caused by the low regularity of the function $E_{r,z}$, we actually consider a smoothed version of $E_{r,z}$. To do so, we replace the absolute value (and its derivative the sign function) responsible for the low regularity in a small interval $[-\varepsilon, \varepsilon]$ with $\varepsilon > 0$, by a cosine function $x \mapsto -\frac{2\varepsilon}{\pi} \cos(\frac{\pi}{2\varepsilon}x) + \varepsilon$. This new function is of class $C^2$ on $\mathbb{R}$.

**5. Numerical results.** In this section, we present the numerical results obtained with the offline and online algorithms presented above. The code for generating the figures in the following can be found at https://github.com/dussong/NonLinearReducedBasisOT. We focus on a system with two nuclei, i.e. $M = 2$, and with charges $z = (0.8, 1.1)$ and $r = (-r, r)$ for $r \in \mathbb{R}_+$. For the training set, the $r$'s are equally distributed on the interval $[0.5, 3]$ with $\#M_{tr} = 251$. On Figure 5.1, we plot the exact solutions called snapshots for three examples, namely $r = 0.5, 1.5, 2.5$.

**5.1. Offline phase.** We now present the results obtained by running the offline algorithm presented in Section 4.1. We present in Figure 5.2 the 7 first selected snapshots. We observe that the two first selected snapshots correspond to the extreme parameters 0.5 and 2.5, then the next ones are relatively well distributed across the parameter space. In Figure 5.3, we plot the decrease of the projection error in $\mathrm{MW}_2$-norm over the training set. We provide both the mean error on the training set and the maximum error. We observe
<figure><img src="image_3.png" /><figcaption>FIG. 5.2. First seven elements selected in the reduced basis in the offline phase.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIG. 5.3. Decay of the projection error in the offline phase.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIG. 5.4. Decay of the projection error in $W_2$-distance in the offline phase for $MW_2$-distance and $W_2$-distance greedy selection.</figcaption></figure>

that this projection error decreases very fast and seems exponentially decreasing. Moreover, we gain about two orders of magnitude between 2 and 15 added snapshots on the mean error.

In terms of computational cost, the most expensive part is the listing of the vertices of the constraint space (4.3), which increases exponentially with the number of selected snapshots. However, the code can be trivially parallelized, and is indeed running on multicores. Also, in the future, a global optimization solver
<figure><img src="image_6.png" /><figcaption>FIG. 5.5. Example of projections on bases with 2,3,5, and 8 elements for r = 1.266.</figcaption></figure>

such as Gurobi could be used instead of the listing of the vertices, possibly loosing the global optimality of
the found minimizer but gaining a lot in computational efficiency.

In Figure 5.4 we compare a greedy selection of snapshots using the MW₂-distance and the W₂-distance,
measuring all projection errors in W₂-distance. We observe that the selection in MW₂-distance clearly out-
performs the W₂-distance-based selection, showing that the MW₂-distance is more adapted to the considered
partial differential equation.

In Figure 5.5, we provide a few examples of projection on the reduced basis for a snapshot with parameter $r = 1.266$, which is in $M_z$ but not in the training set $M_{tr}$. We observe that the projections cannot be visually distinguished from the exact solution already when the reduced basis contains only 5 elements.

**5.2. Online phase.** In this section we provide results on the online optimization algorithm. First, recall that the energy minimization problem (4.6) is a global optimization problem so that Algorithm 4.3 may not necessarily return the global minimizer of the problem, possibly overestimating the presented error results compared to the exact ones. In practice, the parameters of Algorithm 4.3 were chosen as follows. We used a $L = 2000$ elements Sobol sequence covering the set $B_N = [-2, 2]^N \cap \Omega_N(m^1, \dots, m^N)$ as starting points $\lambda_1, \dots, \lambda_L$.

In Figure 5.6 we plot an example of energy landscape $\lambda \mapsto E_{\mathbf{r},\mathbf{z}}(\bar{\text{Bar}}_{\text{MW}_2}^{\lambda}(m^1, m^2))$ heatmap with $\mathbf{z} = (0.8, 1.1)$ and $\mathbf{r} = (-r, r)$ for $r = 2.15$, where $m^1$ and $m^2$ are the first two selected snapshots in the offline phase (see Figure 5.2). The white part corresponds to the outside of the domain $\Omega_2$. We already observe several local minima. Moreover, the energy is nonsmooth due to the absolute values appearing in the formula, see Remark 4.4 for more information on the smoothing technique.

We now provide the plot of the error in energy as a function of the number of selected snapshots in Figure 5.7. We provide both the maximum error and the mean error over a test set of 51 equally distributed elements for $r \in [0.5, 3]$. We observe that the energy maximum error decreases by three orders of magnitude from 2 to 8 snapshots, which is particularly encouraging. Adding more elements in the reduced basis does not seem to improve significantly the results. This may either be due to the increasing difficulty of solving the global optimization problem in larger dimension, but also to the smoothing of the energy functional that is used to avoid convergence problems.

Finally, we provide extrapolation examples. On Figures. 5.8 and 5.9, we show the projection of the
solution on the reduced basis for 2,3,5, and 8 snapshots. We observe that 8 snapshots seems sufficient to
obtain a satisfactory barycenter projection of the exact solution onto the reduced basis. More generally, we
<figure><img src="image_7.png" /><figcaption>FIG. 5.6. Heatmap of an energy functional (for $r = 2.15$) at barycenters between first two selected elements.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIG. 5.7. Decay of the energy error in the online phase for 51 equally distributed elements for $r \in [0.5, 3]$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 5.8. Extrapolation: example of energy projections for $r = 0.3$.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>FIG. 5.9. Extrapolation: example of energy projections for $r = 3.5$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIG. 5.10. Extrapolation: decay of the energy error over a set of 17 equally distributed elements for $r = 0, \dots, 0.48$.</figcaption></figure>

plot on Figure 5.10 the decay of the online error on 17 equally distributed elements ranging from $r = 0$ to $r = 0.48$. On Figure 5.11, we plot the energy online error on 21 equally distributed elements with $r = 3.05, \dots, 4$. We observe that we obtain very accurate results with only a few snapshots in the reduced basis, although the solutions are not in the parameter range of the training set, showing the nice extrapolation capabilities of the method.

**6. Proofs.** The aim of this section is to gather the proofs of the theoretical results stated in Section 3.

**6.1. Preliminary lemma.** Before going into the statement of the different theorems, we provide the following basic lemma, which bounds the error between a function with a lack of regularity at a few points and its best piecewise polynomial approximation. It will be used to obtain upper bound of Kolmogorov $n$-widths in Theorems 3.1 and 3.5.

**LEMMA 6.1.** Let $f$ be a real-valued function defined over a compact interval $I = [a, b]$, and let $a = x_0 < \dots < x_n = b$ be a mesh of maximal size $h$ on $I$. Suppose moreover that $f$ is of class $C^{p+1}$ on intervals $I_k = [x_k, x_{k+1}]$ except a few, named $I_{k_1}, \dots, I_{k_q}$ where it is of class $C^{p-1}$ with $f^{(p-1)}$ absolutely continuous.
<figure><img src="image_12.png" /><figcaption>FIG. 5.11. Extrapolation example: decay of the energy error over a set of 21 equally distributed elements for $r \in [3.05, \dots, 4]$.</figcaption></figure>

Then, defining

$$V_{n,p} = \operatorname{span}\{x \mapsto x^i 1_{I_k}(x)\}_{\substack{i=0,\dots,p \\ k=0,\dots,n-1}},$$

a vector space of dimension $n(p+1)$, we have the following projection error

$$\|f - P_{V_{n,p}} f\|_{L^2(I)} \le \left( \sqrt{b-a} \frac{\|f^{(p+1)}\|_{L^\infty(I \setminus (I_{k_1} \cup \dots \cup I_{k_q}))}}{(p+1)!} + \sqrt{q} \frac{\|f^{(p)}\|_{L^\infty(I_{k_1} \cup \dots \cup I_{k_q})}}{p!} \right) h^{p+\frac{1}{2}}.$$

*Proof.* We denote by $f_{n,p}$ the element of $V_{n,p}$ obtained as the $p$-th order Taylor polynomial at $x_k$ on each interval $I_k$ for $k \notin \{k_1, \dots, k_q\}$, i.e.

$$(6.1) \qquad \forall x \in I_k, \quad f_{n,p}(x) = \sum_{i=0}^{p} \frac{f^{(i)}(x_k)}{i!} (x - x_k)^i,$$

and as the $(p-1)$-th Taylor polynomial at $x_k$ on $I_k$ for the others intervals $I_k$, which is the same as in (6.1) but with a sum up to $p-1$. On intervals $I_k$, for $k \notin \{k_1, \dots, k_q\}$, the remainder can be written in the Lagrange form

$$f(x) - f_{n,p}(x) = \frac{f^{(p+1)}(\xi)}{(p+1)!} (x - x_k)^{p+1},$$

where $\xi \in I_k$. This yields

$$(6.2) \qquad \|f - f_{n,p}\|_{L^2(I_k)}^2 \le h \left( \frac{\|f^{(p+1)}\|_{L^\infty(I_k)}}{(p+1)!} \right)^2 h^{2(p+1)}.$$

On the remaining intervals $I_k$ for $k \in \{k_1, \dots, k_q\}$, as $f^{(p-1)}$ is absolutely continuous, $f^{(p)}$ is defined almost everywhere and we can use the integral form of the remainder and write

$$f(x) - f_{n,p}(x) = \int_{x_k}^{x} \frac{f^{(p)}(t)}{(p-1)!} (x-t)^{p-1} dt.$$

It then follows that for all $k \in \{k_1, \dots, k_q\}$

$$(6.3) \qquad \|f - f_{n,p}\|_{L^2(I_k)}^2 \le h \left( \frac{\|f^{(p)}\|_{L^\infty(I_k)}}{p!} \right)^2 h^{2p}.$$
Combining (6.2) and (6.3), we have

$$ \begin{aligned} \|f - f_{n,p}\|_{L^2(I)}^2 &\le \sum_{k \notin \{k_1, \dots, k_q\}} h \left( \frac{\|f^{(p+1)}\|_{L^\infty(I_k)}}{(p+1)!} \right)^2 h^{2(p+1)} + \sum_{k \in \{k_1, \dots, k_q\}} h \left( \frac{\|f^{(p)}\|_{L^\infty(I_k)}}{p!} \right)^2 h^{2p} \\ &\le \left( (b-a) \left( \frac{\|f^{(p+1)}\|_{L^\infty(I \setminus (I_{k_1} \cup \dots \cup I_{k_q}))}}{(p+1)!} \right)^2 + q \left( \frac{\|f^{(p)}\|_{L^\infty(I_{k_1} \cup \dots \cup I_{k_q})}}{p!} \right)^2 \right) h^{2p+1}, \end{aligned} $$

from which we easily obtain the result. □

## 6.2. Proof of Theorem 3.1.

*Proof of Theorem 3.1. Step 1:* We first prove the lower bound of (3.2). Since $L^2(-R, R)$ can be seen as a subset of $L^2(\mathbb{R})$ (by extending functions by 0 out of $(-R, R)$), it immediately holds that

$$ (6.4) \qquad d_{2,n}(\mathcal{M}_z^R, L^2(\mathbb{R})) \ge d_{2,n}(\mathcal{M}_z^R, L^2(-R, R)). $$

Let us then prove that there exists a constant $\tilde{c}_R > 0$ such that

$$ \tilde{c}_R n^{-\frac{3}{2}} \le d_{2,n}(\mathcal{M}_z^R, L^2(-R, R)). $$

We denote by $K$ the kernel

$$ K: (-R, R)^2 \ni (x, y) \mapsto \int_{-R}^{R} u_r(x) u_r(y) dr, $$

and introduce the integral operator defined by

$$ T_K: \quad \begin{array}{ccc} L^2(-R, R) & \to & L^2(-R, R) \\ \varphi & \mapsto & \int_{(-R, R)} K(\cdot, y) \varphi(y) dy. \end{array} $$

The operator $T_K$ is compact since $K \in L^2((-R, R)^2)$, self-adjoint because of the symmetry $K(x, y) = K(y, x)$, and non-negative. Indeed, let $f \in L^2(-R, R)$. We have

$$ \begin{aligned} \langle T_K f, f \rangle_{L^2(-R, R)} &= \int_{(-R, R)^2} \left( \int_{-R}^{R} u_r(x) u_r(y) dr \right) f(x) f(y) dx dy \\ &= \int_{-R}^{R} \langle u_r, f \rangle_{L^2(-R, R)}^2 dr \ge 0. \end{aligned} $$

Thus, from the spectral theorem, there exists a Hilbert basis $(\varphi_k)_{k \in \mathbb{N}^*}$ and a non-increasing sequence of non-negative real numbers $(\sigma_k)_{k \in \mathbb{N}^*}$ going to 0 as $k$ goes to $+\infty$ satisfying

$$ \forall k \in \mathbb{N}^*, T_K \varphi_k = \sigma_k \varphi_k. $$

Moreover, from [4, (1.46)], we can link the $L^2$ Kolmogorov $n$-width to the eigenvalues via the following formula

$$ (6.5) \qquad d_{2,n}(\mathcal{M}_z^R, L^2(-R, R)) = \sqrt{\sum_{k=n+1}^{+\infty} \sigma_k}. $$

Then, we will use the following intermediate lemma, which is proved below.
LEMMA 6.2. The spectrum of $T_K$ is equal to the set $\bigcup_{l \in \mathbb{N}^*} \{\lambda_l, \mu_l\}$, where for all $l \in \mathbb{N}^*$,

$$ \lambda_l = \frac{4z^4}{(2a_l^2 + z^3)^2} \quad \text{and} \quad \mu_l = \frac{4z^4}{(2b_l^2 + z^3)^2}, $$

with $a_l \in \left(\frac{(l-1)\pi}{R}, \frac{(l-1)\pi}{R} + \frac{\pi}{2R}\right)$ the $l$-th positive zero of the function

$x \mapsto x \sin(Rx) - z \cos(Rx)$ and $b_l \in \left(\frac{\pi}{2R} + \frac{(l-1)\pi}{R}, \frac{l\pi}{R}\right)$ the $l$-th positive zero of the function $x \mapsto x \cos(Rx) + z \sin(Rx)$. For all $l \in \mathbb{N}^*$, let us denote by

$$ \varphi_l : x \mapsto \cos(a_l x) \quad \text{and} \quad \psi_l : x \mapsto \sin(b_l x). $$

Then, $\varphi_l$ (respectively $\psi_l$) is an eigenvector of $T_K$ with corresponding eigenvalue $\lambda_l$ (respectively $\mu_l$). In addition, it holds that $\{\varphi_l, \psi_l\}_{l \in \mathbb{N}^*}$ is an orthogonal basis of $L^2(-R, R)$.

An immediate consequence of Lemma 6.2 is that for all $l \in \mathbb{N}^*$,

$$ \sigma_{2l-1} = \lambda_l \quad \text{and} \quad \sigma_{2l} = \mu_l. $$

In particular, for all $k \in \mathbb{N}^*$, $\sigma_k \ge \lambda_k$ and it thus holds that

$$ d_{2,n}(M_z^R, L^2(-R, R))^2 = \sum_{k=n+1}^{+\infty} \sigma_k \ge \sum_{k=n+1}^{+\infty} \lambda_k. $$

In particular, it can be easily seen that the sequence $(\lambda_k)_{k \in \mathbb{N}^*}$ is decreasing and that there exists a constant $\tilde{c}_R > 0$ such that $\lambda_k \ge \tilde{c}_R k^{-4}$. Therefore, combining (6.4) and (6.5), we obtain that

$$ d_{\infty,n}(M_z^R, L^2(\mathbb{R})) = d_{2,n}(M_z^R, L^2(\mathbb{R})) \ge \tilde{c}_R n^{-\frac{3}{2}}, $$

which proves the lower bound of (3.2).

**Step 2:** We now prove the upper bound of (3.1). For $n \in \mathbb{N}^*$, let us take $r \in [-R, R]$ and $x_0 = -R < \dots < x_n = R$ the equidistant subdivision of $I = [-R, R]$. We also define the $2n$-dimensional vector space

$$ V_n = \operatorname{span} \{x \mapsto x^i 1_{I_k}(x)\}_{\substack{i=0,1 \\ k=0,\dots,n-1}}, $$

where $I_k = [x_k, x_{k+1}]$, and denote by $k_r$ the index such that $r \in I_{k_r}$. We also define the vector space

$$ V = \operatorname{span} \{x \mapsto e^{zx} 1_{(-\infty, -R)}, x \mapsto e^{-zx} 1_{(R, +\infty)}\}. $$

Using these vector spaces, it is clear from the shape of the solution $u_r$ (see (2.5)) that

$$ \|u_r - P_{V \oplus V_n} u_r\|_{L^2(\mathbb{R})} = \|u_r - P_{V_n} u_r\|_{L^2(I)}. $$

Moreover, applying Lemma 6.1 with $p = 1$ and $f = u_r$ which is twice differentiable over the intervals $I_k$ for $k \ne k_r$ and absolutely continuous on $I_{k_r}$, we have

$$ \|u_r - P_{V_n} u_r\|_{L^2(I)} \le (2R)^{\frac{3}{2}} \left( \frac{\sqrt{2R}}{2} \|u_r''\|_{L^\infty(I \setminus I_{k_r})} + \|u_r'\|_{L^\infty(I_{k_r})} \right) n^{-\frac{3}{2}}. $$

And since for $x \in \mathbb{R}$,

$$ |u_r'(x)| = \frac{z^2}{2} e^{-z|x-r|} \le \frac{z^2}{2} \quad \text{and} \quad |u_r''(x)| = \frac{z^3}{2} e^{-z|x-r|} \le \frac{z^3}{2}, $$
there exists a positive constant $C_R$ depending on $R$ such that

$$ \|u_r - P_{V_n} u_r\|_{L^2(I)} \le C_R n^{-\frac{3}{2}}. $$

We can conclude by writing that

$$ d_{\infty, 2n+2}(\mathcal{M}_z^R, L^2(\mathbb{R})) \le \sup_{r \in [-R, R]} \|u_r - P_{V \oplus V_n} u_r\|_{L^2(\mathbb{R})} \le C_R n^{-\frac{3}{2}}. $$

**Step 3:** The two remaining bounds are easily deduced remarking that for any $n$-dimensional subspace $V_n \subset L^2(\mathbb{R})$, we have

$$ \left( \int_{-R}^{R} \|u_r - P_{V_n} u_r\|_{L^2(\mathbb{R})}^2 dr \right)^{\frac{1}{2}} \le \sqrt{2R} \sup_{r \in [-R, R]} \|u_r - P_{V_n} u_r\|_{L^2(\mathbb{R})}. $$

which concludes the proof.

We provide now the proof of Lemma 6.2.

*Proof of Lemma 6.2.* Using the definition of the kernel $K$ and rearranging the integrals we remark that, for all $k, l \in \mathbb{N}^*$,

$$ \langle \varphi_k, T_K \varphi_l \rangle = \int_{-R}^{R} \left( \int_{-R}^{R} u_r(x) \cos(a_k x) dx \right) \left( \int_{-R}^{R} u_r(y) \cos(a_l y) dy \right) dr \\ = \langle T_S \cos(a_k), T_S \cos(a_l) \rangle, $$

where $T_S$ is the compact and self-adjoint operator defined by

$$ (6.6) \qquad \begin{array}{ccc} L^2(-R, R) & \to & L^2(-R, R) \\ T_S: & \varphi & \mapsto \left( r \mapsto \int_{-R}^{R} u_r(x) \varphi(x) dx \right). \end{array} $$

The self-adjointness of $T_S$ stems from the fact that $u_r(x) = u_x(r)$ for all $(x, r) \in (-R, R)$. Similarly, it holds that for all $k, l \in \mathbb{N}^*$,

$$ \langle \varphi_k, T_K \psi_l \rangle = \langle T_S \varphi_k, T_S \psi_l \rangle \quad \text{and} \quad \langle \psi_k, T_K \psi_l \rangle = \langle T_S \psi_k, T_S \psi_l \rangle. $$

For $k \in \mathbb{N}^*$, we now compute $T_S \cos(a_k \cdot)$:

$$ (6.7) \qquad T_S \cos(a_k \cdot)(r) = \frac{z}{2} \left( \int_{-R}^{r} e^{z(x-r)} \cos(a_k x) dx + \int_{r}^{R} e^{-z(x-r)} \cos(a_k x) dx \right) $$

We compute the two integrals using two integrations by parts. First

$$ \begin{aligned} I^{-}(r) &:= \int_{-R}^{r} e^{z(x-r)} \cos(a_k x) dx \\ &= \frac{1}{a_k} \left[ \sin(a_k x) e^{z(x-r)} \right]_{-R}^{r} - \frac{z}{a_k} \int_{-R}^{r} e^{z(x-r)} \sin(a_k x) dx \\ &= \frac{\sin(a_k r)}{a_k} + \frac{\sin(a_k R)}{a_k} e^{-z(R+r)} + \frac{z}{a_k^2} \left[ \cos(a_k x) e^{z(x-r)} \right]_{-R}^{r} - \frac{z^2}{a_k^2} I^{-}(r) \\ &= \frac{\sin(a_k r)}{a_k} + \frac{z}{a_k^2} \cos(a_k r) + \frac{a_k \sin(a_k R) - z \cos(a_k R)}{a_k^2} e^{-z(R+r)} - \frac{z^2}{a_k^2} I^{-}(r) \\ &= \frac{\sin(a_k r)}{a_k} + \frac{z}{a_k^2} \cos(a_k r) - \frac{z^2}{a_k^2} I^{-}(r), \end{aligned} $$
noting that $a_k \sin(a_k R) - z \cos(a_k R) = 0$. In the same manner we can consider the other integral $I^+(r) := \int_r^R e^{-z(x-r)} \cos(a_k x)dx = I^-((-r))$ to find

$$I^+(r) = -\frac{\sin(a_k r)}{a_k} + \frac{z}{a_k^2} \cos(a_k r) - \frac{z^2}{a_k^2} I^+(r).$$

Hence, continuing from (6.7),

$$T_S \cos(a_k \cdot)(r) = \frac{z}{2} \left( \frac{2z}{a_k^2} \cos(a_k r) - \frac{z^2}{a_k^2} T_S \cos(a_k \cdot)(r) \right),$$

which means that

$$T_S \cos(a_k \cdot)(r) = \frac{2z^2}{2a_k^2 + z^3} \cos(a_k r).$$

Similarly, for $k \in \mathbb{N}^*$, we now compute $T_S \sin(b_k \cdot)$:

$$(6.8) \qquad T_S \sin(b_k \cdot)(r) = \frac{z}{2} \left( \int_{-R}^{r} e^{z(x-r)} \sin(b_k x)dx + \int_{r}^{R} e^{-z(x-r)} \sin(b_k x)dx \right)$$

We compute the two integrals using two integrations by parts. First

$$
\begin{align*}
J^{-}(r) &:= \int_{-R}^{r} e^{z(x-r)} \sin(b_k x) dx \\
&= \frac{1}{b_k} \left[ -\cos(b_k x) e^{z(x-r)} \right]_{-R}^{r} + \frac{z}{b_k} \int_{-R}^{r} e^{z(x-r)} \cos(b_k x) dx \\
&= -\frac{\cos(b_k r)}{b_k} + \frac{\cos(b_k R)}{b_k} e^{-z(R+r)} + \frac{z}{b_k^2} \left[ \sin(b_k x) e^{z(x-r)} \right]_{-R}^{r} - \frac{z^2}{b_k^2} J^{-}(r) \\
&= -\frac{\cos(b_k r)}{b_k} + \frac{z}{b_k^2} \sin(b_k r) + \frac{b_k \cos(b_k R) + z \sin(b_k R)}{b_k^2} e^{-z(R+r)} - \frac{z^2}{b_k^2} J^{-}(r) \\
&= -\frac{\cos(b_k r)}{b_k} + \frac{z}{b_k^2} \sin(b_k r) - \frac{z^2}{b_k^2} J^{-}(r),
\end{align*}
$$

noting that $b_k \cos(b_k R) + z \sin(b_k R) = 0$. In the same manner we can consider the other integral $J^+(r) := \int_r^R e^{-z(x-r)} \sin(b_k x)dx = -J^-((-r))$ to find

$$J^+(r) = \frac{\cos(b_k r)}{b_k} + \frac{z}{b_k^2} \sin(b_k r) - \frac{z^2}{b_k^2} J^+(r).$$

Hence, continuing from (6.8)

$$T_S \sin(b_k \cdot)(r) = \frac{z}{2} \left( \frac{2z}{b_k^2} \sin(b_k r) - \frac{z^2}{b_k^2} T_S \sin(b_k \cdot)(r) \right),$$

which means that

$$T_S \sin(b_k \cdot)(r) = \frac{2z^2}{2b_k^2 + z^3} \sin(b_k r).$$

Hence, for all $k, l \in \mathbb{N}^*$, the functions $\varphi_k$ and $\psi_l$ are eigenvectors of $T_S$ with respective distinct eigenvalues $\sigma_k = \frac{2z^2}{2a_k^2 + z^3}$ and $\tau_l = \frac{2z^2}{2b_k^2 + z^3}$. Thus, denoting by $\lambda_k = \sigma_k^2$ and by $\mu_k = \tau_k^2$ for all $k \in \mathbb{N}^*$, we obtain that
$\bigcup_{k \in \mathbb{N}^*} \{\lambda_k, \mu_k\} \subset \sigma(T_K)$ where $\sigma(T_K)$ is the spectrum of $T_K$. It can also be easily checked that $\{\varphi_k, \psi_k\}_{k \in \mathbb{N}^*}$ forms an orthogonal family of functions of $L^2(-R, R)$.

It remains to prove that $\text{Span}\{\phi_k, \psi_k\}_{k \in \mathbb{N}^*}$ is dense in $L^2(-R, R)$. From [13, Theorem 2], it is clear that $(\psi_k)_{k \in \mathbb{N}^*}$ is an orthogonal basis of the set of odd functions of $L^2(-R, R)$. From [26], it holds similarly that $(\phi_k)_{k \in \mathbb{N}^*}$ is an orthogonal basis of the set of even functions of $L^2(-R, R)$, hence the desired result. In particular, we then have that $\bigcup_{k \in \mathbb{N}^*} \{\lambda_k, \mu_k\} = \sigma(T_K)$. $\square$

### 6.3. Proof of Theorem 3.5.

*Proof of Theorem 3.5.* First note that from the definition of the solution (2.2) alongside its special form in this case (2.6), we can explicitly compute the cumulative distribution function as well as inverse cumulative distribution function as

$$ \text{cdf}_{u_r}(x) = \begin{cases} \frac{1}{2} \cosh(\zeta_r r) e^{\zeta_r x}, & x < -r, \\ \frac{1}{2} (1 + e^{-\zeta_r r} \sinh(\zeta_r x)), & -r \le x < r, \\ 1 - \frac{1}{2} \cosh(\zeta_r r) e^{-\zeta_r x}, & x \le -r, \end{cases} $$

and

$$ (6.9) \qquad \text{icdf}_{u_r}(s) = \begin{cases} \frac{1}{\zeta_r}(\ln(2s) - \ln(\cosh(\zeta_r r))), & 0 < s < s_r, \\ \frac{1}{\zeta_r}\arcsin h(e^{\zeta_r r}(2s-1)), & s_r \le s < 1-s_r, \\ -\frac{1}{\zeta_r}(\ln(2(1-s)) - \ln(\cosh(\zeta_r r))), & 1-s_r < s \le 1, \end{cases} $$

with

$$ s_r = \frac{1}{2} \cosh(\zeta_r r) e^{-\zeta_r r} = \frac{1}{4} (1 + e^{-2\zeta_r r}) = \frac{1}{4z} \zeta_r, $$

thanks to (2.7).

Now, for $n \in \mathbb{N}^*$, we introduce $\frac{1}{4} = s_0 < \dots < s_n = \frac{3}{4}$ the equidistant subdivision of $I := [\frac{1}{4}, \frac{3}{4}]$, and $V$, $V_n$ the vector spaces respectively defined by

$$ (6.10) \qquad V = \operatorname{span} \left\{ s \mapsto \ln(2s) \mathbf{1}_{(0, \frac{1}{4})}(s), s \mapsto \mathbf{1}_{(0, \frac{1}{4})}(s), s \mapsto \ln(2(1-s)) \mathbf{1}_{(\frac{3}{4}, 1)}(s), s \mapsto \mathbf{1}_{(\frac{3}{4}, 1)}(s) \right\}, $$

and

$$ (6.11) \qquad V_n = \operatorname{span} \left\{ s \mapsto s^i \mathbf{1}_{I_k}(s) \right\}_{\substack{i=0,1,2 \\ k=0,\dots,n-1}}, $$

where $I_k$ is the interval $[s_k, s_{k+1}]$. Since

$$ (6.12) \qquad d_{\infty, 3n+4}(\mathcal{T}_R^-, L^2(0,1)) \le \sup_{r \in [0, R]} \|\text{icdf}_{u_r} - P_{V \oplus V_n} \text{icdf}_{u_r}\|_{L^2(0,1)}, $$

we are interested in estimating the error $\|\text{icdf}_{u_r} - P_{V \oplus V_n} \text{icdf}_{u_r}\|_{L^2(0,1)}$ for all $r \in [0, R]$.

First, for $s \in (0, \frac{1}{4}) \cup (\frac{3}{4}, 1)$, it is clear that for all $r \in [0, R]$, the error $|\text{icdf}_{u_r}(s) - P_{V \oplus V_n} \text{icdf}_{u_r}(s)|$ is equal to 0 from the definition of $V$.

On the remaining interval $I$, we use Lemma 6.1 with $p = 2$ and $f = \text{icdf}_{u_r}$ which is three times differentiable on the intervals $I_k$ for $k \notin \{k_r, \tilde{k}_r\}$ where $k_r, \tilde{k}_r$ are the indices such that $s_r \in I_{k_r}$, and $1 - s_r \in I_{\tilde{k}_r}$, and $f'$ is absolutely continuous on the intervals $I_{k_r}$ and $I_{\tilde{k}_r}$. We have

$$ (6.13) \qquad \|\text{icdf}_{u_r} - P_{V \oplus V_n} \text{icdf}_{u_r}\|_{L^2(I)} \le \left( \frac{1}{6\sqrt{2}} \|\text{icdf}_{u_r}^{(3)}\|_{L^\infty(I \setminus (I_{k_r} \cup I_{\tilde{k}_r}))} \frac{\sqrt{2}}{2} \|\text{icdf}_{u_r}^{(2)}\|_{L^\infty(I_{k_r} \cup I_{\tilde{k}_r})} \right) (2n)^{-\frac{5}{2}}. $$
From (6.9), we have

$$ \text{icdf}_{u_r}^{(2)}(s) = \begin{cases} -\frac{1}{\zeta_r s^2}, & \frac{1}{4} \le s < s_r, \\ -\frac{4}{\zeta_r} \frac{2s-1}{(e^{-2\zeta_r r} + (2s-1)^2)^{\frac{3}{2}}}, & s_r \le s \le \frac{1}{2}, \end{cases} $$

and

$$ \text{icdf}_{u_r}^{(3)}(s) = \begin{cases} \frac{2}{\zeta_r s^3}, & \frac{1}{4} \le s < s_r, \\ \frac{8}{\zeta_r} \frac{2(2s-1)^2 - e^{-2\zeta_r r}}{(e^{-2\zeta_r r} + (2s-1)^2)^{\frac{5}{2}}}, & s_r \le s \le \frac{1}{2}. \end{cases} $$

Given these formulas, we can easily bound $\text{icdf}_{u_r}^{(2)}$ and $\text{icdf}_{u_r}^{(3)}$ on $I$. On $[\frac{1}{4}, s_r]$, we have

$$ |\text{icdf}_{u_r}^{(2)}(s)| = \frac{1}{\zeta_r s^2} \le \frac{16}{\zeta_r} \le \frac{16}{z}, \quad \text{and} \quad |\text{icdf}_{u_r}^{(3)}(s)| = \frac{2}{\zeta_r s^3} \le \frac{128}{\zeta_r} \le \frac{128}{z}, $$

since $\zeta_r \ge z$ by (2.7). On $[s_r, \frac{1}{2}]$, we have

$$ |\text{icdf}_{u_r}^{(2)}(s)| = \frac{4}{\zeta_r} \frac{1-2s}{((2s-1)^2 + e^{-2\zeta_r r})^{\frac{3}{2}}} \le \frac{4}{\zeta_r} e^{3\zeta_r r} (1-2s) \le \frac{2}{\zeta_r} e^{3\zeta_r r} \le \frac{2}{z} e^{3\zeta_R R}, $$

and

$$ |\text{icdf}_{u_r}^{(3)}(s)| = \frac{8}{\zeta_r} \frac{|2(2s-1)^2 - e^{-2\zeta_r r}|}{((2s-1)^2 + e^{-2\zeta_r r})^{\frac{5}{2}}} \le \frac{8}{\zeta_r} e^{5\zeta_r r} \le \frac{8}{z} e^{5\zeta_R R} $$

since

$$ |2(2s-1)^2 - e^{-2\zeta_r r}| \le \max(e^{-2\zeta_r r}, 2(2s_r - 1)^2 - e^{-2\zeta_r r}) \le 1. $$

The estimate for the remaining part $[\frac{1}{2}, \frac{3}{4}]$ follows noting that $\text{icdf}(s-1/2) = -\text{icdf}(s+1/2)$. It then follows from (6.13) and the above bounds for the second and third derivatives of $\text{icdf}_{u_r}$ that

$$ (6.14) \qquad \|\text{icdf}_{u_r} - P_{V \oplus V_n} \text{icdf}_{u_r}\|_{L^2(0,1)} \le \frac{1}{4z} \left( \frac{2}{3} e^{5\zeta_R R} + e^{3\zeta_R R} \right) n^{-\frac{5}{2}} \le C e^{5\zeta_R R} n^{-\frac{5}{2}}, $$

for $e^{3\zeta_R R} \ge 8$ and $e^{5\zeta_R R} \ge 16$, with $C > 0$ a real constant independent of $R$. We conclude the proof by combining (6.12) and (6.14). $\square$

## 6.4. Proof of Theorem 3.6.

*Proof of Theorem 3.6.* Let $R > 0$. For any $r \ge R$, we define the functions $g_r$ and $h_r$ as

$$ g_r = \frac{1}{2} \text{cdf}_{S_{\zeta_r, -r}}, $$

$$ h_r = \frac{1}{2} \text{cdf}_{S_{\zeta_r, r}}, $$

the cumulative distribution functions of the two Slater distributions of $u_r$, which clearly yields

$$ \text{cdf}_{u_r} = g_r + h_r. $$

Moreover, we have

$$ g_r(x) = \begin{cases} \frac{1}{4} e^{\zeta_r (x+r)}, & x < -r, \\ \frac{1}{2} \left(1 - \frac{1}{2} e^{-\zeta_r (x+r)}\right), & x \ge -r, \end{cases} \qquad h_r(x) = \begin{cases} \frac{1}{4} e^{\zeta_r (x-r)}, & x < -r, \\ \frac{1}{2} \left(1 - \frac{1}{2} e^{-\zeta_r (x-r)}\right), & x \ge -r, \end{cases} $$

and

$$ g_r^{-1}(s) = \begin{cases} \frac{1}{\zeta_r} \ln(4s) - r, & s \in (0, \frac{1}{4}), \\ -\frac{1}{\zeta_r} \ln(2(1-2s)) - r, & s \in [\frac{1}{4}, \frac{1}{2}), \end{cases} \qquad h_r^{-1}(s) = \begin{cases} \frac{1}{\zeta_r} \ln(4s) + r, & s \in (0, \frac{1}{4}), \\ -\frac{1}{\zeta_r} \ln(2(1-2s)) + r, & s \in [\frac{1}{4}, \frac{1}{2}). \end{cases} $$
Consider now the vector space $V_R$

$$ (6.15) \qquad V_R = \operatorname{span} \left\{ \mathbb{1}_{(0, \frac{1}{2})}, g_R^{-1}, \mathbb{1}_{(\frac{1}{2}, 1)}, \left(\frac{1}{2} + h_R\right)^{-1} \right\}. $$

Our aim is to prove that there exists a positive constant $C$ independent of $r$ and $R$ such that for all $r \ge R$, the following bound for the projection error holds:

$$ \| \mathrm{icdf}_{u_r} - P_{V_R} \mathrm{icdf}_{u_r} \|_{L^2(0,1)} \le C R e^{-\frac{1}{2} \zeta_R R}, $$

which will yield the desired result. It roughly states that solutions $u_r$ with a large position parameter $r$ are close to the mean of the two solutions with $M = 1$ centered in $-r$ and $r$ with a fixed charge $z$.

By parity with respect to $s = 1/2$ of $\mathrm{icdf}_{u_r}$ and remarking that $V_R$ contains functions on $(0, \frac{1}{2})$ and their exact translations on $(\frac{1}{2}, 1)$, we only need to focus on the left interval $(0, \frac{1}{2})$. Now, as $g_r^{-1}$ is in $V_R$, we bound the projection by

$$ (6.16) \qquad \begin{aligned} \| \mathrm{icdf}_{u_r} - P_{V_R} \mathrm{icdf}_{u_r} \|_{L^2(0, \frac{1}{2})}^2 & \le \| \mathrm{icdf}_{u_r} - g_r^{-1} \|_{L^2(0, \frac{1}{2})}^2 \\ & = \int_0^{\frac{1}{2}} |\mathrm{icdf}_{u_r}(s) - g_r^{-1}(s)|^2 ds. \end{aligned} $$

Next, we consider separately the integrals on $(0, g_r(0)]$ and $(g_r(0), \frac{1}{2}]$. To estimate the first part, we remark that $g_r^{-1} - \mathrm{icdf}_{u_r}$ is positive and increasing on the interval $(0, \frac{1}{2}]$. Hence

$$ \int_0^{g_r(0)} |\mathrm{icdf}_{u_r}(s) - g_r^{-1}(s)|^2 ds \le |\mathrm{icdf}_{u_r}(g_r(0))| \int_0^{g_r(0)} (g_r^{-1}(s) - \mathrm{icdf}_{u_r}(s)) ds. $$

Using inverses, we have that

$$ \int_0^{g_r(0)} (g_r^{-1}(s) - \mathrm{icdf}_{u_r}(s)) ds \le \int_{-\infty}^0 (\mathrm{cdf}_{u_r}(x) - g_r(x)) dx = \int_{-\infty}^0 h_r(x) dx = \frac{1}{4\zeta_r} e^{-\zeta_r r}, $$

and

$$ |\mathrm{icdf}_{u_r}(g_r(0))| = \frac{1}{\zeta_r} \arcsin h \left( \frac{1}{2} \right). $$

From this, it follows that

$$ \int_0^{g_r(0)} |\mathrm{icdf}_{u_r}(s) - g_r^{-1}(s)|^2 ds \le \frac{1}{4\zeta_r^2} \arcsin h \left( \frac{1}{2} \right) e^{-\zeta_r r}. $$

We now bound the second part:

$$ \begin{aligned} \int_{g_r(0)}^{\frac{1}{2}} |\mathrm{icdf}_{u_r}(s) - g_r^{-1}(s)|^2 ds & \le 2 \int_{g_r(0)}^{\frac{1}{2}} |\mathrm{icdf}_{u_r}(s) + r|^2 ds + \frac{2}{\zeta_r^2} \int_{g_r(0)}^{\frac{1}{2}} |\ln(2(1-2s))|^2 ds \\ & \le 2r^2 \left(\frac{1}{2} - g_r(0)\right) + \frac{1}{2\zeta_r^2} \int_0^{e^{-\zeta_r r}} |\ln(t)|^2 dt \\ & \le \frac{r^2}{2} e^{-\zeta_r r} + \frac{1}{2\zeta_r^2} (\zeta_r^2 r^2 + 2\zeta_r r + 2) e^{-\zeta_r r} \\ & = \left(r^2 + \frac{r}{\zeta_r} + \frac{2}{\zeta_r^2}\right) e^{-\zeta_r r}. \end{aligned} $$

Combining these two estimates, we bound the integral in (6.16) by

$$ \int_0^{\frac{1}{2}} |\mathrm{icdf}_{u_r}(s) - g_r^{-1}(s)|^2 ds \le Cr^2 e^{-\zeta_r r}, $$
with $C$ a positive real constant independent $r$ and $R$. From this, we easily deduce that

$$ \|\mathrm{icdf}_{u_r} - P_{V_R}\mathrm{icdf}_{u_r}\|_{L^2(0,1)} \le CRe^{-\frac{1}{2}\zeta_R R}, $$

for $R$ large enough such that for all $r \ge R$, we have $r^2e^{-\zeta_r r} \le R^2e^{-\zeta_R R}$. □

### 6.5. Proof of Theorem 3.7.

*Proof of Theorem 3.7.* We consider the vector spaces $V$ and $V_n$ as respectively defined in (6.10) and (6.11) and the vector space $V_R$ as defined in (6.15).

First of all, fixing $\varepsilon > 0$ there exists a constant $C_\varepsilon > 0$ such that

$$ \forall r \ge R, \quad \|\mathrm{icdf}_{u_r} - P_{V_R}\mathrm{icdf}_{u_r}\|_{L^2(0,1)} \le CRe^{-\frac{\varepsilon}{2}R} \le C_\varepsilon e^{-(\frac{\varepsilon}{2}-\varepsilon)R}, $$

from Theorem 3.6 and (2.7). Also, from (6.14) and (2.7), we have

$$ \forall r \in [0, R], \quad \|\mathrm{icdf}_{u_r} - P_{V \oplus V_n}\mathrm{icdf}_{u_r}\|_{L^2(0,1)} \le Ce^{10zR}n^{-\frac{5}{2}}. $$

Hence,

$$ \begin{aligned} d_{2,3n+4}(\mathcal{T}, L^2(0,1)) &\le \sup_{r \in \mathbb{R}^+} \|\mathrm{icdf}_{u_r} - P_{(V \oplus V_n)+V_R}\mathrm{icdf}_{u_r}\|_{L^2(0,1)} \\ &\le \max \left( \sup_{r \in [0,R]} \|\mathrm{icdf}_{u_r} - P_{V \oplus V_n}\mathrm{icdf}_{u_r}\|_{L^2(0,1)}; \sup_{r \ge R} \|\mathrm{icdf}_{u_r} - P_{V_R}\mathrm{icdf}_{u_r}\|_{L^2(0,1)} \right) \\ &\le \max \left( C e^{10zR} n^{-\frac{5}{2}}; C_\varepsilon e^{-(\frac{\varepsilon}{2}-\varepsilon)R} \right). \end{aligned} $$

Choosing $R_n$ satisfying $Ce^{10zR_n}n^{-\frac{5}{2}} = C_\varepsilon e^{-(\frac{\varepsilon}{2}-\varepsilon)R_n}$, that is

$$ R_n = \frac{\ln \frac{C_\varepsilon}{C} + \frac{5}{2} \ln n}{\frac{21}{2} z - \varepsilon}, $$

we obtain the result. □

### 6.6. Proof of Theorem 3.8.

*Proof of Theorem 3.8.* Let $m^1$ and $m^2$ be two symmetric mixtures of $K = 2$ elements, with parameters $r^1 = 1, r^2 = 2, \zeta^1 = 1$ and $\zeta^2 = 2$, that is

$$ m^1 = \frac{1}{2} (S_{\zeta^1, -r^1} + S_{\zeta^1, r^1}) \quad \text{and} \quad m^2 = \frac{1}{2} (S_{\zeta^2, -r^2} + S_{\zeta^2, r^2}). $$

Let us denote by $m_1^1 := S_{\zeta^1, -r^1}$, $m_2^1 := S_{\zeta^1, r^1}$, $m_1^2 := S_{\zeta^2, -r^2}$ and $m_2^2 := S_{\zeta^2, r^2}$ so that $m^1 = \frac{1}{2}(m_1^1 + m_2^1)$ and $m^2 = \frac{1}{2}(m_1^2 + m_2^2)$. Since

$$ W_2(m_1^1, m_1^2)^2 = W_2(m_2^1, m_2^2)^2 = (r^1 - r^2)^2 + 2 \left( \frac{1}{\zeta^1} - \frac{1}{\zeta^2} \right)^2 = \frac{3}{2}, $$

and

$$ W_2(m_1^1, m_2^1)^2 = W_2(m_2^1, m_2^1)^2 = (r^1 + r^2)^2 + 2 \left( \frac{1}{\zeta^1} - \frac{1}{\zeta^2} \right)^2 = \frac{19}{2}, $$

it is easy to determine the only solution $w^*(m; \bar{\lambda})$ of (2.20) with $\bar{\lambda} = (1/2, 1/2)$ and $m = (m^1, m^2)$. Indeed, in this case, the distance $\mathrm{MW}_2(m^1, m^2)$ simplifies to

$$ \begin{array}{l} \mathrm{MW}_2(m^1, m^2) = 2 \qquad \min_{\substack{w_{11}, w_{12} \in \mathbb{R}_+ \\ w_{11} + w_{12} = \frac{1}{2}}} w_{11} W_2(m_1^1, m_1^2)^2 + w_{12} W_2(m_2^1, m_2^2)^2 = \min_{\substack{w_{11}, w_{12} \in \mathbb{R}_+ \\ w_{11} + w_{12} = \frac{1}{2}}} 3w_{11} + 19w_{12}, \end{array} $$
which is clearly attained for $w_{11} = \frac{1}{2}$ and $w_{12} = 0$. Hence, the unique solution $w^*(m; \bar{\lambda}) = \begin{pmatrix} \frac{1}{2} & 0 \\ 0 & \frac{1}{2} \end{pmatrix}$.

Then the barycenters between $m^1$ and $m^2$ for a weight

$$ \lambda = (\lambda_1, \lambda_2) \in \Omega_2(m^1, m^2) := \left\{ (\lambda_1, \lambda_2) \in \mathbb{R}^2 : \lambda_1 + \frac{1}{2}\lambda_2 > 0 \right\} $$

(since for all $(\lambda_1, \lambda_2) \in \mathbb{R}^2$, $\frac{\lambda_1}{\zeta_1} + \frac{\lambda_2}{\zeta_2} = \lambda_1 + \frac{1}{2}\lambda_2$) are equal to

$$ \overline{\mathrm{Bar}}_{\mathrm{MW}_2}^{\lambda}(m^1, m^2) = \frac{1}{2} (S_{-r\lambda, \zeta\lambda} + S_{r\lambda, \zeta\lambda}), $$

where

$$ r^{\lambda} = \lambda_{1}r^{1} + \lambda_{2}r^{2} = \lambda_{1} + 2\lambda_{2} \quad \text{and} \quad \zeta^{\lambda} = \left[ \frac{\lambda_{1}}{\zeta^{1}} + \frac{\lambda_{2}}{\zeta^{2}} \right]^{-1} = \left[ \lambda_{1} + \frac{1}{2}\lambda_{2} \right]^{-1}. $$

Let $u_r \in \mathcal{M}_z$. To show that $u_r$ is indeed an approximate mixture Wasserstein barycenter of $m^1$ and $m^2$, we find $\lambda \in \mathbb{R}^2$ such that $r^\lambda = r$ and $\zeta^\lambda = \zeta_r$, which amounts to solving the linear system

$$ \begin{cases} \lambda_1 + 2\lambda_2 = r \\ \lambda_1 + \frac{1}{2}\lambda_2 = \zeta_r^{-1}, \end{cases} $$

the unique solution $(\lambda_1, \lambda_2)$ of which is given by $\lambda_1 = \frac{1}{3}(-r + 4\zeta_r^{-1})$ and $\lambda_2 = \frac{2}{3}(r - \zeta_r^{-1})$) and can be easily checked to belong to $\Omega_2(m^1, m^2)$. $\square$

**7. Conclusion.** In this work, we have focused our attention on a one-dimensional parametrized toy problem, which is insightful in terms of the difficulties faced by standard model order reduction methods to accelerate the resolution of parametrized electronic structure problems. We proved that the linear Kolmogorov n-width of solution sets for this equation decays at a slow algebraic rate with respect to n. We proved that modified Kolmogorov widths, based on optimal transport tools, in particular Wasserstein mixture distances, decay much faster. Motivated by this result, we proposed a modified greedy algorithm, precisely based on mixture Wasserstein distances and corresponding barycenters, which gives highly encouraging results. Our aim is now to export the ideas and concepts of the present work in order to build efficient reduced order models for realistic parametrized electronic structure problems in a forthcoming work.

The nonlinear reduced basis strategy presented in this work can, in principle, be adapted to application domains beyond electronic structure calculations. In particular, one could consider atomic measures different from Slater distributions with localized support—for example, Wigner distributions—in order to construct nonlinear reduced order models for problems defined on bounded subsets of $\mathbb{R}^d$, as is often the case in mechanical engineering applications. However, the present approach has some important limitations: (i) Outputs of interest must be functions with constant sign (at least up to some simple transformations) to relate them to probability distributions. (ii) One cannot prescribe complex boundary conditions on the solutions of the reduced order model, which may be a bottleneck towards the extension of the present methodology to realistic engineering applications.

**Acknowledgments.** The authors thank Alexandre Nou for pointing out the Paley–Wiener results linked to the proof presented in Appendix. They also thank Christoph Ortner for interesting discussions.

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreements EMC2 No 810367 and HighLEAP No 101077204). This work was supported by the French ‘Investissements d’Avenir’ program, project Agence Nationale de la Recherche (ISITE-BFC) (contract ANR-15-IDEX-0003), as well as the ANR project NUMERIQ (ANR-24-CE46-2255). GD was also supported by the Ecole des Ponts-ParisTech. GD acknowledges the support of the region Bourgogne Franche-Comté. VE acknowledges support from the ANR project COMODO (ANR-19-CE46-0002).
REFERENCES

[1] M. AGUEH AND G. CARLIER, *Barycenters in the wasserstein space*, SIAM J. Math. Anal., 43 (2011), pp. 904-924.

[2] P. C. ÁLVAREZ-ESTEBAN, E. DEL BARRIO, J. A. CUESTA-ALBERTOS, AND C. MATRÁN, *A fixed-point approach to barycenters in wasserstein space*, J. Math. Anal. Appl., 441 (2016), pp. 744-762.

[3] Y. CHEN, T. T. GEORGIOU, AND A. TANNENBAUM, *Optimal transport for gaussian mixture models*, IEEE Access, 7 (2018), pp. 6269-6278.

[4] A. COHEN AND R. DeVORE, *Approximation of high-dimensional parametric PDEs*, Acta Numer., 24 (2015), pp. 1-159.

[5] A. COHEN, C. FARHAT, A. SOMACAL, AND Y. MADAY, *Nonlinear compressive reduced basis approximation for PDEs*, Hal preprint hal-04031976, (2023).

[6] J. DELON AND A. DESOLNEUX, *A Wasserstein-Type distance in the space of gaussian mixture models*, SIAM J. Imaging Sci., 13 (2020), pp. 936-970.

[7] R. A. DeVore, *Nonlinear approximation*, Acta numerica, 7 (1998), pp. 51-150.

[8] M.-H. DO, J. FEYDY, AND O. MULA, *Approximation and structured prediction with sparse wasserstein barycenters*, arXiv preprint arXiv:2302.05356, (2023).

[9] G. DUSSON, V. EHRLACHER, AND N. NOUAIME, *A wasserstein-type metric for generic mixture models, including location-scatter and group invariant measures*, arXiv preprint arXiv:2301.07963, (2023).

[10] V. EHRLACHER, D. LOMBARDI, O. MULA, AND F.-X. VIALARD, *Nonlinear model reduction on metric spaces. application to one-dimensional conservative PDEs in wasserstein spaces*, Esaim Math. Model. Numer. Anal., 54 (2020), pp. 2159-2197.

[11] C. A. FLOUDAS AND V. VISWESWARAN, *Quadratic optimization*, in Handbook of Global Optimization, R. Horst and P. M. Pardalos, eds., Springer US, Boston, MA, 1995, pp. 217-269.

[12] W. GANGBO AND A. ŚWIECH, *Optimal maps for the multidimensional Monge-Kantorovich problem*, Commun. Pure Appl. Math., (1998).

[13] J. HAMMERSLEY, *A non-harmonic Fourier series*, Acta Mathematica, 89 (1953), pp. 243-260.

[14] J. S. HESTHAVEN, G. ROZZA, AND B. STAMM, *Certified Reduced Basis Methods for Parametrized Partial Differential Equations*, SpringerBriefs in Mathematics, Springer, 2016.

[15] A. IOLLO AND T. TADDEI, *Mapping of coherent structures in parameterized flows by learning optimal transportation with gaussian models*, J. Comput. Phys., 471 (2022), p. 111671.

[16] R. MILANI, A. QUARTERONI, AND G. ROZZA, *Reduced basis method for linear elasticity problems with many parameters*, Comput. Methods Appl. Mech. Eng., 197 (2008), pp. 4812-4829.

[17] M. NONINO, F. BALLARIN, G. ROZZA, AND Y. MADAY, *Overcoming slowly decaying kolmogorov n-width by transport maps: application to model order reduction of fluid dynamics and fluid-structure interaction problems*, arXiv preprint arXiv:1911.06598, (2019).

[18] M. OHLBERGER AND S. RAVE, *Reduced basis methods: Success, limitations and future challenges*, arXiv preprint arXiv:1511.02021, (2015).

[19] D. H. PHAM, *Galerkin method using optimized wavelet-Gaussian mixed bases for electronic structure calculations in quantum chemistry*, PhD thesis, Université Grenoble Alpes, June 2017.

[20] A. QUARTERONI, A. MANZONI, AND F. NEGRI, *Reduced Basis Methods for Partial Differential Equations: An Introduction*, Springer, Aug. 2015.

[21] A. QUARTERONI, G. ROZZA, AND A. MANZONI, *Certified reduced basis approximation for parametrized partial differential equations and applications*, J. Math. Ind., 1 (2011), p. 3.

[22] D. RIM, B. PEHERSTORFER, AND K. T. MANDLI, *Manifold approximations via transported subspaces: Model reduction for transport-dominated problems*, SIAM Journal on Scientific Computing, 45 (2023), pp. A170-A199.

[23] F. ROMOR, G. STABILE, AND G. ROZZA, *Non-linear manifold Reduced-Order models with convolutional autoencoders and reduced Over-Collocation method*, J. Sci. Comput., 94 (2023), p. 74.

[24] G. ROZZA, C. N. NGUYEN, A. T. PATERA, AND S. DEPARIS, *Reduced basis methods and a posteriori error estimators for heat transfer problems*, ASME 2009 Heat Transfer Summer Conference collocated with the InterPACK09 and 3rd Energy Sustainability Conferences, (2010), pp. 753-762.

[25] V. N. TEMLYAKOV, *Nonlinear kolmogorov widths*, Mathematical Notes, 63 (1998), pp. 785-795.

[26] S. VERBLUNSKY, *On the roots of a transcendental equation, occurring in the theory of trigonometric series*, Math. Z., 61 (1954), pp. 324-335.