# Nonlinear compressive reduced basis approximation for PDE's

Albert Cohen, Charbel Farhat, Agustín Somacal, Yvon Maday

► To cite this version:

Albert Cohen, Charbel Farhat, Agustín Somacal, Yvon Maday. Nonlinear compressive reduced basis approximation for PDE's. 2023. hal-04031976

HAL Id: hal-04031976

https://hal.science/hal-04031976

Preprint submitted on 16 Mar 2023

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 compressive reduced basis approximation for PDE's

A. Cohen, C. Farhat, Y. Maday, A. Somacal

March 16, 2023

## Abstract

Linear model reduction techniques design offline low-dimensional subspaces that are tailored to the approximation of solutions to a parameterized partial differential equation, for the purpose of fast online numerical simulations. These methods, such as the POD or Reduced Basis (RB) methods, are very effective when the family of solutions has fast-decaying Karhunen-Loève eigenvalues or Kolmogorov widths, reflecting the approximability by finite-dimensional linear spaces. On the other hand, they become ineffective when these quantities have a slow decay, in particular for families of solutions to hyperbolic transport equations with parameter-dependent shock positions. The objective of this work is to explore the ability of nonlinear model reduction to circumvent this particular situation. To this end, we first describe particular notions of non-linear widths that have a substantially faster decay for the aforementioned families. Then, we discuss a systematic approach for achieving better performance via a non-linear reconstruction from the first coordinates of a linear reduced model approximation, thus allowing us to stay in the same “classical” framework of projection-based model reduction. We analyze the approach and report on its performance for a simple and yet instructive univariate test case.

## 1 Introduction

The approximation of the solution of a parameterized partial differential equation (PDE) : given $\mu$, find $u$ solution to

$$ \mathcal{D}(u; \mu) = 0 $$

can benefit from the a priori analysis of the set of all generated solutions when the parameter $\mu$ is varied, that is,

$$ \mathcal{K} := \{u_{\mu} : \mu \in \mathcal{P}\}, $$

where $u_{\mu}$ is the solution for the given value $\mu = (\mu_1, \dots, \mu_d)$ of the parameter vector that ranges in some set $\mathcal{P} \subset \mathbb{R}^d$. The set $\mathcal{K}$ is also referred to as the *solution manifold*, since it may be thought of as a parameterized $d$-dimensional manifold typically immersed in a Hilbert space $X$, where the solution to the PDE is well defined. In what follows, the norm in $X$ and the scalar product are respectively denoted by $\|\cdot\|_X$ and $\langle\cdot, \cdot\rangle_X$.

Assuming $\mathcal{K}$ to be compact in $X$, its *Kolmogorov m-width* defined as

$$ d_m(\mathcal{K})_X = \inf_{\dim(X_m) \le m} \max_{v \in \mathcal{K}} \min_{w \in X_m} \|v - w\|_X, \quad (1) $$

describes how well the set can be approximated by an ideally selected (and usually out of reach) $m$-dimensional space. If $d_m$ has a certain rate of decay as $m \to \infty$, it is possible to practically construct low-dimensional spaces $X_m$ that perform with the same approximation rate, by pre-computing offline a reduced basis consisting of $m$ solutions associated with a well-chosen set of parameters. If $d_m$ has a fast rate of decay, the RB method yields an approximation of the solution for any parameter based on an algebraic system involving very few unknowns. We refer to [QMN16], [HRS15] or [MP20] for a presentation of RB methods.

If the parameters $\mu_i$ are considered as random variables and thus $u_{\mu}$ is an $X$-valued random variable, a stochastic counterpart to these concepts is described by the principal component analysis in the Hilbert space $X$, that is, the spectral analysis of the covariance operator

$$ v \mapsto \mathbb{E}(\langle v, u_{\mu} \rangle_X u_{\mu}), $$
once $u_\mu$ has been recentered so that $E(u_\mu) = 0$. Denoting by $\sigma_1 \ge \sigma_2 \ge \dots$ the sequence of positive eigenvalues of this compact self-adjoint operator, and by $e_1, e_2, \dots$ the Karhunen-Loève orthonormal basis of eigenfunctions, it is well known that

$$ \kappa_m^2 := \min_{\dim(X_m) \le m} \mathbb{E} \left( \min_{w \in X_m} \|u_\mu - w\|_X^2 \right) = \sum_{n>m} \sigma_n, $$

and that the minimizing space is spanned by $e_1, \dots, e_m$. This is the starting point to the proper orthogonal decomposition (POD) method which amounts to replacing the aforementioned eigenfunctions by approximations computed offline, based on a sufficiently large set of training solutions.

These linear reduced modeling methods have penetrated industrial applications, a guarantee of their success. However, there are still cases where these approaches have difficulties to overcome, namely, when the Kolmogorov width $d_m$ or eigenvalues $\sigma_m$ do not decay fast.

This is in particular what happens for transport type problems. Even in the conceptually simple case of constant speed translation of an initial condition given by a step function, where the only parameter is the position of the discontinuity, it is well known that with $X = L^2$, the numbers $d_m$ and $\kappa_m$ decay slowly like $O(m^{-1/2})$. In other words, for a target precision of $\varepsilon$, the basis is of prohibitive size $O(\varepsilon^{-2})$.

For families of such functions, substantial gain can be expected when searching for nonlinear reduced models. Prominent examples of nonlinear approximation include rational fractions, finite elements on adaptive grids of fixed cardinality, n-term approximations in a basis or dictionary, and neural networks, see [DeV98, DHP21] for a general treatment. In these methods, the “coordinates” describing the approximation to a function $u$ are typically nonlinear functionals applied to $u$, and the reconstruction map from such parameters is also nonlinear. In the frame of model reduction, we refer to [AZF12, GFTBM21, BCD$^{+}$21] that considers libraries of affine reduced models, [BF22] that uses quadratic manifolds, and [LC20, FDM21, GGJW22, PDG$^{+}$22, PWL23, BFM23] for neural network based approaches, see also [BIR18, CMS19, BSU20, HMC22], and [Peh22] for an overview on these nonlinear approaches.

Interestingly, it appears that an efficient approach to nonlinear model reduction is to maintain linear functionals for computing the coordinates while performing reconstruction in a well-chosen nonlinear way. This state of affair is in particular illustrated by the development of compressed sensing in the last two decades, where signals are reconstructed from linear measurements by nonlinear methods promoting sparsity, such as $\ell^1$ minimization.

In this note, we begin by substantiating this idea more precisely in §2, by recalling and comparing certain notions of linear and nonlinear m-widths. We present in §3 a general approach that consists in taking as linear functionals the first components in a linear reduced model (RB or POD) that has been learned offline; and also use the offline stage to learn a computationally tractable nonlinear map that reconstructs the missing components from these first ones to reach a better accuracy. One key aspect lies in the type of nonlinear maps that is allowed. This approach is analyzed in §4, in the case of a simple univariate model of step functions; it is illustrated by numerical tests for this model in §5.

## 2 Linear and nonlinear notions of m-widths

Generally speaking, the process of dimensionality reduction can be described by a pair of continuous mappings, the encoder

$$ E: X \to \mathbb{R}^m, $$

and the decoder

$$ D: \mathbb{R}^m \to X. $$

The maximum distorsion of the encoding procedure over $\mathcal{K}$ is given by the quantity

$$ \max_{v \in \mathcal{K}} \|v - D(E(v))\|_X. $$

Then, for a general Banach space $X$ and a compact set $\mathcal{K} \subset X$, we can define various notion of widths

$$ \inf_{D, E} \max_{v \in \mathcal{K}} \|v - D(E(v))\|_X, $$

by optimizing the choice of $E$ and $D$, under specific restrictions:
• If $D$ and $E$ are both assumed to be linear, one obtains the approximation numbers

$$a_m(\mathcal{K})_X := \inf_L \max_{v \in \mathcal{K}} \|v - Lv\|_X,$$

where the infimum is taken over operators $L$ of rank at most $m$.

• If only $D$ is assumed to be linear, one obtains the already mentioned Kolmogorov width $d_m(\mathcal{K})_X$. When $X$ is a general Banach space, the inequality

$$d_m(\mathcal{K})_X \le a_m(\mathcal{K})_X,$$

can be strict. Equality obviously holds in the case when $X$ is a Hilbert space since best approximation in a subspace of $X$ of dimension $m$ is achieved by linear orthogonal projection.

• The sensing numbers $s_m(\mathcal{K})_X$ correspond to the reciprocal situation, where $E$ is assumed to be linear and $D$ is assumed to be nonlinear. In other words, they can be defined as

$$s_m(\mathcal{K})_X := \inf_{D, \lambda_1, \dots, \lambda_m} \max_{v \in \mathcal{K}} \|v - D(\lambda_1(v), \dots, \lambda_m(v))\|_X,$$

where the infimum is taken over all choice of linear functionals $\lambda_1, \dots, \lambda_m \in X'$ and decoding map $D$. These number are closely related to the Gelfand width classically defined as

$$d^m(\mathcal{K})_X := \inf_{\lambda_1, \dots, \lambda_m} \max\{\|v\|_X : v \in \mathcal{K}, \lambda_1(v) = \dots = \lambda_m(v) = 0\}.$$

It is easily checked that $s_m(\mathcal{K})_X = d^m(\mathcal{K})_X$ in the case where $\mathcal{K}$ is convex and centrally symmetric; and that

$$s_m(\mathcal{K})_X \le d^m(\mathcal{K} - \mathcal{K})_X \le 2s_m(\mathcal{K})_X,$$

for a general compact set $\mathcal{K}$ and $\mathcal{K} - \mathcal{K}$ is a notation for the set $\{u : u = v - w, v \in \mathcal{K}, w \in \mathcal{K}\}$.

• Finally, the nonlinear width or manifold width $\delta_m(\mathcal{K})_X$ is defined when no other assumption but continuity is made on the operators $E$ and $D$. For numerical stability purpose, it is interesting to tame this notion by imposing that $D$ and $E$ are both Lipschitz continuous, that is

$$\|D(a) - D(b)\|_X \le \gamma \|a - b\|_m \quad \text{and} \quad \|E(v) - E(w)\|_m \le \gamma \|v - w\|_X, \quad a, b \in \mathbb{R}^m, v, w \in X,$$

for some fixed $\gamma > 1$, with $\| \cdot \|_m$ an arbitrary norm on $\mathbb{R}^m$. The resulting infimized quantity $\delta_m(\mathcal{K})_X$ is referred to as the stable width.

The last two notions of width $s_m$ and $\delta_m$ (or $\delta_m^\gamma$) are natural to describe the expected performance of optimal nonlinear model reduction, since the manifold is approximated by the set $D(\mathbb{R}^m)$ – which is no longer a linear space. However, the sensing numbers take the view that encoding can be restricted to simple linear measurements.

As already mentioned, the quantities $d_m$ and $a_m$ typically decay slowly for families of piecewise smooth functions, which reflects the fact that they cannot be well approximated efficiently by linear spaces. A substantial gain in the rate of decay can be expected however when considering the nonlinear widths $\delta_m$ and $\delta_m^\gamma$. Interestingly, it appears that this substantial optimal gain is already present when considering the sensing numbers $s_m$.

As a basic example, consider the two-parameter family of univariate step functions

$$\mathcal{K} := \{u := \chi_{[a, a+\ell]} : a \in \mathbb{R}, \ell > 0\}.$$

Clearly, the parameters $(a, \ell)$ are not linear functionals of $u$. However, any $u \in \mathcal{K}$ can be exactly reconstructed from two linear functionals, namely, the first moments

$$\lambda_k(u) = \int x^k u(x) dx, \quad k = 0, 1.$$

Indeed, $\lambda_0 = \ell$ and $\lambda_1 = \frac{1}{2}\ell(2a+\ell)$, so that $a$ and $\ell$ can be exactly recovered from such data. Therefore, one has $s_m(\mathcal{K})_X = 0$ for any $m > 2$ and for any Banach space $X$.
At a more general level, it was proved in [CDPW22] that when $X$ is a Hilbert space, then both $s_m(K)_X$ and $\delta_m^\gamma(K)_X$ are tied to the so-called entropy numbers $\varepsilon_m(K)_X$ defined as the smallest value of $\varepsilon > 0$ such that $K$ can be covered by $2^m$ balls of radius $\varepsilon$. More precisely, it was shown that, on the one hand, for any $s > 0$, one has a Carl-type inequality

$$ \sup_{m>0} m^s \varepsilon_m(K)_X \le C_s \sup_{m>0} m^s \delta_m^\gamma(K)_X, $$

where $C_s$ depends on $(s, \gamma)$, and that on the other hand, there exists $c > 0$ depending on $\gamma > 1$ such that

$$ \delta_{cm}^\gamma(K)_X \le 3\varepsilon_m(K)_X, \quad m \ge 1. $$

In the proof of this last inequality, the $\gamma$-stable encoding-decoding pair $(E, D)$ which is constructed has actually a linear $E$. In turn, one also has

$$ s_{cm}(K)_X \le 3\varepsilon_m(K)_X, \quad m \ge 1. $$

One consequence of these results is that $s_m(K)_X$, $\delta_m^\gamma(K)_X$ and $\varepsilon_m(K)_X$ share the same algebraic rates of decay.

**Remark 1.** An additional aspect of non-linear dimensionality reduction is the notion of adaptivity, which means that the measurements $E(u) = (E_1(u), \dots, E_m(u))$ are chosen incrementally, that is, the functional $E_m$ is picked depending on the value of $E_1(u), \dots, E_{m-1}(u)$. This allows the definition of similar notions of adaptive sensing numbers and non-linear widths. Our next described approach is not of this form, since we use linear functionals that are pre-defined through the standard POD or RB analysis.

## 3 Non-linear compressive Reduced Basis approximation

In this contribution, we thus intend to deal with situations where:

* The Kolmogorov widths $d_m(K)_X$, or the singular values $\sigma_n$, decay slowly.

* The sensing numbers $s_m(K)_X$, and stable non-linear widths $\delta_m^\gamma(K)_X$, decay much faster.

In other words, a target accuracy $\varepsilon > 0$ can be reached by $d_N(K)_X$ or $\kappa_N$, however with a dimension $N = N(\varepsilon)$ much larger than the value of $n = n(\varepsilon)$, such that $s_n(K)_X$ reaches the same accuracy.

Since the optimal linear functionals in the definition of $s_n(K)_X$ are usually out of reach and could be computationally unpractical to apply, **we take the view of fixing these measurements to be a small number $n$ of components in the offline computed (orthonormalized) RB or POD basis** $(e_j)_{j=1,...,N}$ for some $N \gg n$. Typically, we choose the $n$ first ones, that is,

$$ \lambda_j(v) = \langle v, e_j \rangle_X, \quad j = 1, \dots, n. $$

Intuitively, it is expected that in the situation where $s_n(K)_X$ is very small, then the unknown component $(\lambda_j(v))_{j=n+1,...,N}$ should be somehow dependent, up to a small error, of the $n$ first ones that carry most of the relevant information. This idea was at first presented in [BFM23]. Here, we formalize it and study its validity in detail on a simple step function model, and propose a general numerical strategy that we test on this model.

Our objective is thus to predict from these first components the extra components $\lambda_k(v)$ for $k = n+1, \dots, N$ that are needed to approximate the functions $v \in K$ with target accuracy $\varepsilon$. We are thus interested to construct $N-n$ functions $\psi_k: \mathbb{R}^n \to \mathbb{R}$ so that

$$ \tilde{\lambda}_k(v) := \psi_k(\lambda_1(v), \dots, \lambda_n(v)), $$

is a very accurate approximation to $\lambda_k(v)$ for $k = n+1, \dots, N$ and can be fastly computed.

Let us stress that $\psi_k$ should typically be a nonlinear function. Indeed consider the ideal case of the PCA basis computed after having recentered the variable $u_\mu$. Then the variables

$$ z_j = \lambda_j(u_\mu), $$
are uncorrelated and centered, such that, for any $k > n$,

$$ \min_{\alpha_1, \dots, \alpha_n} \mathbb{E}(|z_k - \sum_{j=1}^{n} \alpha_j z_j|^2) = \mathbb{E}(|z_k|^2). $$

Thus, the best choice of a linear function would be the null one that does not deliver any information.

On the other hand, the best choice of a nonlinear function in this mean square sense, that is, minimizing $\mathbb{E}(|z_k - \psi(z_1, \dots, z_n)|^2)$ over all functions $\psi$, is given by the conditional expectation

$$ \psi_k^*(z_1, \dots, z_n) = \mathbb{E}(z_k | z_1, \dots, z_n), $$

which is out of reach and should be approximated by a computationally tractable function.

Our practical approach to the construction of $\psi_k$ is by learning it in a second step of the offline stage, after the basis $(e_j)_{j=1,\dots,N}$ has been identified. Having in mind the above mean square loss, one typical approach is to select $\psi_k$ within a sufficiently rich class $\mathcal{F}$ of nonlinear functions by empirical risk minimization : with $(u^i)_{i=1,\dots,M}$ a training set of random snapshots $u^i = u_{\mu^i}$, we define

$$ \psi_k := \underset{\lambda_k(u^i) \in \mathcal{F}}{\operatorname{argmin}} \left\{ \sum_{i=1}^{m} |\lambda_k(u^i) - \psi(\lambda_1(u^i), \dots, \lambda_n(u^i))|^2 : \psi \in \mathcal{F} \right\}. $$

A critical aspect in this approach lies in the choice of the class $\mathcal{F}$, which could be, for example, the set of:

* Quadratic functions, as in [BF22] or [GWW23].

* Polynomials of some higher degree $d > 2$.

* Neural networks with a given architecture, as proposed in [BFM23] (see also [LC20] where an autoencoder-based approximation was proposed, which was in a way a pioneering idea but unfortunately not computationally tractable one).

This class should be able to approximate correctly the ideal but out of reach $\psi_k^*$ by a computationally tractable function $\psi_k \in \mathcal{F}$. Another difficulty with this approach is the fact that when the number $n$ of informative components is chosen to be not very small, one faces a regressing problem in large dimension, for which classical methods such as splines or polynomials are known to suffer from the curse of dimensionality.

For these reasons, we have also considered in our numerical tests regression methods based on trees (CART) and random forests, that are both universally consistent and able to tackle large-dimensional problems. These methods seem to deliver the best numerical results for the considered problems.

# 4 Analysis of a model framework : periodic step functions

In order to investigate the aforementioned questions, we place ourselves in a framework where the Karhunen-Loève basis is explicitely known. Specifically, we work in the Hilbert space

$$ X = L^2(0, 1), $$

and consider a randomly parameterized family such that

$$ \mathbb{E}(u_{\mu}(x)) = \bar{u}, $$

independently of $x \in [0, 1]$ and such that

$$ \mathbb{E}((u_{\mu}(x) - \bar{u})(u_{\mu}(y) - \bar{u})) = R(x - y), $$

where $R$ is an even and 1-periodic function. In other words, $u_{\mu}$ is a periodic stationary process, its covariance operator coincides with the convolution operator by $R$, and therefore its Karhunen-Loève basis is exactly given by the basis of the Fourier series on $[0, 1]$ (see e.g. [PUP02]).
More specifically, we consider a simple model of periodic stationary step functions by introducing the three-parameter family

$$u_{\mu}(x) := \begin{cases} b & \text{for } x \in (a, a + \ell) \pmod{1} \\ 0 & \text{for } x \in (a + \ell, a) \pmod{1} \end{cases}, \quad \mu = (a, \ell, b), \qquad (2)$$

that is, $u_\mu = b\chi_{[a,a+\ell]}$ in a 1-periodic sense.

Here $a, \ell$, and $b$ are assumed to have independent uniform distributions. Taking the base point $a$ to be uniformly distributed over $[0, 1]$, it is easily checked that the process is periodic stationary. In addition, we take the height $b$ to be uniformly distributed in $[0, 1]$ and the length $\ell$ to be uniformly distributed in $[\ell_{\min}, 1 - \ell_{\min}]$ for some $0 < \ell_{\min} < \frac{1}{2}$.

The best linear approximation of dimension $m = 2n+1$ is thus given by the truncation up to $k \le n$ of the Fourier expansion

$$u_{\mu} = \sum_{k \in \mathbb{N}} \alpha_k \cos(2\pi k x) + \sum_{k \in \mathbb{N}^*} \beta_k \sin(2\pi k x) \qquad (3)$$

where

$$\left\{ \begin{array}{lcl} \alpha_0 = \alpha_0(a, \ell, b) & = & b\ell \\ \alpha_k = \alpha_k(a, \ell, b) & = & b\frac{\sin(2\pi k(a+\ell)) - \sin(2\pi ka)}{2\pi k} = b\frac{\sin(\pi k\ell)\cos(\pi k(2a+\ell))}{\pi k} \\ \beta_k = \beta_k(a, \ell, b) & = & b\frac{\cos(2\pi k(a+\ell)) - \cos(2\pi ka)}{2\pi k} = -b\frac{\sin(\pi k\ell)\sin(\pi k(2a+\ell))}{\pi k} \end{array} \right. \qquad (4)$$

Clearly $\sigma_0 = \mathbb{E}(|\alpha_0|^2) = \mathbb{E}(\ell^2)\mathbb{E}(b^2) = \frac{1}{9}((1 - \ell_{\min})^3 - \ell_{\min}^3)$. It is also easily checked that the eigenvalues associated to the functions $x \mapsto \cos(2\pi kx)$ and $x \mapsto \sin(2\pi kx)$ are the same and are given by

$$\sigma_k = \mathbb{E}(|\alpha_k|^2) = \mathbb{E}(|\beta_k|^2) = \frac{c}{k^2},$$

for some $c = c(\ell_{\min}) > 0$. It follows that the best linear approximation has a mean-square error $\kappa_m^2$ behaving like $m^{-1}$. Note that, for the corresponding manifold

$$\mathcal{K} := \{u_{\mu} : a \in [0, 1], \ell \in [\ell_{\min}, 1 - \ell_{\min}], b \in [0, 1]\},$$

one obviously has $d_m(\mathcal{K})_X \ge \kappa_m$, since a worst case error dominates the average error. On the other hand, it is also readily seen that the worst case approximation by Fourier series behaves like $m^{-1/2}$, and therefore

$$d_m(\mathcal{K})_X \sim \kappa_m \sim m^{-1/2}.$$

We also consider the two-parameter family $\mathcal{K}_{\ell_0}$ obtained by freezing the value $\ell = \ell_0$ and the one-parameter family $\mathcal{K}_{b_0, \ell_0}$ obtained by freezing in addition the value $b = b_0$. It is easily checked that one has the same behaviour $m^{-1/2}$ for $\kappa_m$ and $d_m$ after such restrictions.

The one-parameter family $\mathcal{K}_{b_0, \ell_0}$ can be encoded by the data of $a$, so that its nonlinear width satisfies

$$\delta_m(\mathcal{K}_{b_0, \ell_0})_X = 0, \quad m > 1.$$

It is easily seen that the data of only one Fourier coefficient is not sufficient to characterize the elements of this family. Indeed, $\alpha_0(a, \ell, b)$ is independent of $a$, $\alpha_k(a, \ell, b) = \alpha_k(1/2 - a - \ell, \ell, b)$, and $\beta_k(a, \ell, b) = \beta_k(1/4 - a - \ell, \ell, b)$ for $k \ne 0$.

On the other hand, the recovery of $a$ can be done through the data of the two coefficients $\alpha_1$ and $\beta_1$ since

$$a = -\frac{\ell}{2} - \frac{1}{2\pi} \arctan \left( \frac{\beta_1}{\alpha_1} \right) \pmod{1} \qquad (5)$$

Similarly, any element in the two-parameter family $\mathcal{K}_{\ell_0}$, parametrized by $a$ and $b$, can be recovered from the data of these two coefficients, since one also has

$$b = \frac{\pi}{\sin(\pi \ell_0)} (\alpha_1^2 + \beta_1^2)^{1/2}.$$

When more coefficients are available, we note that there is not a unique reconstruction map : for example from the three coefficients $\alpha_0, \alpha_1$, and $\beta_1$, we can also recover $b$ according to

$$b = \frac{\alpha_0}{\ell_0}.$$
Finally, in the case of the three parameter family $\mathcal{K}$, exact recovery of $(a, b, \ell)$ can be obtained by solving the nonlinear system

$$ \left\{ \begin{array}{l} a + \frac{\ell}{2} = -\frac{1}{\pi}\arctan\left(\frac{\beta_1}{\alpha_1}\right) \\ \pi b \sin(\pi \ell) = \left(\alpha_1^2 + \beta_1^2\right)^{1/2} \\ b \ell = \alpha_0 \end{array} \right. \qquad (6) $$

however the exact recovery map does not anymore have an explicit form.

These exact recovery procedures induce for all $k > 1$ an exact recovery map $\psi_k^*$ such that

$$ \alpha_k = \psi_k^*(\alpha_0, \alpha_1, \beta_1). $$

and similarly an exact recovery map $\tilde{\psi}_k^*$ such that

$$ \beta_k = \tilde{\psi}_k^*(\alpha_0, \alpha_1, \beta_1), $$

In this very simple case, the success of the learning strategy outlined in the previous section therefore depends on how these maps can be approximated by the family $\mathcal{F}$.

A simple intuition can be given when looking at the particular case of $\tilde{\psi}_k^*$ for the one-parameter family $\mathcal{K}_{b_0, \ell_0}$, when $b_0 = 1$ and $\ell_0 = \frac{1}{2}$. Then, we find that

$$ \beta_k = b \frac{\sin(\pi k/2) \cos(2\pi k a)}{\pi k}, $$

which is null for even values of $k$, and for odd values $k = 2j + 1$ satisfies

$$ \beta_k = \frac{b}{k\pi} (-1)^j T_k\left(\frac{\pi}{b}\beta_1\right), $$

where $T_k$ is the Chebychev polynomial of degree $k$. Therefore for such values, the optimal reconstruction is exact and given by

$$ \tilde{\psi}_k^*(x, y, z) = \frac{b}{k\pi} (-1)^j T_k\left(\frac{\pi}{b}z\right). $$

Clearly, the function $\tilde{\psi}_k^*$ cannot be well approximated by polynomials of moderate dimensions for large values of $k$. On the other hand, it is well known that the derivative of $T_k$ has maximal norm of order $k$ over $[-1, 1]$, and this implies that the functions $\tilde{\psi}_k^*$ are Lipschitz continuous with Lipschitz constant bounded independently of $k > 0$.

This property holds in more generality from the following argument: the derivative of the arctan function being upper bounded by 1, the recovery of $a$ from $\alpha_1(a, \ell_0, b_0)$ and $\beta_1(a, \ell_0, b_0)$ is stable as

$$ \frac{da}{d\alpha_1} = \frac{1}{2\pi} \frac{1}{1 + \left[\frac{\beta_1}{\alpha_1}\right]^2} \frac{\beta_1}{\alpha_1^2} = -\frac{1}{2\pi} \frac{\beta_1}{\alpha_1^2 + \beta_1^2} \qquad (7) $$

and

$$ \frac{da}{d\beta_1} = -\frac{1}{2\pi} \frac{1}{1 + \left[\frac{\beta_1}{\alpha_1}\right]^2} \frac{1}{\alpha_1^2} = -\frac{1}{2\pi} \frac{\alpha_1}{\alpha_1^2 + \beta_1^2} \qquad (8) $$

are both bounded since, by construction,

$$ [\alpha_1]^2 + [\beta_1]^2 = b_0^2 \frac{\sin^2(\pi \ell_0)}{\pi^2}. $$

Hence, an error in the values of $\alpha_1$ or $\beta_1$ will induce an error of comparable size on $a$. The same hold for the determination of $b$ in the case of $\mathcal{K}_{\ell_0}$.

On the other hand, it is readily seen from the definition of Fourier coefficients that the two maps

$$ \mu \mapsto \alpha_k(u_\mu) \quad \text{and} \quad \mu \mapsto \beta_k(u_\mu), $$

are Lipschitz continuous with Lipschitz constants bounded independently of $k > 0$. In turn, the stable recovery of $a$ and $b$ from $\alpha_0, \alpha_1$ and $\beta_1$, induces recovery maps

$$ (\alpha_1(a, \ell_0, b_0), \beta_1(a, \ell_0, b_0)) \mapsto (\alpha_k(a, \ell_0, b_0), \beta_k(a, \ell_0, b_0)) \qquad (9) $$
and

$$ (\alpha_0(a, \ell_0, b), \alpha_1(a, \ell_0, b), \beta_1(a, \ell_0, b)) \mapsto (\alpha_k(a, \ell_0, b), \beta_k(a, \ell_0, b)) \quad (10) $$

that are Lipschitz continuous with Lipschitz constants bounded independently of $k > 0$.

This state of affairs explains that universally consistent methods such as random forests are well adapted for the joint approximation of $\psi_k^*$ and $\tilde{\psi}_k^*$, while approaches based on low order polynomials are doomed to fail. This is confirmed by the numerical tests presented in the next section.

In the perspective of recovering more general piecewise smooth functions, we expect that the low-order components are affected by the smooth pieces in addition to the jumps, while the high-order components are only affected by the jumps. Thus it is interesting to address the question of the recovery of the parameters of the step function from a few components $\alpha_k$ and $\beta_k$ for larger values of $k$.

This task appears to be more involved and requiring more coefficients. For example, when recovering $a$ as in (5), we find

$$ a = -\frac{\ell}{2} - \frac{1}{2\pi} \arctan \left[ \frac{\beta_k}{\alpha_k} \right] \pmod{1/k}. \quad (11) $$

One possibility to lift the indeterminacy (mod $1/k$) is to combine the information coming from (11) and

$$ a = -\frac{\ell}{2} - \frac{1}{2\pi} \arctan \left[ \frac{\beta_{k+1}}{\alpha_{k+1}} \right] \pmod{1/(k+1)} \quad (12) $$

since $a = a'$, (mod $1/k$), $a = a'$, (mod $1/(k+1)$), and $a, a' \in (0, 1)$ imply $a = a'$.

Thus we can in principle recover the parameters $a$ and $b$ out of 4 coefficients of arbitrary high frequencies ($k, k+1$). However we also observe that the stability of this recovery is deteriorated since $\frac{da}{d\alpha_k}$ and $\frac{da}{d\beta_k}$ increase linearly with $k$. We may hope to improve the stability by using a larger number of coefficient values with indexes ($k, k+1, \dots, k+d$) for some $d > 1$.

# 5 Numerical illustrations

In this section, we investigate the ability of different methods to learn mappings that use different amount $m$ of components, namely

$$ m = 2, \quad (\alpha_1, \beta_1) \quad \mapsto \quad (\alpha_k, \beta_k) $$

$$ m = 3, \quad (\alpha_0, \alpha_1, \beta_1) \quad \mapsto \quad (\alpha_k, \beta_k) $$

$$ m = 5, \quad (\alpha_0, \alpha_1, \beta_1, \alpha_2, \beta_2) \quad \mapsto \quad (\alpha_k, \beta_k) $$

for each of the three families $K_{b_0, \ell_0}$ (Figures 1 and 2), $K_{\ell_0}$ (Figures 3 and 4), $K$ (Figures 5 and 6), for all $2 \le k \le 500$. In these Figures, the average recovery error for the $\alpha_k$ and $\beta_k$ are presented in a symmetric manner, on the left and right side of the x-axis respectively.

The learning methods are

* linear regression: $\mathcal{F}$ is the set of linear functions.

* quadratic regression: $\mathcal{F}$ is the set of polynomials of total degree 2.

* quartic regression: $\mathcal{F}$ is the set of polynomials of total degree 4.

* decision tree [HTF09]

* random forest [Bre01] [GEW06]

For all the numerical illustrations we used Python 3.8 and scikit-learn 1.2.2 [PVG+11] for the implementation of each of the regression methods described above. For more information, the code can be found at https://github.com/agussomacal/NonLinearRBA4PDEs

As can be expected, linear regression give the same results as the null forecast, and quadratic and quartic regression give the same (bad) result here, with some improvement over the null forecast only for very small value of $k$ in certain cases (see Figures 3 and 4).

In contrast, decision tree and random forest are well suited as can be seen for the one parameter family on Figure 1, with improved results on Figure 2 obtained from a larger training samples (10 000
rather than 1 000). This reflects the universal consistency of these methods that are guaranteed to converge towards the regression function as the number of sample tends to $+\infty$.

The same also holds for the two parameter family, as seen on Figures 3 and 4: the problem is slightly more involved but nevertheless decision tree and random forest manage to obtain a fair (resp. good) approximation after a learning phase of 1 000 (resp. 10 000) training samples.

<figure><img src="image_2.png" /><figcaption>Figure 1: In this figure we plot the error obtained from different recovery methods for the family $K_{b=1, \ell=.5}$ where we recover all coefficients $\alpha_k$ and $\beta_k$ in (3) for $2 \le k \le 500$ from 2 (left) 3 (center) and 5 (right) Fourier coefficients with different approaches: linear, quadratic, quartic, tree and random forests. Note that linear, quadratic, quartic are superposed and do not improve over the trivial recovery of the missing modes by value 0. The learning phase is based on 1 000 training samples. The x-axis represents (in a log scale) the index $k$ of $\alpha_k$ and the index $k$ of $\beta_k$ and the y-axis the mean-square reconstruction error on the mode.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2: Same test as Figure 1 with 10 000 training samples.</figcaption></figure>

The numerical results for the three parameter family are displayed on Figure 5 for the range $\ell \in [0.4, 0.6]$, and on Figure 6 for the range $\ell \in [0.01, 0.99]$, that is $\ell_{\min} = 0.4$ and $0.01$ respectively. One first observation is that all methods fail in the case of $m = 2$ known component since they are insufficient to characterize an element of $K$. Secondly we observe that the performance deteriorates as $\ell$ is allowed to be very small in which case the exact recovery becomes less stable in view of the multiplicative coupling between $b$ and $\ell$ in the last two equations of the system (6). These last results reveal difficulties for these too simple non-linear recovery methods to achieve a satisfactory recovery. We expect that more involved approaches such as deep neural networks can improve this state of affair.

Finally, in the case of the two parameter family ($\ell$ fixed), we study the recovery error of $(\alpha_k, \beta_k)$ for $k > 10 + d$ when using information from high frequency coefficients $(\alpha_j, \beta_j)$ for $j = 10, 11, \dots, 10 + d$. As explained in the end of the last section, the exact recovery is feasible for $d = 1$, yet less stable and thus more difficult to learn. This is confirmed by Figure 7, where we see an improvement when using a larger value of $d$ and a larger training sample.
<figure><img src="image_4.png" /><figcaption>Figure 3: Same test as Figure 1 for the two parameter family $\mathcal{K}_{\ell=.5}$, using 1 000 training samples.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: Same test as Figure 3 for the two parameter family $\mathcal{K}_{\ell=.5}$, using 10 000 training samples.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: Same test as Figure 1 for the three parameter family $\mathcal{K}$, using 10 000 training and $\ell_{\min} = 0.4$.</figcaption></figure>

## References

[AZF12] David Amsallem, Matthew J Zahr, and Charbel Farhat. Nonlinear model order reduction based on local reduced-order bases. *International Journal for Numerical Methods in Engineering*, 92(10):891–916, https://doi.org/10.1002/nme.4371, 2012.

[BCD+21] Andrea Bonito, Albert Cohen, Ronald Devore, Diane Guignard, Peter Jantsch, and Guergana Petrova. Nonlinear methods for model reduction. *ESAIM: Mathematical Modelling and Numerical Analysis*, 55(2):507–531, March 2021, https://doi.org/10.1051/m2an/2020057.
<figure><img src="image_7.png" /><figcaption>Figure 6: Same test as Figure 1 for the three parameter family $\mathcal{K}$, using 10 000 samples and $\ell_{\min} = 0.01$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: Recovery of components $k$ for $k > 10 + d$ for the two parameter family $\mathcal{K}_{\ell=5}$ using random forest and components $j$ for $j = 10, 11, \dots, 10 + d$ with $d = 4$ (left) and $d = 9$ (right), and various numbers of training samples.</figcaption></figure>

[BF22] Joshua Barnett and Charbel Farhat. Quadratic approximation manifold for mitigating the kolmogorov barrier in nonlinear projection-based model order reduction. *Journal of Computational Physics*, 464:111348, https://doi.org/10.1016/j.jcp.2022.111348, 2022.

[BFM23] Joshua L Barnett, Charbel Farhat, and Yvon Maday. Mitigating the kolmogorov barrier for the reduction of aerodynamic models using neural-network-augmented reduced-order models. In *AIAA SCITECH 2023 Forum*, page 0535, https://doi.org/10.2514/6.2023-0535, 2023.

[BIR18] Florian Bernard, Angelo Iollo, and Sébastien Riffaud. Reduced-order model for the bgk equation based on pod and optimal transport. *Journal of Computational Physics*, 373:545-570, https://doi.org/10.1016/j.jcp.2018.07.001, 2018.

[Bre01] Leo Breiman. Random forests. *Machine learning*, 45:5-32, 2001.

[BSU20] Felix Black, Philipp Schulze, and Benjamin Unger. Projection-based model reduction with dynamically transformed modes. *ESAIM: Mathematical Modelling and Numerical Analysis*, 54(6):2011-2043, https://doi.org/10.48550/arXiv.1912.11138, 2020.

[CDPW22] Albert Cohen, Ronald DeVore, Guergana Petrova, and Przemysław Wojtaszczyk. Optimal stable nonlinear approximation. *Foundations of Computational Mathematics*, 22(3):607-648, https://doi.org/10.1007/s10208-021-09494-z, 2022.
[CMS19] Nicolas Cagniart, Yvon Maday, and Benjamin Stamm. Model order reduction for problems with large convection effects. *Contributions to partial differential equations and applications*, pages 131–150, https://doi.org/10.1007/978-3-319-78325-3_10, 2019.

[DeV98] Ronald DeVore. Nonlinear approximation. *Acta numerica*, 7:51–150, https://doi.org/10.1017/S0962492900002816, 1998.

[DHP21] Ronald DeVore, Boris Hanin, and Guergana Petrova. Neural network approximation. *Acta Numerica*, 30:327–444, https://doi.org/10.1017/S0962492921000052, 2021.

[FDM21] Stefania Fresca, Luca Dede, and Andrea Manzoni. A comprehensive deep learning-based approach to reduced order modeling of nonlinear time-dependent parametrized pdes. *Journal of Scientific Computing*, 87(2):1–36, https://doi.org/10.1007/s10915-021-01462-7, 2021.

[GEW06] Pierre Geurts, Damien Ernst, and Louis Wehenkel. Extremely randomized trees. *Machine Learning*, 63:3–42, 04 https://doi.org/10.1007/s10994-006-6226-1, 2006.

[GFTBM21] Sebastian Grimberg, Charbel Farhat, Radek Tezaur, and Charbel Bou-Mosleh. Mesh sampling and weighting for the hyperreduction of nonlinear petrov-galerkin reduced-order models with local reduced-order bases. *International Journal for Numerical Methods in Engineering*, 122(7):1846–1874, https://doi.org/10.1002/nme.6603, 2021.

[GGJW22] Anthony Gruber, Max Gunzburger, Lili Ju, and Zhu Wang. A comparison of neural network architectures for data-driven reduced-order modeling. *Computer Methods in Applied Mechanics and Engineering*, 393:114764, https://doi.org/10.1016/j.cma.2022.114764, 2022.

[GWW23] Rudy Geelen, Stephen Wright, and Karen Willcox. Operator inference for non-intrusive model reduction with quadratic manifolds. *Computer Methods in Applied Mechanics and Engineering*, 403:115717, https://doi.org/10.1016/j.cma.2022.115717, 2023.

[HMC22] Willy Haik, Yvon Maday, and Ludovic Chamoin. Assimilation de données variationnelle et hybride pour le contrôle thermique en temps réel du procédé de fabrication additive. In *CSMA: 15ème Colloque National en Calcul des Structure*, https://hal.science/hal-03607001, 2022.

[HRS15] Jan S Hesthaven, Gianluigi Rozza, and Benjamin Stamm. Certified reduced basis methods for parametrized partial differential equations, volume 590. Springer, https://doi.org/10.1007/978-3-319-22470-1, 2015.

[HTF09] Trevor Hastie, Robert Tibshirani, and Jerome Friedman. *The Elements of Statistical Learning: Data Mining, Inference, and Prediction, Second Edition* (Springer Series in Statistics). 02 https://doi.org/10.1007/978-0-387-84858-7, 2009.

[LC20] Kookjin Lee and Kevin T Carlberg. Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders. *Journal of Computational Physics*, 404:108973, https://doi.org/10.48550/arXiv.1812.08373, 2020.

[MP20] Yvon Maday and A Patera. Reduced basis methods. *Model Order Reduction; Benner, P., Grivet-Talocia, S., Quarteroni, A., Rozza, G., Schilders, W., Silveira, L., Eds*, pages 139–179, https://doi.org/10.1515/9783110671490-004, 2020.

[PDG+22] Davide Papapicco, Nicola Demo, Michele Girfoglio, Giovanni Stabile, and Gianluigi Rozza. The neural network shifted-proper orthogonal decomposition: A machine learning approach for non-linear reduction of hyperbolic equations. *Computer Methods in Applied Mechanics and Engineering*, 392:114687, https://doi.org/10.1016/j.cma.2022.114687, 2022.

[Peh22] Benjamin Peherstorfer. Breaking the kolmogorov barrier with nonlinear model reduction. *Notices of the American Mathematical Society*, 69(5):725–733, https://doi.org/10.1090/noti2475, 2022.
[PUP02] Athanasios Papoulis and S Unnikrishna Pillai. *Probability, random variables and stochastic processes*. McGraw-Hill, ISBN 0-07-366011-6, 2002.

[PVG+11] Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, et al. Scikit-learn: Machine learning in python. *Journal of machine learning research*, 12(Oct):2825–2830, 2011.

[PWL23] Zhichao Peng, Min Wang, and Fengyan Li. A learning-based projection method for model order reduction of transport problems. *Journal of Computational and Applied Mathematics*, 418:114560, https://doi.org/10.1016/j.cam.2022.114560, 2023.

[QMN16] Alfio Quarteroni, Andrea Manzoni, and Federico Negri. *Reduced basis methods for partial differential equations: an introduction*, volume 92. Springer, https://doi.org/10.1007/978-3-319-15431-2, 2016.