# AMStraMGRAM: Adaptive Multi-cutoff Strategy
## Modification for ANaGRAM

Nilo Schwencke, Cyriaque Rousselot, Alena Shilova, Cyril Furtlehner

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

To cite this version:

Nilo Schwencke, Cyriaque Rousselot, Alena Shilova, Cyril Furtlehner. AMStraMGRAM: Adaptive Multi-cutoff Strategy Modification for ANaGRAM. 2025. hal-05302122

HAL Id: hal-05302122

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

Preprint submitted on 13 Oct 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# AMSTRAMGRAM: ADAPTIVE MULTI-CUTOFF STRATEGY
## MODIFICATION FOR ANAGRAM

A PREPRINT

**Nilo Schwencke**

LISN-Université Paris-Saclay-INRIA-Saclay
nilo.schwencke@protonmail.com

**Cyriaque Rousselot**

LISN-Université Paris-Saclay-INRIA-Saclay
cyriaque.rousselot@inria.fr

**Alena Shilova**

LISN-Université Paris-Saclay-INRIA-Saclay
alena.shilova@inria.fr

**Cyril Furtlehner**

LISN-Université Paris-Saclay-INRIA-Saclay
cyril.furtlehner@inria.fr

October 13, 2025

## ABSTRACT

Recent works have shown that natural gradient methods can significantly outperform standard optimizers when training physics-informed neural networks (PINNs). In this paper, we analyze the training dynamics of PINNs optimized with ANaGRAM, a natural-gradient-inspired approach employing singular value decomposition with cutoff regularization. Building on this analysis, we propose a multi-cutoff adaptation strategy that further enhances ANaGRAM's performance. Experiments on benchmark PDEs validate the effectiveness of our method, which allows to reach machine precision on some experiments. To provide theoretical grounding, we develop a framework based on spectral theory that explains the necessity of regularization and extend previous shown connections with Green's functions theory.

**Keywords** Physics-Informed Neural Networks · Natural Gradient · Optimization · Partial Differential Equations · Neural Tangent Kernel

## 1 Introduction

Physics-informed neural networks (PINNs) have recently emerged as a promising alternative for the numerical solution of partial differential equations (PDEs) (Raissi et al., 2019). By leveraging neural networks as universal function approximators (Leshno et al., 1993), PINNs replace traditional mesh-based discretizations with sampling-based collocation methods, enabling a straightforward extension to high-dimensional domains. This mesh-free formulation not only circumvents the “curse of dimensionality” inherent in grid-based approaches, but also allows continuous evaluation of the solution throughout the domain without explicit mesh generation (Cuomo et al., 2022).

Despite these advantages, achieving low training error with PINNs remains a major challenge (Wang et al., 2023; Urbán et al., 2025; Kiyani et al., 2025; De Ryck et al., 2024). Open questions include how to select and distribute collocation points, how to balance the PDE residual against boundary-condition penalties, and which optimization strategies most effectively minimize the composite loss (Krishnapriyan et al., 2021; Wang et al., 2021; McClenny & Braga-Neto, 2022).

A different line of research has recently reexamined PINNs from the perspective of functional geometry (Müller & Zeinhofer, 2023, 2024; Jnini et al., 2024), providing a mathematically principled view of the training dynamics. In this vein, the ANaGRAM algorithm (Schwencke & Furtlehner, 2025) applies a natural-gradient update (Amari, 1998; Ollivier, 2015), based on a reinterpretation and generalization of the neural tangent kernel (NTK; Jacot et al. (2018)) as the kernel of the projection onto the neural network's tangent space. This leads to a notion of the empirical natural gradient that projects the true functional gradient onto the empirical tangent space, yielding significantly faster convergence and lower errors compared to standard optimizers on PDE benchmarks.
Nevertheless, while ANaGRAM improves over standard optimizers, it still falls short of the accuracy attained by classical mesh-based methods, such as the finite element method (Grossmann et al., 2024). Moreover, its final performance is highly affected by the way the pseudo-inverse of the feature matrix is computed. In particular, ANaGRAM sets a fixed level of cutoff: a value below which the singular values of the feature matrix are ignored, i.e. it controls how much loss signal is incorporated into an update. ANaGRAM’s cutoff is currently chosen manually, as no automatic selection procedure has been proposed.

In this paper, we study the performance and training dynamics of ANaGRAM, with a particular focus on the role of the chosen cutoff. Typically, the training loss of ANaGRAM exhibits the slow convergence at the early iterations followed by a sudden drop at the end of the training – similar behavior is shown by the eNGD method (Müller & Zeinhofer, 2023). We discover that it is closely connected to what we further refer as the flattening phenomenon, which we define and characterize using the reconstruction error: a novel metric that measures how much of the loss signal is lost by different choices of cutoffs. Relying on the adaptive multi-cutoff strategy, our new algorithm AMStraMGRAM manages to capitalize on this phenomenon, resulting in a significant improvement (of several orders of magnitude) on various PDE benchmarks. To complement our empirical findings, we also present a functional-analytic view linking cutoff (and ridge regularization) to (generalized) Green operator theory, clarifying why cutoff regularization is essential and not just a mere fix to stabilize training.

# 2 Problem Statement

## 2.1 Differential Operators and Physics-Informed Neural Networks (PINNs)

Let $\Omega \subset \mathbb{R}^d$ be a domain. We introduce two operators, $D$ and $B$, defined on a Hilbert space $\mathcal{H}$ of real-valued functions, acting respectively on $\Omega$ and on its boundary $\partial\Omega$:

$$D: \left\{ \begin{array}{rcl} \mathcal{H} & \to & L^2(\Omega \to \mathbb{R}, \mu) \\ u & \mapsto & D[u] \end{array} \right. , \qquad B: \left\{ \begin{array}{rcl} \mathcal{H} & \to & L^2(\partial\Omega \to \mathbb{R}, \sigma) \\ u & \mapsto & B[u] \end{array} \right. . \qquad (1)$$

Here, $D$ denotes a differential operator, while $B$ represents a boundary operator. A function $u \in \mathcal{H}$ is said to be a classical solution to the Partial Differential Equation (PDE) associated with $D$ and $B$ if it satisfies

$$\left\{ \begin{array}{ll} D(u) = f \in L^2(\Omega \to \mathbb{R}, \mu), & \text{in } \Omega, \\ B(u) = g \in L^2(\partial\Omega \to \mathbb{R}, \sigma), & \text{on } \partial\Omega, \end{array} \right. \qquad (2)$$

A physics-informed neural network (PINN) approximates the solution $u$ by a parametric model $u_\theta$, where $u_\theta$ is a neural network with parameters $\theta \in \mathbb{R}^p$. The learning objective is to minimize the empirical loss

$$\ell_{D,B}(\theta) := \frac{1}{2S_D} \sum_{i=1}^{S_D} (D[u_\theta](x_i^D) - f(x_i^D))^2 + \frac{1}{2S_B} \sum_{i=1}^{S_B} (B[u_\theta](x_i^B) - g(x_i^B))^2 . \qquad (3)$$

## 2.2 PINNs Optimizers

Training PINNs is notoriously challenging. Issues such as spectral bias, where networks struggle to learn high-frequency components, and the difficulty of balancing residual and boundary loss terms—often with vastly different magnitudes— result in unsatisfactory performance of standard deep learning optimizers (Wang et al., 2021; De Ryck et al., 2024; Krishnapriyan et al., 2021; Liu et al., 2024).

To mitigate these challenges, researchers have proposed various strategies. These include adaptive sampling approaches that focus on regions with high error (Krishnapriyan et al., 2021), dynamic loss weighting schemes (McClenny & Braga-Neto, 2022), and architectural modifications (Wang et al., 2024). Another promising line of research has focused on modifying the optimizers. In particular, two main branches of optimization approaches for PINNs have emerged:

(i) **Second-Order Methods.** These methods, based on Quasi-Newton techniques, particularly the BFGS algorithm (Nocedal & Wright, 1999, Chapter 6) and its memory-efficient approximation L-BFGS (Liu & Nocedal, 1989), address some of the training difficulties by considering the curvature of the loss landscape. This curvature arises from the non-linearities of both the neural network and the differential operators (Rathore et al., 2024). Recently, Urbán et al. (2025) extended this approach by modifying the self-scaled BFGS (SSBFGS; Al-Baali, 1998) and self-scaled Broyden (SSBroyden; Al-Baali & Khalfan, 2005), along with other computational enhancements such as point resampling (Wu et al., 2023) and boundary condition enforcement (Wang et al., 2023), achieving state-of-the-art results (Kiyani et al., 2025).
(ii) **Natural Gradient Methods.** In contrast to second-order methods, natural gradient methods are first-order techniques<fnref n="1" /> that provide a principled way to incorporate the geometry and metric structure of the problem space. Initially introduced in the context of information geometry by Amari (1998) and later extended by Ollivier (2015), these methods were introduced for PINNs by Müller & Zeinhofer (2023). In subsequent work, Schwencke & Furtlehner (2025) connected these methods to kernel methods, yielding an efficient implementation they linked to Green's function theory (Duffy, 2015).

## 2.3 Natural Gradient Methods for PINNs

As a preliminary observation highlighted in Schwencke & Furtlehner (2025, Section 4.1), PINNs can be interpreted as a quadratic regression problem. This viewpoint arises naturally once the parametric model $u_{\theta}$ is replaced with the following compound model:

$$ (D, B) \circ u: \left\{ \begin{array}{rcl} \mathbb{R}^P & \to & \mathcal{H} \\ \boldsymbol{\theta} & \mapsto & u_{\theta} \end{array} \right. \to \begin{array}{c} \mathrm{L}^2(\Omega, \mu) \times \mathrm{L}^2(\partial\Omega, \sigma) \\ (D[u_{\theta}], B[u_{\theta}]) \end{array}. \quad (4) $$

For ease of exposition, and without loss of generality, we restrict attention to regression in $\mathrm{L}^2(\Omega, \mu)$. Given $f \in \mathrm{L}^2(\Omega, \mu)$, we define the associated empirical loss

$$ \ell(\boldsymbol{\theta}) := \frac{1}{2S} \sum_{i=1}^{S} (u_{\theta}(x_i) - f(x_i))^2, \quad (5) $$

which can be seen as a discretization of the functional loss

$$ \mathcal{L}(u) := \frac{1}{2} \|u - f\|_{\mathrm{L}^2(\Omega, \mu)}^2, \quad u \in \mathrm{L}^2(\Omega, \mu). \quad (6) $$

The natural gradient approach seeks to compute the optimal update direction in function space and then pull it back to parameter space. A single Fréchet derivative of the functional loss Equation (6) yields $\nabla \mathcal{L}_{|u} = u - f$. The key insight is that admissible updates are constrained to the tangent space of the parametric model,

$$ T_{\theta} \mathcal{M} := \mathrm{Im}(du_{\theta}) = \mathrm{Span} (\partial_{p} u_{\theta} : 1 \le p \le P) \subset \mathcal{H}, \quad (7) $$

where $\mathcal{M} := \mathrm{Im}(u) = \{u_{\theta} : \theta \in \mathbb{R}^P\} \subset \mathcal{H}$ is the manifold of functions parametrized by $\theta$. Thus, the optimal update in function space is the projection of $\nabla \mathcal{L}_{|u}$ onto the tangent space (cf. Figure 4),

$$ u_{\theta_{t+1}} \leftarrow u_{\theta_t} - \eta \Pi_{T_{\theta_t} \mathcal{M}} (\nabla \mathcal{L}_{u_{\theta_t}}); \quad \theta_{t+1} \leftarrow \theta_t - \eta \mathrm{d}u_{\theta_t}^{\dagger} (\Pi_{T_{\theta_t} \mathcal{M}} (\nabla \mathcal{L}_{u_{\theta_t}})), \quad (8) $$

where the second equation is simply the pullback of the functional update to parameter space. We prove in Section H.1 that this update is equivalent to the Gram–matrix formulation:

$$ \theta_{t+1} \leftarrow \theta_t - \eta G_{\theta_t}^{\dagger} \nabla \ell(\theta_t); \quad G_{\theta_t, p, q} := \langle \partial_p u_{\theta_t}, \partial_q u_{\theta_t} \rangle_{\mathrm{L}^2(\Omega, \mu)}. \quad (9) $$

## 2.4 ANaGRAM: Empirical Natural Gradient

The $O(P^3)$ complexity of matrix inversion in Equation (9) renders a direct implementation prohibitively expensive. ANaGRAM (Schwencke & Furtlehner, 2025) circumvents this by exploiting a motivated approximation. The key observation is that the update can be expressed in terms of the empirical feature matrix $\hat{\phi} \in \mathbb{R}^{S \times P}$ and the empirical functional residuals $\hat{\mathcal{L}}_{\theta} \in \mathbb{R}^S$:

$$ \theta_{t+1} \leftarrow \theta_t - \eta \hat{\phi}^{\dagger} \nabla \mathcal{L}_{\theta_t}; \quad \hat{\phi}_{i,p} := \partial_p u_{\theta}(x_i); \quad (\nabla \mathcal{L}_{\theta})_i := u_{\theta}(x_i) - f(x_i). \quad (10) $$

Here, the pseudo-inverse is computed via singular value decomposition (SVD): $\hat{\phi}^{\dagger} = \hat{U}\hat{\Delta}^{\dagger}\hat{V}^T$ with $\hat{\phi} = \hat{V}\hat{\Delta}\hat{U}^T$, where $\hat{U} \in \mathbb{R}^{P \times r_{svd}}$, $\hat{\Delta} \in \mathbb{R}^{r_{svd} \times r_{svd}}$, $\hat{V} \in \mathbb{R}^{S \times r_{svd}}$, and $r_{svd} = \min(P, S)$. This reduces computational cost to $O(\min(PS^2, P^2S))$, which is tractable in practice. A comparable complexity was later obtained by Guzmán-Cordero et al. (2025) using a Cholesky factorization approach.

For further details on the derivation of the empirical natural gradient, we refer to Schwencke & Furtlehner (2025). In what follows, we adopt a slight abuse of notation by omitting the explicit dependence on $\theta$ whenever it is clear from context. When iteration indices matter, we explicitly write $t$ to emphasize the connection to $\theta_t$.

<footnote n="1">¹contrary to a widespread misconception, which arises from their analogy in the context of information theory</footnote>
## 2.5 Regularization

As discussed in Section G.1, the type of problem we consider is ill-conditioned, which necessitates the use of regularization. We distinguish between two main regularization schemes: (i) ridge regression, which consists in adding a factor $\alpha^2 I_d$ (or, according to conventions, $\alpha^{-2} I_d$) to the Gram matrix $G_\theta$ in Equation (9) (or its approximation $\hat{G}_\theta$), thereby making it invertible or (ii) cutoff regularization, a scheme that applies a binary threshold (used in ANaGRAM):

$$ \hat{\Delta}_{t,i} = \begin{cases} \hat{\Delta}_{t,i}^{-1}, & \text{if } \hat{\Delta}_{t,i} \ge \alpha, \\ 0, & \text{otherwise.} \end{cases} \qquad (11) $$

Here $\alpha$ denotes the cutoff threshold. This regularization is the focus of our analysis in Section 3. For completeness, we provide a geometric interpretation of each scheme in Section G. We further show that cutoff regularization extends previously established connections between natural gradient methods and Green's function theory (Schwencke & Furtlehner, 2025). In particular, we obtain:

**Theorem 1.** The generalized Green's function of the operator $D$ in the regularized space $H_{D, \mathcal{H}_0}^\alpha$ is given, for all $x, y \in \Omega$, by

$$ g_D(x, y) := D[k_D(x, \cdot)](y), \qquad (12) $$

where $H_{D, \mathcal{H}_0}^\alpha$ is a regularized space with reproducing kernel $k_D$, defined in Section G.4.

# 3 Insights on ANaGRAM's Training Dynamics

In this section, we will look at relevant quantities of interest to understand this empirical phenomenon.

## 3.1 Reconstruction Error of Functional Gradient

Let $\theta \in \mathbb{R}^P$, the empirical feature matrix $\hat{\phi} \in \mathbb{R}^{S \times P}$, and the empirical functional gradient $\widehat{\nabla L} \in \mathbb{R}^S$ as defined in Equation (10). Let us consider various empirical tangent spaces formed by taking different ranges of right singular vectors of $\hat{\phi} = \hat{U}\hat{\Delta}\hat{V}^T$, i.e. $\hat{T}_N^M \mathcal{M} = \text{Span}(\hat{V}_{t,i} : M \le i \le N)$. For $1 \le N \le r_{svd}$, reconstruction error measures how much information from the functional gradient signal is lost when considering only first $N$ components in SVD (the error caused by the projection onto the empirical tangent space $\hat{T}_N^0 \mathcal{M}$) is defined as follows

$$ \text{RCE}_N^S := \frac{1}{\sqrt{S}} \| \hat{V} \Pi_N^0 \hat{V}^T \widehat{\nabla L} - \widehat{\nabla L} \|_{\mathbb{R}^S} = \frac{1}{\sqrt{S}} \| \Pi_{\hat{T}_N^0 \mathcal{M}}^{\perp} \widehat{\nabla L} - \widehat{\nabla L} \|, \qquad (13) $$

where we define $\Pi_N^M \in \mathbb{R}^{r_{svd} \times r_{svd}}$ as a projection operator onto $\hat{T}_N^M \mathcal{M}$:

$$ \Pi_N^M = \sum_{p=M+1}^{N} e^{(p)} e^{(p)^T}, \qquad (14) $$

with $(e^{(p)})_{1 \le p \le r_{svd}}$ being the canonical basis of $\mathbb{R}^{r_{svd}}$.

**Proposition 1.** $RCE_N^S$ is a non-increasing function of $N$, i.e. for all $1 \le M, N \le r_{svd}$:

$$ M \le N \implies RCE_M^S \ge RCE_N^S. \qquad (15) $$

Furthermore, assuming that $(x_i)_{i=1}^S$ are i.i.d sampled from $\mu$, we have $\mu$-almost surely

$$ \lim_{S \to \infty} RCE_N^S = \| \nabla \mathcal{L}_{u_\theta} - \Pi_{\hat{T}_N^0 \mathcal{M}}^{\perp} \nabla \mathcal{L}_{u_\theta} \|_{L^2(\Omega, \mu)} = \| \Pi_{[\hat{T}_N^0 \mathcal{M}]^{\perp}}^{\perp} \nabla \mathcal{L}_{u_\theta} \|_{L^2(\Omega, \mu)}, \qquad (16) $$

where $\hat{T}_N^M \mathcal{M} = \text{Span}(V_{t,i} : M \le i \le N)$, while $(V_{t,i})_{1 \le i \le r_{svd}}$ are the right singular-vectors of the differential $du_\theta$ ordered in a decreasing order according to their associated singular values.

**Remark 1.** Note that $\hat{V}_{t,i} \in \mathbb{R}^S$ for $i \in 1, \dots, N$, the right singular vectors of $\hat{\phi}$, can be seen as discretized versions of $V_{t,i}$ from Proposition 1. Indeed, a weak convergence holds, i.e. $\forall h \in \mathcal{H}$, $\frac{1}{S} \sum_{j=1}^S \hat{V}_{t,i,j} h_j = \frac{1}{S} \sum_{j=1}^S V_{t,i}(x_j) h(x_j) \xrightarrow{S \to \mathcal{C}} \langle V_{t,i}, h \rangle_{L^2}$.

Proof of Proposition 1 can be found in Section H.3. From Proposition 1 RCE is related to the concept of expressivity bottleneck illustrated in Verbockhaven et al. (2024), and measures what part of the learning signal is not captured by truncating at $N$ components for natural gradient computation. Therefore, this metric allows us to explicitly estimate and compare different cutoff choices. Note that this metric incurs no additional computational cost since ANaGRAM already computes the required SVD.
## 3.2 Empirical Observations: Flattening

Here we illustrate the evolution of training loss and reconstruction error, where Figure 1 schematically outlines key stages of ANaGRAM's training dynamics. The plot of a real experiment is provided in Section E.

<figure><img src="image_3.png" /><figcaption>Figure 1: ANaGRAM training dynamics. Legend (top) and four key phases: (a) initial evolution, (b) reconstruction-singular value intersection passes target precision, (c) emergence of the flattening regime, (d) complete flattening yielding final loss level. Despite changing scale, target precision is constant and fixed across all plots. The number of ANaGRAM's retained components $r_{cutoff}$ is at intersection of precision line with singular values curve.</figcaption></figure>

Let $\alpha$ is a cutoff level (also referred to as precision) and $r_{cutoff}$ denote the number of components retained by the cutoff, i.e., $r_{cutoff}(t) = \max\{j : \hat{\Delta}_{t,j} \ge \alpha\}$. In Figure 1, we observe different stages of the training. First, the reconstruction error is above the wanted precision (Figure 5a). As the training progresses, the training loss drops and the reconstruction error drops until reaching the cutoff precision (Figure 5b). Eventually, the reconstruction error drops below the cutoff threshold (Figure 5c). During this phase, the training loss (corresponding to the RCE for 0 component (green line in the figure)) is not decreasing a lot.

Then, a phenomenon that we call "flattening" occurs: once the reconstruction error is small compared to the cutoff precision value, reconstruction error *flattens* over the interval $[N_{flat}, r_{cutoff}]$, where $N_{flat}$ is the smallest number such as

$$RCE_{N_{flat}}^{S} - RCE_{r_{cutoff}}^{S} \approx 0. \qquad (17)$$
Eventually, the phenomenon propagates toward low numbers of retained components (Figure 5e) and $N_{flat} = 0$. Reconstruction error is now constant for all retained components and the training ends with training loss at cutoff precision. We refer a reader to Section H.3 to have a more theoretical insight on what is happening during the flattening. Remark 2. This phenomenon sheds light on the sharp drop in training loss observed near the end of optimization, as reported in Schwencke & Furtlehner (2025). By combining Equations (5), (10) and (13) and using that $\Pi_0^0 = 0$, we obtain

$$RCE_0^{S^2} \stackrel{13}{=} \frac{1}{S} \left\| \hat{V} \Pi_0^0 \hat{V}^T \hat{\nabla} \mathcal{L} - \hat{\nabla} \mathcal{L} \right\|_{\mathbb{R}^S}^2 = \frac{1}{S} \left\| \hat{\nabla} \mathcal{L} \right\|_{\mathbb{R}^S}^2 \stackrel{10}{=} \frac{1}{S} \sum_{i=1}^S (u_{\theta}(x_i) - f(x_i))^2 \stackrel{5}{=} \ell(\theta). \quad (18)$$

Thus, the last iteration of flattening is **directly responsible for the sudden drop of train loss** at the end of the training. Remark 3. We see that for higher precision than the cutoff value ($N > r_{cutoff}$), the RCE is still decreasing as we increase the number of components kept. This indicates that there is still information to capture in the functional eigenspace composed of components associated to lower eigenvalues, see also Section H.3.

The final interesting observation is that

$$RCE_0^S - RCE_{r_{cutoff}}^S \approx 0 \quad \Leftrightarrow \quad \Pi_{r_{cutoff}}^0 \hat{V}^T \hat{\nabla} \mathcal{L} \approx 0. \quad (19)$$

Thus, the flattening phenomenon means that the projection of the signal onto the first $r_{cutoff}$ components retained by the cutoff is negligible. In other words, the optimization has extracted all the *usable* signal from these components at this cutoff level.

## 3.3 Incomplete Flattening and Adaptive Strategies

In practice, for some experiments we observe that the flattening may remain incomplete with $\lim_{t \to \infty} N_{flat} = N_{flat}^\infty > 0$: the system remains in a state similar to that shown in Figure 1c and never (at least not within a reasonable number of iterations) reaches the configuration illustrated in Figure 1d. A natural question arises: *what happens if we adjust the cutoff to retain exactly $N_{flat}^\infty$ components?*

If we try this trick in practice (see Figure 6), then a single natural gradient step with an adjusted cutoff can be enough to get immediate and complete flattening ($N_{flat} = 0$) and eventually dramatically reduce training loss. This abrupt flattening when restricting cutoff to low number of feature is typically accompanied by a learning rate found by the line search to be very close to one. A possible explanation is that this may represent an iteration in the *lazy training* regime (NTK and the feature matrix are nearly constant), where we regress linearly (and thus fast) based on learned features. This hypothesis should be further explored in future work.

This empirical insight motivates the use of an adaptive algorithm: by dynamically adjusting cutoffs, we can hope to accelerate convergence and achieve higher precision.

# 4 Algorithmic Design: Exploiting Flattening

Building upon the empirical analysis presented in Section 3, we develop a principled algorithm that controls and exploits the flattening phenomenon identified in ANaGRAM's training dynamics. Our approach is based on tracking the relationship between reconstruction error and singular values to automatically determine well-adapted cutoff in order to reach the target precision (error) $\epsilon$ at the end of the training. This well-adapted cutoff should vary from one iteration to another to adjust to the currently learned weights and training dynamics in such a way to avoid early flattening (if flattening happens too early, the training stagnates at higher values of losses) and when intersection between RCE and singular values goes below the target precision $\epsilon$, we enforce the flattening, so that the final training loss also drops to $\epsilon$.

## 4.1 Adaptive Cutoff Strategy

In what follows, we suggest an adaptive cutoff rank $r_{cutoff}$ that indicates how much components of $\hat{\Delta}$ are retained for the next update of ANaGRAM. Our algorithm operates by dynamically selecting cutoff ranks based on the relationship between reconstruction error and singular values:

$$r_{cutoff}(t) = \begin{cases} r_{\text{int}}(t) := \max\{j: \text{RCE}_j^S(t) \le \hat{\Delta}_{t,j}\} & \text{if } \text{RCE}_{r_{\text{int}}(t)}^S(t) > \epsilon \text{ (intersection rank),} \\ r_{\epsilon}(t) := \max\{j: \text{RCE}_j^S(t) \ge \epsilon\} & \text{if } \text{RCE}_{r_{\text{int}}(t)}^S(t) \le \epsilon \text{ (precision rank).} \end{cases} \quad (20)$$

The algorithm terminates when $r_{\epsilon}(t) = 0$, indicating that the reconstruction error $RCE_0^S$ that is equal to the training error is indeed below the predefined precision threshold.
For ease of presentation, we provide only the core elements of AMStraMGRAM in Algorithm 1 consisting in adaptively choosing, which $r_{cutoff}$ to apply for $\hat{\Delta}$ at each update of ANaGRAM. The final algorithm is explained in Section C. Final Algorithm 2 addresses some irregularities observed in evolution of RCE and singular values that we explain in more details in Section C.4.

<figure><img src="image_4.png" /><figcaption>Figure 2: **Dynamics of the adaptive multi-cutoff strategy in AMStraMGRAM.** Progression from (a) initial exploration, (b) intersection reaches precision, (c) flattening onset, to (d) converged state. Red arrows (when present) indicate the retained rank dynamics (pointing right – increasing, pointing left – decreasing). Legends are shown below.</figcaption></figure>

## 4.2 Geometrical Interpretation of the Adaptive Strategy

The algorithm exploits the geometric relationship between the empirical tangent space and the functional gradient. By tracking the intersection, we maximize the projection of the functional gradient onto the empirical tangent space while staying out of flattening. Once the intersection reaches the precision level, we exploit the flattening phenomenon to achieve prescribed precision.

According to Proposition 1, the reconstruction error $\text{RCE}_N^S$ measures how much of the functional gradient signal remains to be captured by the first $N$ components. The intersection point thus represents the good balance between signal capture and phase transition.
# 5 Experimental Results

We first compare in Table 1 our method implemented in JAX<fnref n="2" /> with the ANaGRAM method (Schwencke & Furtlehner, 2025) on the benchmark problems presented in their paper, with modified datasets. As we see, for every equation, we perform better.

<table><caption>Table 1: Performance comparison between AMStraMGRAM (our method) and ANaGRAM Schwencke & Furtlehner (2025). The adaptive strategy demonstrates significant improvements across all benchmark problems, with L2 error improvements of up to 8 orders of magnitude.</caption><thead><tr><th rowspan="2">Experiment</th><th colspan="2">Mean Squared Error (MSE)</th><th colspan="2">L<sub>2</sub> Error</th></tr><tr><th>Ours</th><th>ANaGRAM</th><th>Ours</th><th>ANaGRAM</th></tr></thead><tbody><tr><td>Heat Equation</td><td><b>6.29e-29</b> ± <b>6.78e-30</b></td><td>8.56e-11 ± 7.05e-11</td><td><b>2.32e-14</b> ± <b>1.14e-14</b></td><td>1.28e-06 ± 1.75e-06</td></tr><tr><td>Laplace 2D</td><td><b>1.46e-28</b> ± <b>1.87e-29</b></td><td>4.27e-13 ± 4.66e-13</td><td><b>2.24e-15</b> ± <b>2.52e-16</b></td><td>3.49e-09 ± 3.58e-09</td></tr><tr><td>Laplace 5D</td><td><b>2.04e-08</b> ± <b>1.16e-08</b></td><td>6.37e-08 ± 7.01e-08</td><td><b>2.12e-05</b> ± <b>8.15e-06</b></td><td>4.00e-05 ± 2.93e-05</td></tr><tr><td>Allen-Cahn</td><td><b>3.19e-11</b> ± <b>2.37e-11</b></td><td>2.19e-04 ± 4.16e-04</td><td><b>5.87e-05</b> ± <b>6.25e-06</b></td><td>4.32e-03 ± 5.93e-03</td></tr></tbody></table>

We then compare our method with the baseline methods from Urbán et al. (2025) on the benchmark problems presented in their paper. Note that in our case we do not need to enforce boundary constraints. The methodology of sampling is also slightly different, as we sample the data from a fixed grid, following the methodology of Schwencke & Furtlehner (2025), while in Urbán et al. (2025) they perform batching of randomly sampled points.

<table><caption>Table 2: Performance comparison between AMStraMGRAM (our method) and baseline Urbán et al. (2025) methods. Our method demonstrates improvements across benchmark problems, without requiring enforcement of boundary constraints.</caption><thead><tr><th rowspan="2">Experiment</th><th colspan="2">Mean Squared Error (MSE)</th><th colspan="2">L<sub>2</sub> Error</th></tr><tr><th>Ours</th><th>SSBroyden*</th><th>Ours</th><th>SSBroyden*</th></tr></thead><tbody><tr><td>One-dimensional Burgers (1DB)</td><td><b>2.99e-12</b> ± <b>9.26e-13</b></td><td>2.92e-10 ± 1.45e-10</td><td><b>1.5e-06</b> ± <b>9.43e-7</b></td><td>1.59e-06 ± 1.02e-6</td></tr><tr><td>Non-Linear Poisson (k=1)</td><td><b>8.51e-24</b> ± <b>2.24e-24</b></td><td>3.03e-16 ± 3.82e-16</td><td>6.81e-10 ± 1.41e-09</td><td><b>9.29e-12</b> ± <b>5.85e-12</b></td></tr><tr><td>Allen-Cahn (AC)</td><td>3.19e-11 ± 2.37e-11</td><td><b>6.42e-12</b> ± <b>5.52e-12</b></td><td>5.87e-05 ± 6.25e-06</td><td><b>3.94e-06</b> ± <b>1.72e-06</b></td></tr></tbody></table>

# 6 Limitations

Despite its effectiveness, AMStraMGRAM can exhibit overfitting, particularly in problems with sharp features like the Allen-Cahn equation. The algorithm drives the training error to machine precision on the sampled points, but the learned function may develop high-frequency oscillations between them, especially in regions of high curvature where the approximation is the most challenging. These artifacts, visible as “overfitting lines” in Figure 3, are an imprint of the sampling lattice (see regions around $x = \pm 0.5$). They arise because the SVD cutoff effectively projects the update onto a low-rank subspace of the tangent space. This subspace is often aligned with the grid axes, leading to anisotropic smoothing that perfectly fits the data on the grid lines but interpolates poorly in the under-sampled regions between them. Once the flattening phase begins, the training enters a quasi-linear regime that can “lock in” these geometric artifacts.

This phenomenon highlights that while our method significantly improves on ANaGRAM, the quality of the final solution remains fundamentally limited by the sampling strategy. Mitigating such overfitting requires co-designing the sampler and the optimizer. Potential remedies include adaptive sampling, where new collocation points are added in regions of high reconstruction error, or curriculum-based approaches that progressively refine the sampling grid.

<footnote>* refer to method from Urbán et al. (2025) with adaptive sampling and hard constraint enforcement on boundary conditions.</footnote>

<footnote n="2">²https://anonymous.4open.science/r/AMStraMGRAM-8D1B/</footnote>
<figure><img src="image_5.png" /><figcaption>Figure 3: Allen-Cahn overfitting: residual lines align with sampling lines. Low-rank (post-cutoff) tangent projections fit exactly on sampled fibers while interpolation between them inherits weakly constrained oscillations in regions of steep interface curvature.</figcaption></figure>

# 7 Conclusion

In this work, we have introduced AMStraMGRAM, an adaptive multi-cutoff strategy that enhances the ANaGRAM natural gradient method for training PINNs. Our work provides an analytical framework to explain ANaGRAM's convergence behavior, uncovering a *flattening* phenomenon that clarifies its training dynamics. The proposed algorithm automatically adjusts cutoff regularization. Notably, AMStraMGRAM exhibits “overfitting” as demonstrated in Allen-Cahn experiments. These results underscore the potential of natural gradient optimization for PINNs while highlighting the critical role of sampling strategies in realizing their full accuracy.

Future research will focus on integrating residual-based methods to further stabilize training, establishing rigorous convergence guarantees for our adaptive cutoff scheme, and extending the approach to higher-dimensional PDEs and complex geometries. Exploring the interplay between network architecture and optimization—as well as further developing sampling techniques—will be essential to address the fundamental challenge of balancing optimization power with data representation. Ultimately, our findings suggest that with careful algorithmic design, PINNs can achieve the precision required for practical scientific computing, paving the way for mesh-free methods in computational science.
References

Ben Adcock and Daan Huybrechs. Frames and Numerical Approximation. *SIAM Review*, 61(3):443-473, January 2019. ISSN 0036-1445, 1095-7200. doi:10.1137/17M1114697.

Ben Adcock and Daan Huybrechs. Frames and numerical approximation II: Generalized sampling, July 2020.

M. Al-Baali. Numerical Experience with a Class of Self-Scaling Quasi-Newton Algorithms. *Journal of Optimization Theory and Applications*, 96(3):533-553, March 1998. ISSN 0022-3239, 1573-2878. doi:10.1023/A:1022608410710.

Mehiddin Al-Baali and Humaid Khalfan. Wide interval for efficient self-scaling quasi-Newton algorithms. *Optimization Methods and Software*, 20(6):679-691, December 2005. ISSN 1055-6788, 1029-4937. doi:10.1080/10556780410001709448.

Shun-Ichi Amari. Natural gradient works efficiently in learning. *Neural computation*, 10(2):251-276, 1998.

Yurij M. Berezansky, Zinovij G. Sheftel, and Georgij F. Us. *Functional Analysis. Vol. II*, volume 86 of *Operator Theory Advances and Applications*. Birkhäuser, 1996.

Salvatore Cuomo, Vincenzo Schiano Di Cola, Fabio Giampaolo, Gianluigi Rozza, Maziar Raissi, and Francesco Piccialli. Scientific Machine Learning Through Physics-Informed Neural Networks: Where we are and What's Next. *Journal of Scientific Computing*, 92(3):88, July 2022. ISSN 1573-7691. doi:10.1007/s10915-022-01939-z.

Tim De Ryck, Florent Bonnet, Siddhartha Mishra, and Emmanuel de Bézenac. An operator preconditioning perspective on training in physics-informed machine learning, May 2024.

Dean G. Duffy. *Green's Functions with Applications*. Chapman and Hall/CRC, 2015.

Tamara G Grossmann, Urszula Julia Komorowska, Jonas Latz, and Carola-Bibiane Schönlieb. Can physics-informed neural networks beat the finite element method? *IMA Journal of Applied Mathematics*, 89(1):143-174, January 2024. ISSN 0272-4960. doi:10.1093/imamat/hxaе011.

Andrés Guzmán-Cordero, Felix Dangel, Gil Goldshlager, and Marius Zeinhofer. Improving Energy Natural Gradient Descent through Woodbury, Momentum, and Randomization, May 2025.

Arthur Jacot, Franck Gabriel, and Clément Hongler. Neural tangent kernel: Convergence and generalization in neural networks. *Advances in neural information processing systems*, 31, 2018.

Anas Jnini, Flavio Vella, and Marius Zeinhofer. Gauss-Newton Natural Gradient Descent for Physics-Informed Computational Fluid Dynamics, February 2024.

Andreas Kirsch. *An Introduction to the Mathematical Theory of Inverse Problems*, volume 120 of *Applied Mathematical Sciences*. Springer International Publishing, Cham, 2021. ISBN 978-3-030-63342-4 978-3-030-63343-1. doi:10.1007/978-3-030-63343-1.

Elham Kiyani, Khemraj Shukla, Jorge F. Urbán, Jérôme Darbon, and George Em Karniadakis. Which Optimizer Works Best for Physics-Informed Neural Networks and Kolmogorov-Arnold Networks?, April 2025.

Rainer Kress. *Linear Integral Equations*, volume 82 of *Applied Mathematical Sciences*. Springer, New York, NY, 2014. ISBN 978-1-4614-9592-5 978-1-4614-9593-2. doi:10.1007/978-1-4614-9593-2.

Aditi Krishnapriyan, Amir Gholami, Shandian Zhe, Robert Kirby, and Michael W Mahoney. Characterizing possible failure modes in physics-informed neural networks. In *Advances in Neural Information Processing Systems*, volume 34, pp. 26548-26560. Curran Associates, Inc., 2021.

Moshe Leshno, Vladimir Ya Lin, Allan Pinkus, and Shimon Schocken. Multilayer feedforward networks with a nonpolynomial activation function can approximate any function. *Neural networks*, 6(6):861-867, 1993.

Dong C. Liu and Jorge Nocedal. On the limited memory BFGS method for large scale optimization. *Mathematical Programming*, 45(1):503-528, August 1989. ISSN 1436-4646. doi:10.1007/BF01589116.

Songming Liu, Chang Su, Jiachen Yao, Zhongkai Hao, Hang Su, Youjia Wu, and Jun Zhu. Preconditioning for Physics-Informed Neural Networks, February 2024.

Levi McClenny and Ulisses Braga-Neto. Self-Adaptive Physics-Informed Neural Networks using a Soft Attention Mechanism, April 2022.

Johannes Müller and Marius Zeinhofer. Achieving high accuracy with PINNs via energy natural gradient descent. In *International Conference on Machine Learning*, pp. 25471-25485. PMLR, 2023.

Johannes Müller and Marius Zeinhofer. Position: Optimization in SciML Should Employ the Function Space Geometry. In *Forty-First International Conference on Machine Learning*, February 2024.

Jorge Nocedal and Stephen J. Wright. *Numerical Optimization*. Springer, 1999.
Yann Ollivier. Riemannian metrics for neural networks I: Feedforward networks. *Information and Inference: A Journal of the IMA*, 4(2):108–153, 2015.

Vern I. Paulsen and Mrinal Raghupathi. *An Introduction to the Theory of Reproducing Kernel Hilbert Spaces*, volume 152. Cambridge university press, 2016.

M. Raissi, P. Perdikaris, and G.E. Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. *Journal of Computational Physics*, 378:686–707, February 2019. ISSN 00219991. doi:10.1016/j.jcp.2018.10.045.

Pratik Rathore, Weimu Lei, Zachary Frangella, Lu Lu, and Madeleine Udell. Challenges in Training PINNs: A Loss Landscape Perspective, February 2024.

Nilo Schwencke and Cyril Furtlehner. ANaGRAM: A natural gradient relative to adapted model for efficient PINNs learning. In *The Thirteenth International Conference on Learning Representations*, 2025.

Jorge F. Urbán, Petros Stefanou, and José A. Pons. Unveiling the optimization process of physics informed neural networks: How accurate and competitive can PINNs be? *Journal of Computational Physics*, 523:113656, February 2025. ISSN 0021-9991. doi:10.1016/j.jcp.2024.113656.

Manon Verbockhaven, Sylvain Chevallier, and Guillaume Charpiat. Growing tiny networks: Spotting expressivity bottlenecks and fixing them optimally. *arXiv preprint arXiv:2405.19816*, 2024.

Sifan Wang, Yujun Teng, and Paris Perdikaris. Understanding and mitigating gradient flow pathologies in physics-informed neural networks. *SIAM Journal on Scientific Computing*, 43(5):A3055–A3081, 2021.

Sifan Wang, Shyam Sankaran, Hanwen Wang, and Paris Perdikaris. An Expert's Guide to Training Physics-informed Neural Networks, August 2023.

Sifan Wang, Bowen Li, Yuhan Chen, and Paris Perdikaris. PirateNets: Physics-informed Deep Learning with Residual Adaptive Networks. *Journal of Machine Learning Research*, 25(402):1–51, 2024. ISSN 1533-7928.

Chenxi Wu, Min Zhu, Qinyang Tan, Yadhu Kartha, and Lu Lu. A comprehensive study of non-adaptive and residual-based adaptive sampling for physics-informed neural networks. *Computer Methods in Applied Mechanics and Engineering*, 403:115671, 2023.
## A Illustration of Natural Gradient

<figure><img src="image_6.png" /><figcaption>Figure 4: Illustration of the orthogonal projection of the functional gradient onto the tangent space. While the ideal update direction would be the functional gradient $\nabla \mathcal{L}_{|u_{\theta}|}$ (shown in blue), our model constrains us to follow directions within the tangent space $T_{\theta}M$ (shown as a green plane). The optimal feasible direction is thus the orthogonal projection $\Pi_{T_{\theta}M} (\nabla \mathcal{L}_{|u_{\theta}|})$ (shown in red).</figcaption></figure>

## B Our vocabulary

* **Domain** ($\Omega$).

* **Boundary** ($\partial\Omega$).

* **Differential operators** ($D, B$).

* **Cutoff** ($\alpha_t$). A threshold below which the components of the matrix $\hat{\Delta}$ are truncated, i.e. $\hat{\Delta} \leftarrow \begin{cases} \hat{\Delta} & \text{if } \hat{\Delta} \ge \alpha_t, \\ 0 & \text{else.} \end{cases}$

* **Full rank** ($r_{svd}$). A full rank of feature matrix $\hat{\phi}$ that we assume, without loss of generality, to be equal to $\min(P, S)$.

* **Rank** ($r_{cutoff}$). A number of $\hat{\Delta}$ components that are retained when computing a pseudo-inverse of $\hat{\Delta}$ in ANaGRAM. Depending on a current regime of the training and a desired effect, it can be set at $r_{int}$ or $r_e$.

* **Flattening**. The phenomenon described in Section 3.2, when reconstruction error starts to stabilize for a range of possible ranks.

* **Flat cutoff** ($N_{flat}$). A number of components that corresponds to the beginning of flattening in reconstruction error curve.

* **Feature matrix** ($\hat{\phi} \in \mathbb{R}^{P \times S}$). It is defined by a jacobian $\partial_p u_{\theta}(x_i)$, which is used in an ANaGRAM's update to "project" a functional gradient onto parameter space of $\theta$.
* **Precision** (ε). A hyperparameter of AMStraMGRAM that prescribes a target error level that the algorithm should achieve.

* **Intersection rank** ($r_{\text{int}}$). Defined in Equation (20), roughly speaking it corresponds to a number of components at which reconstruction error and singular values curves are intersecting.

* **Precision rank** ($r_{\epsilon}$). Defined in Equation (20), it corresponds to a number of components at which reconstruction error curve and precision level are intersecting.

* **Functional gradient** ($\nabla \mathcal{L}$). A Frechet derivative of squared $L^2$ loss $\mathcal{L}$, its negative gives the "ideal" update direction in non-parametric case.

* **Empirical functional gradient** ($\nabla \mathcal{L} \in \mathbb{R}^S$). A vector obtained by evaluating $\nabla \mathcal{L}$ on some finite number of samples $x_i \in \Omega$, for $i \in 1, \dots, S$.

* **Parametric model** ($u_{\theta}$). A function parametrized with $\theta$ that serves to approximate a solution to a problem (regression or PDE). Typically, it is a neural network, where $\theta$ are its full set of weights.

* **Differential of the model** ($du_{\theta}$). Defined as $du_{\theta}(h) = \sum_{p=1}^{P} h_p \frac{\partial u}{\partial \theta_p} = \lim_{\epsilon \to 0} \frac{u_{\theta+\epsilon h} - u_{\theta}}{\epsilon}$. It measures how much $u_{\theta}$ changes in a given direction $h$.

* **Tangent space** ($T_{\theta}M$). Image of a differential of the model, giving a space of possible updates for a model $u_{\theta}$.

* **SVD components of** $\hat{\phi}(\hat{U}, \hat{\Delta}, \hat{V})$. In particular, $\hat{\phi} = \hat{U}\hat{\Delta}\hat{V}^T$, where $\hat{U} \in \mathbb{R}^{P \times S}$ is a left singular vector matrix, $\hat{\Delta} \in \mathbb{R}^{r_{\text{svd}} \times r_{\text{svd}}}$ is a diagonal matrix with singular values on a diagonal ordered in a decreasing order and $\hat{V}$ is a right singular vector matrix.

* **Functional singular vectors** ($V_{t,i}$). Right singular vectors of the differential $du_{\theta}$.

* **Empirical tangent space** ($T_{\theta}^M M$). A subspace of tangent space $T_{\theta}M$, restricted to a span of the right functional singular vectors $V_{t,i}$ corresponding to a range of components from $M$ to $N$, i.e. $\text{Span}(V_{t,i}: 1 \le M \le N \le N)$.

* **Discretized empirical tangent space** ($\tilde{T}_{\theta}^M M$). A version of $T_{\theta}^M M$ discretized on a set of samples $\{x_i\}_{i=1}^S$ coming from $\Omega$.

* **Reconstruction error** ($RCE_N^S$). A measure identifying the portion of the functional gradient signal that is lost when restricting $\nabla \mathcal{L}$ to $\tilde{T}_{\theta}^0 M$.

* **Feature development phase**. The early phase in the training, during which high volatility is observed in both quantities of interest with high sensitivity to the choice of $r_{\text{cutoff}}$.

* **Flattening phase**. The later phase in the training, during which reconstruction error starts to flatten for some values of $N$, at the same time singular values dominate over reconstruction error for all retained components, resulting in a drop of training loss.

# C Practical Implementation Considerations

While the principled algorithm discussed in the main paper and summarized in Algorithm 1 provides a sound framework, empirical observations reveal that additional mechanisms are necessary for robust performance across diverse PDE problems. This section describes additional modifications to make the algorithm more practical.

## C.1 The Dual Cutoff Strategy: Addressing Empirical Challenges

Our experiments reveal that the single cutoff approach, while theoretically elegant, suffers from numerical instabilities and incomplete convergence in practice. We observed three critical issues:

1. **Ignition failure:** The intersection between reconstruction error and singular values sometimes fails to evolve, preventing the algorithm from reaching lower error values.

2. **Retreating dynamics:** The intersection rank may decrease during training, disrupting convergence.

3. **Incomplete flattening:** Without additional stabilization, the flattening phenomenon may not complete, leading to suboptimal final accuracy.

To address these challenges, we introduce a dual cutoff strategy inspired by the staged design of rocket launches:
## C.2 Three-Phase Training Dynamics

### C.2.1 Ignition Phase

We initialize two cutoffs:

* **Minimum cutoff ($r_{min}$):** Set at the intersection point $r_{int}(t)$

* **Maximum cutoff ($r_{max}$):** Set at the "elbow" of the singular value curve (see algorithm 4)

The algorithm performs two natural gradient steps per iteration, one with each cutoff. If the intersection position remains static after both updates, we increment $r_{max}$ by one to promote exploration of additional gradient components.

This phase ends when $r_{min}$ reaches $r_{max}$—an event we term **liftoff**.

### C.2.2 Ascent Phase

During ascent, both cutoffs track the moving intersection, but with a stability mechanism:

$$r_{max}(t) = \max(r_{max}(t-1), r_{int}(t)). \qquad (21)$$

This monotonicity constraint prevents the intersection rank from falling to zero, which would disrupt training dynamics.

### C.2.3 Stage Separation and Precision Locking

When $RCE_{r_{int}}^S(t) \le \epsilon$, we trigger **stage separation**:

* $r_{min}$ is fixed at the precision level: $r_{min} = r_{\epsilon}(t)$

* $r_{max}$ continues tracking the intersection to maintain stability

The algorithm continues until $r_{min} = 0$ (**booster return**), indicating complete convergence. The final algorithm that combines all three stages is mentioned in Algorithm 2.
<caption>**Algorithm 1:** Sketch of the Adaptive MultiCutoff Strategy for ANaGRAM (AMStraMGRAM)</caption>

**Input:** $u_{\theta}: \mathbb{R}^P \to L^2(\Omega, \mu)$, $\theta_0 \in \mathbb{R}^P$, $f \in L^2(\Omega, \mu)$, $(x_i) \in \Omega^S$, $\epsilon > 0$, $T_{\text{max}} \in \mathbb{N}$
// Initialization
1 & $t \leftarrow 0$
2 & $\hat{\phi}_0 \leftarrow (\partial_p u_{\theta_0}(x_i))_{i,p}$ for $i \in 1, \dots, S$ and $p \in 1, \dots, P$
3 & $\hat{U}_0, \hat{\Delta}_0, \hat{V}_0^T \leftarrow \text{SVD}(\hat{\phi}_0)$
4 & $\widehat{\nabla} \mathcal{L}_0 \leftarrow (u_{\theta_0}(x_i) - f(x_i))_i$ for $i \in 1, \dots, S$
5 & Compute $(\text{RCE}_j^S)$ for all $j \in 1, \dots, r_{\text{svd}}$ following Equation (13)
6 & **repeat**
7 & & & // Compute adaptive ranks
8 & & & Compute $r_{\text{int}}$ and $r_{\epsilon}$ using expressions from Equation (20)
9 & & & // Determine a final cutoff rank
10 & & & if $\text{RCE}_{r_{\text{int}}}^S > \epsilon$ then
11 & & & & $r_{\text{cutoff}} \leftarrow r_{\text{int}}$
12 & & & & // Track intersection
13 & & & else
14 & & & $r_{\text{cutoff}} \leftarrow r_{\epsilon}$
15 & & & // Lock on precision
16 & & & // Natural gradient step
17 & & & Set $\hat{\Delta}_t \leftarrow \begin{cases} \hat{\Delta}_{t,i} & \text{if } i \le r_{\text{cutoff}}, \\ 0 & \text{else}; \end{cases}$
18 & & & // Update for next iteration
19 & & & $\hat{\phi}_{t+1} \leftarrow (\partial_p u_{\theta_{t+1}}(x_i))_{i,p}$
20 & & & $\hat{U}_{t+1}, \hat{\Delta}_{t+1}, \hat{V}_{t+1}^T \leftarrow \text{SVD}(\hat{\phi}_{t+1})$
21 & & & $\widehat{\nabla} \mathcal{L}_{t+1} \leftarrow (u_{\theta_{t+1}}(x_i) - f(x_i))_i$
22 & & & Recompute $\text{RCE}_j^S$ for all $j \in 1, \dots, r_{\text{svd}}$ following Equation (13)
23 & & & $t \leftarrow t + 1$
24 & & & **until** $r_{\epsilon} = 0$ or $t \ge T_{\text{max}}$
25 & & & **Output:** $\theta_t$
### C.3 Complete Practical Algorithm

**Algorithm 2:** AMStraMGRAM : Adaptive Multicutoff Strategy Modification for ANaGRAM

**Input:** $\bullet$ $u: \mathbb{R}^P \to L^2(\Omega, \mu)$ // neural network architecture

* $\theta_0 \in \mathbb{R}^P$ // initialization of the neural network

* $f \in L^2(\Omega, \mu)$ // target function of the quadratic regression

* $(x_i) \in \Omega^S$ // a batch in $\Omega$

* $\epsilon > 0$ // precision level of the optimization

1. **begin Initialization**
2. &nbsp;&nbsp;&nbsp;&nbsp;`$\lambda \leftarrow False$` // Liftoff indicator
3. &nbsp;&nbsp;&nbsp;&nbsp;`$\hat{\phi}_{\theta_0} \leftarrow (\partial_p u_{\theta_0}(x_i))_{1 \le i \le S, 1 \le p \le P}$` // Computed via auto-differentiation
4. &nbsp;&nbsp;&nbsp;&nbsp;`$\hat{U}_{\theta_0}, \hat{\Delta}_{\theta_0}, \hat{V}_{\theta_0}^t \leftarrow \text{SVD}(\hat{\phi}_{\theta_0})$`
5. &nbsp;&nbsp;&nbsp;&nbsp;`$\nabla \mathcal{L}_{\theta_0} \leftarrow (u_{\theta_0}(x_i) - f(x_i))_{1 \le i \le S}$`
6. &nbsp;&nbsp;&nbsp;&nbsp;`$\text{RCE}_0^S \leftarrow \text{ReconstructionErrors}(\hat{V}_{\theta_0}^t, \nabla \mathcal{L}_{\theta_0})$`
7. &nbsp;&nbsp;&nbsp;&nbsp;`$r_{\text{max}0} \leftarrow \text{FindElbow}((1, \dots, r_{\text{svd}}), \hat{\Delta}_{\theta_0})$`
8. **repeat**
9. &nbsp;&nbsp;&nbsp;&nbsp;10. &nbsp;&nbsp;&nbsp;&nbsp;`$r_{1t} \leftarrow \# \{\text{RCE}_{0j}^S \le \hat{\Delta}_{\theta_j} : 1 \le j \le r_{\text{svd}}\}$`
10. &nbsp;&nbsp;&nbsp;&nbsp;11. &nbsp;&nbsp;&nbsp;&nbsp;`$r_{2t} \leftarrow \# \{\text{RCE}_{0j}^S \ge \epsilon : 1 \le j \le r_{\text{svd}}\}$`
11. &nbsp;&nbsp;&nbsp;&nbsp;12. `/* with # standing for the cardinal */`
12. &nbsp;&nbsp;&nbsp;&nbsp;13. &nbsp;&nbsp;&nbsp;&nbsp;`$r_{\text{min}t} \leftarrow \min(r_{1t}, r_{2t})$`
13. &nbsp;&nbsp;&nbsp;&nbsp;14. &nbsp;&nbsp;&nbsp;&nbsp;`$r_{\text{max}t} \leftarrow \max(r_{1t}, r_{\text{max}t-1})$`
14. &nbsp;&nbsp;&nbsp;&nbsp;15. &nbsp;&nbsp;&nbsp;&nbsp;**if not $\lambda_t$ then**
15. &nbsp;&nbsp;&nbsp;&nbsp;16. &nbsp;&nbsp;&nbsp;&nbsp;| &nbsp;&nbsp;&nbsp;&nbsp;**if** $r_{\text{min}t} \ge r_{\text{max}t}$ **then**
16. &nbsp;&nbsp;&nbsp;&nbsp;17. &nbsp;&nbsp;&nbsp;&nbsp;| &nbsp;&nbsp;&nbsp;&nbsp;`$\lambda_t \leftarrow \text{True}$`
17. &nbsp;&nbsp;&nbsp;&nbsp;18. &nbsp;&nbsp;&nbsp;&nbsp;**else if** $r_{\text{min}t-1} = r_{\text{min}t}$ **then**
18. &nbsp;&nbsp;&nbsp;&nbsp;19. &nbsp;&nbsp;&nbsp;&nbsp;| &nbsp;&nbsp;&nbsp;&nbsp;`$r_{\text{max}t} \leftarrow r_{\text{max}t} + 1$`
19. &nbsp;&nbsp;&nbsp;&nbsp;20. &nbsp;&nbsp;&nbsp;&nbsp;**foreach** $r_{\text{cutoff}} \in \{r_{\text{max}t}, r_{\text{min}t}\}$ **do**
20. &nbsp;&nbsp;&nbsp;&nbsp;21. &nbsp;&nbsp;&nbsp;&nbsp;| &nbsp;&nbsp;&nbsp;&nbsp;`$\hat{\Delta}_{\theta_t} \leftarrow (\hat{\Delta}_{\theta_{t,p}} \text{ if } p \ge r_{\text{cutoff}} \text{ else } 0)_{1 \le p \le P}$`
21. &nbsp;&nbsp;&nbsp;&nbsp;22. &nbsp;&nbsp;&nbsp;&nbsp;| &nbsp;&nbsp;&nbsp;&nbsp;`$\nabla \mathcal{L}_{\theta_t} \leftarrow (u_{\theta_t}(x_i) - f(x_i))_{1 \le i \le S}$`
22. &nbsp;&nbsp;&nbsp;&nbsp;23. &nbsp;&nbsp;&nbsp;&nbsp;| &nbsp;&nbsp;&nbsp;&nbsp;`$d_{\theta_t} \leftarrow \hat{V}_{\theta_t} \hat{\Delta}_{\theta_t}^\dagger \hat{U}_{\theta_t}^t \nabla \mathcal{L}_{\theta_t}$`
23. &nbsp;&nbsp;&nbsp;&nbsp;24. &nbsp;&nbsp;&nbsp;&nbsp;| &nbsp;&nbsp;&nbsp;&nbsp;`$\eta_t \leftarrow \arg \min_{\eta \in \mathbb{R}^+} \sum_{1 \le i \le S} (f(x_i) - u_{\theta_t} - \eta d_{\theta_t}(x_i))^2$` // via line search
24. &nbsp;&nbsp;&nbsp;&nbsp;25. &nbsp;&nbsp;&nbsp;&nbsp;| &nbsp;&nbsp;&nbsp;&nbsp;`$\theta_{t+1} \leftarrow \theta_t - \eta_t d_{\theta_t}$`
25. &nbsp;&nbsp;&nbsp;&nbsp;26. &nbsp;&nbsp;&nbsp;&nbsp;| &nbsp;&nbsp;&nbsp;&nbsp;`$\hat{\phi}_{\theta_{t+1}} \leftarrow (\partial_p u_{\theta_{t+1}}(x_i))_{1 \le i \le S, 1 \le p \le P}$` // Computed via auto-differentiation
26. &nbsp;&nbsp;&nbsp;&nbsp;27. &nbsp;&nbsp;&nbsp;&nbsp;| &nbsp;&nbsp;&nbsp;&nbsp;`$\hat{U}_{\theta_{t+1}}, \hat{\Delta}_{\theta_{t+1}}, \hat{V}_{\theta_{t+1}}^t \leftarrow \text{SVD}(\hat{\phi}_{\theta_{t+1}})$` // via line search

### C.4 Empirical Justification for Design Choices

The dual cutoff strategy addresses specific empirical challenges we observed:

**Dual gradient steps:** Without the second cutoff, training dynamics sometimes stagnate. The dual approach provides both stability (via $r_{\text{min}}$) and exploration (via $r_{\text{max}}$).

**Elbow initialization:** The elbow point marks where singular values cease contributing meaningful signal, providing a natural upper bound for exploration.

**Monotonic $r_{\text{max}}$:** Prevents catastrophic retreat of the intersection point, which we observed in complex equations like Allen-Cahn.
**Stage separation timing:** Triggered precisely when the intersection error drops below target precision, ensuring optimal utilization of the flattening phenomenon.

We see in the next section how this practical algorithm successfully improve empirical robustness.

## D Algorithmic details

<table><caption>**Algorithm 3: Find elbow**</caption><thead><tr><th></th><th></th><th></th></tr></thead><tbody><tr><td>1</td><td>Function FindElbow</td><td></td></tr><tr><td></td><td>Input:- ($x_i \in \mathbb{R}^m$ // an increasing sequence of $m \in \mathbb{N}$ points in $\mathbb{R}$</td><td></td></tr><tr><td></td><td>- $\hat{f} \in \mathbb{R}^m$ // a decreasing function evaluated at points ($x_i$)</td><td></td></tr><tr><td></td><td>/* Clockwise normal vector to ($x_m - x_1, \hat{f}_m - \hat{f}_1$) */</td><td></td></tr><tr><td>2</td><td>$\vec{n} \leftarrow (\hat{f}_m - \hat{f}_1, x_1 - x_m) \in \mathbb{R}^2$</td><td></td></tr><tr><td>3</td><td>($s_j$)$_{1 \le j \le m} \leftarrow (\langle \vec{n}, (x_j - x_1, \hat{f}_j - \hat{f}_1) \rangle_{\mathbb{R}^2})_{1 \le j \le m}$</td><td></td></tr><tr><td></td><td>Output: $\arg\max_{1 \le j \le m} s_j$</td><td></td></tr><tr><td>4</td><td>end</td><td></td></tr></tbody></table>

<table><caption>**Algorithm 4: Reconstruction Errors**</caption><thead><tr><th></th><th></th><th></th></tr></thead><tbody><tr><td>1</td><td>Function ReconstructionErrors</td><td></td></tr><tr><td></td><td>Input:- $\hat{V}^t \in \mathbb{R}^{r_{svd}, S}$ // right singular vectors of the Jacobian $\hat{\phi}$</td><td></td></tr><tr><td></td><td>- $\nabla \mathcal{L} \in \mathbb{R}^S$ // Evaluated functional gradient</td><td></td></tr><tr><td>2</td><td>begin Initialization</td><td></td></tr><tr><td>3</td><td>$\hat{\Sigma} \leftarrow 0 \in \mathbb{R}^S$ // cumulative approximation of $\nabla \mathcal{L}$</td><td></td></tr><tr><td>4</td><td>$\text{RCE}^S \leftarrow 0 \in \mathbb{R}^{r_{svd}}$ // cumulated reconstruction errors</td><td></td></tr><tr><td>5</td><td>$\hat{c} \leftarrow \hat{V}^t \nabla \mathcal{L} \in \mathbb{R}^{r_{svd}}$</td><td></td></tr><tr><td>6</td><td>end</td><td></td></tr><tr><td>7</td><td>foreach $j \in (1, \dots, r_{svd})$ do</td><td></td></tr><tr><td>8</td><td>$\hat{\Sigma} \leftarrow \hat{\Sigma} + \hat{c}_j$</td><td></td></tr><tr><td>9</td><td>$\text{RCE}_j^S \leftarrow \|\hat{\Sigma} - \hat{c}_j\|_2$</td><td></td></tr><tr><td>10</td><td>end</td><td></td></tr><tr><td>11</td><td>end</td><td></td></tr></tbody></table>

## E Empirical example of Anagram Training Dynamics

In Figure 5, we analyze ANaGRAM's training on the heat equation with a fixed cutoff threshold $\alpha = 10^{-3}$ and line search for the learning rate. The training loss coincides with $\|\nabla \mathcal{L}\|^2$. We can see the flattening phenomenon to occur on Iteration 120 and completed at 150. As discussed in the main paper, sometimes the flattening can be incomplete, and for many iterations remain without any further progress ($N_{flat}$ never reaching zero). In this case, changing a cutoff threshold results in an immediate and complete flattening for all first components up to $r_{cutoff}$, which is demonstrated in Figure 6 for Iteration 120 of Figure 5.
<figure><img src="image_7.png" /><figcaption>(a) Iteration 0: intersection point between singular values and reconstruction error lies before cutoff.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(b) Iteration 40: intersection point shifts rightward toward cutoff.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(c) Iteration 90: intersection point passes the cutoff threshold.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(d) Iteration 120. Beginning of flattening: reconstruction errors stabilizes at constant level before cutoff.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(e) Iteration 150: Complete flattening. Training loss reaches the flattened reconstruction error level.</figcaption></figure>

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

**Figure 5:** Evolution of quantities of interest during ANaGRAM training on heat equation. The dynamics reveal two distinct phases culminating in reconstruction error flattening.
<figure><img src="image_13.png" /><figcaption>Figure 6: Illustration of “instant flattening” through adaptive cutoff adjustment. A single step with adjusted cutoff completes the flattening process.</figcaption></figure>

# F Deep dive on selected experiments

In this section we look at curves of training and estimations obtained with AMStraMGRAM on benchmark of PDEs.

## F.1 One Dimensional Burgers Equation

<figure><img src="image_14.png" /><figcaption>Figure 7: Training metrics for the One-Dimensional Burgers equation, showing convergence behavior with our adaptive multi-cutoff strategy.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 8: Results for One Dimensional Burgers Equation with cutoff 10<sup>-6</sup>.</figcaption></figure>

## F.2 Heat Equation

<figure><img src="image_16.png" /><figcaption>Figure 9: Convergence results for the Heat equation showing the $L_2$ error over iterations. Our method (AMStraMGRAM) converges faster and reaches a lower final error than ANaGRAM and baselines. Variability across runs is due to differing feature development speed from the random initialization.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 10: Results for the Heat equation (solution cutoff 10<sup>-14</sup>). The error remains uniformly low over the domain, illustrating the effectiveness of the adaptive multi-cutoff strategy.</figcaption></figure>

## F.3 Laplace Equations (L2D and L5D)

For the Laplace equation in 2D, our method also demonstrates remarkable performance improvements over the baselines. The convergence is both faster and reaches a significantly lower error plateau.
<figure><img src="image_18.png" /><figcaption>Figure 11: Convergence results for the Laplace 2D problem, showing the $L_2$ error over iterations. Our method (AMStraMGRAM) achieves both faster convergence and lower final error compared to ANaGRAM and other baseline methods. The observed variance between runs can be explained by different speed of convergence depending on the initialization.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 12: Results for Laplace 2D Equation with cutoff $10^{-6}$.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 13: Convergence results for the Laplace 5D problem, showing the $L_2$ error over iterations. Our method (AMStraMGRAM) achieves faster convergence but not lower final error compared to ANaGRAM and other baseline methods. We see that seeds change the speed of convergence of the algorithm.</figcaption></figure>

## F.4 Non Linear Poisson Equation

To compare ourselves with Urbán et al. (2025), we select $K=1$.
<figure><img src="image_21.png" /><figcaption>Figure 14: Convergence results for the Non Linear Poisson equation, showing the $L_2$ error over iterations. Our method (AMStraMGRAM) achieves both faster convergence and lower final error compared to ANaGRAM and other baseline methods.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure 15: Results for the Nonlinear Poisson equation (cutoff $10^{-4}$).</figcaption></figure>

## F.5 Allen-Cahn Equation

<figure><img src="image_23.png" /><figcaption>Figure 16: Training curves for the Allen-Cahn equation, showing the evolution of loss and error over iterations.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>Figure 17: Results on the Allen-Cahn equation, showing the error distribution (left), model prediction (middle), and true solution (right). The error is mostly present in regions with the "sharpest" transitions, which exemplifies the challenge of accurately capturing sharp interfaces still remains even for our advanced optimization approach.</figcaption></figure>

# G Geometrical interpretation of regularizations

## G.1 Why Regularization is Necessary

We recall that our goal is to solve the operator equation $D[u] = f$ by minimizing the squared residual

$$ \|D[u] - f\|_{L^2(\Omega, \mu)}^2. \qquad (22) $$

For simplicity, assume $D$ is linear. Then the mapping

$$ u \in C^{\infty}(\Omega) \mapsto \|D[u]\|_{L^2(\Omega, \mu)} \qquad (23) $$

defines a semi-norm on $C^{\infty}(\Omega)$. We can “upgrade” this semi-norm into a true norm by introducing the following generalized Sobolev norm:

$$ \|\cdot\|_{\tilde{H}_D} : \begin{cases} C^{\infty}(\Omega \to \mathbb{R}) & \to \mathbb{R}^+ \\ u & \mapsto \sqrt{\|u\|_{L^2(\Omega \to \mathbb{R}, \mu)}^2 + \|D[u]\|_{L^2(\Omega \to \mathbb{R}, \mu)}^2} \end{cases} \qquad (24) $$

Clearly, for any $u$,

$$ \|u\|_{L^2(\Omega, \mu)} \le \|u\|_{\tilde{H}_D}, \qquad (25) $$

which guarantees that $\|\cdot\|_{\tilde{H}_D}$ is definite, i.e. $\|u\|_{\tilde{H}_D} = 0 \iff u = 0$.

Completing $C^{\infty}(\Omega)$ with respect to $\|\cdot\|_{\tilde{H}_D}$ yields a generalized Sobolev space $(\mathcal{H}_D, \|\cdot\|_{\mathcal{H}_D})$. This Hilbert space is the largest subspace of $L^2(\Omega, \mu)$ on which $D$ is continuous. Indeed, for every $u \in \mathcal{H}_D$,

$$ \|D[u]\|_{L^2(\Omega, \mu)} \le \|u\|_{\mathcal{H}_D}. \qquad (26) $$

Since our goal is to solve $D[u] = f$, we need $D$ to be continuously invertible. That is, we need the reverse inequality of Equation (26) to hold (up to a constant $\alpha > 0$). Formally, if $D$ were algebraically invertible (bijective as a mapping), this condition would read:

$$ \begin{align} & \left( \exists \alpha > 0, \forall u \in \mathcal{H}_D, \|u\|_{\mathcal{H}_D} \le \alpha \|D[u]\|_{L^2(\Omega \to \mathbb{R}, \mu)} \right) \notag \\ \iff & \left( \exists \alpha > 0, \forall u \in \mathcal{H}_D, \|D^{-1}[D[u]\|_{\mathcal{H}_D} \le \alpha \|D[u]\|_{L^2(\Omega \to \mathbb{R}, \mu)} \right) . \tag{27} \\ \iff & \left( \exists \alpha > 0, \forall f \in L^2(\Omega \to \mathbb{R}, \mu), \|D^{-1}[f]\|_{\mathcal{H}_D} \le \alpha \|f\|_{L^2(\Omega \to \mathbb{R}, \mu)} \right) \notag \end{align} $$

**Operator ill-conditioning.** Even if $D$ is bijective, Equation (27) may fail to hold, i.e. $D$ can be ill-conditioned. Suppose there exists a subspace $\mathcal{H}_K \subset \mathcal{H}_D$ such that $D$ acts compactly on $\mathcal{H}_K$ with infinite rank. Then $D$ admits a singular value decomposition (Kress, 2014, Theorem 15.16): for $u \in \mathcal{H}_K$,

$$ D[u] = \sum_{n \in \mathbb{N}} e_n \lambda_n \langle v_n, u \rangle_{\mathcal{H}_D}, \qquad (28) $$

with $(v_n)$ orthonormal in $\mathcal{H}_D$, $(e_n)$ orthonormal in $L^2(\Omega, \mu)$, and $\lambda_n \to 0$ as $n \to \infty$.

For Equation (27) to hold, we would need $\inf_n \lambda_n > 0$, contradicting $\lambda_n \to 0$. This is exactly the classical inverse problem setting: $D$ is bijective but ill-conditioned, and regularization is unavoidable. Among the many schemes developed, Tikhonov regularization is the canonical example (Kirsch, 2021).
**Non-bijectivity.** If $D$ is not bijective, two additional issues may occur.

**Non-surjectivity.** If $\mathrm{Im} D$ is a closed subspace, we can still obtain a solution by replacing the target $f$ with its projection $\Pi_{\mathrm{Im} D} f$. Note that minimizing $\|D[u] - f\|_{L^2(\Omega, \mu)}^2$ yields precisely this least-squares solution.

**Non-injectivity.** The lack of injectivity is a much more subtle issue. Since $D$ is linear and continuous, its null space $\mathrm{Ker} D$ is a closed subspace of $\mathcal{H}_D$. In principle, one could restrict the domain of $D$ to $\mathrm{Ker} D^{\perp}$ to make it injective. The problem, however, is that identifying $\mathrm{Ker} D$ is typically just as hard as solving the original problem itself, since it amounts to characterizing all $u \in \mathcal{H}_D$ such that $D[u] = 0$. Therefore, unless one can rely on theoretical results that explicitly describe $\mathrm{Ker} D$, or construct a subspace $\mathcal{H}_0 \subset \mathcal{H}_D$ for which $\mathrm{Ker} D \cap \mathcal{H}_0$ is explicitly known (so that $D$ can be restricted to $\mathcal{H}_0$), it is generally impossible to “get rid of” $\mathrm{Ker} D$ in practice.

On the other hand, if we do not filter out $\mathrm{Ker} D$, this has the unwanted consequence of introducing “spurious” low-energy signals. To be concrete, suppose we approximate our solution in a space $\mathcal{H}_K$ with orthonormal basis $(u_n)_{n \in \mathbb{N}}$. Assume there exists a subsequence $(u_n^S) \notin \mathrm{Ker} D$ converging towards $\mathrm{Ker} D$. Since $\mathrm{Ker} D$ is closed (by continuity of $D$), this means

$$ \lim_{n \to \infty} \| \Pi_{\mathrm{Ker} D} u_n^S - u_n^S \|_{\mathcal{H}_D}^2 = 0. \quad (29) $$

Equivalently, after extraction, this can be rewritten for all $n \in \mathbb{N}$ as

$$ \frac{\| \Pi_{\mathrm{Ker} D} u_n^S \|_{\mathcal{H}_D}^2}{\| u_n^S \|_{\mathcal{H}_D}^2} \ge 1 - 2^{-n}. \quad (30) $$

Now consider normalized vectors $u_n^S / \|u_n^S\|_{\mathcal{H}_D}$. We have

$$ \begin{aligned} 0 < \|D\left[\frac{u_n^S}{\|u_n^S\|_{\mathcal{H}_D}}\right]\|_{\mathcal{H}_D}^2 &= \left\|D\left[\frac{\Pi_{\mathrm{Ker} D \perp} u_n^S + \Pi_{\mathrm{Ker} D} u_n^S}{\|u_n^S\|_{\mathcal{H}_D}}\right]\right\|_{\mathcal{H}_D}^2 \\ &= \left\|D\left[\frac{\Pi_{\mathrm{Ker} D \perp} u_n^S}{\|u_n^S\|_{\mathcal{H}_D}}\right] + D\left[\frac{\Pi_{\mathrm{Ker} D} u_n^S}{\|u_n^S\|_{\mathcal{H}_D}}\right]\right\|_{\mathcal{H}_D}^2 \\ &= \left\|D\left[\frac{\Pi_{\mathrm{Ker} D \perp} u_n^S}{\|u_n^S\|_{\mathcal{H}_D}}\right]\right\|_{\mathcal{H}_D}^2 \\ &\stackrel{(26)}{\le} \left\|\frac{\Pi_{\mathrm{Ker} D \perp} u_n^S}{\|u_n^S\|_{\mathcal{H}_D}}\right\|_{\mathcal{H}_D}^2 \\ &= 1 - \frac{\|\Pi_{\mathrm{Ker} D} u_n^S\|_{\mathcal{H}_D}^2}{\|u_n^S\|_{\mathcal{H}_D}^2} \\ &\stackrel{(30)}{\le} 2^{-n}. \end{aligned} \quad (31) $$

In particular, if (for simplicity) the normalized $(u_n/\|u_n\|_{\mathcal{H}_D})$ are right singular vectors of $D$, then the vectors $(u_n^S/\|u_n^S\|_{\mathcal{H}_D})$ will correspond to singular values vanishing at least as fast as $(2^{-n})$. Crucially, however, these vanishing singular values do not reflect an intrinsic ill-conditioning of $D$, but rather an *artificial* ill-conditioning induced by the choice of approximation space $\mathcal{H}_K$. In other words, the spurious instability arises from how we approximate the operator, not from the operator itself. For more details on this approximation-induced phenomenon, see Adcock & Huybrechs (2019, 2020).

These remarks highlight the *inevitable need for regularization* in practice. In the next section, we will provide a geometric interpretation of the two regularization schemes introduced in Section 2.5, emphasizing how fundamentally different they are in nature.

**Remark 4.** The above discussion becomes even more critical when we restrict ourselves to a finite-dimensional approximation space $\mathcal{H}_{\mathrm{app}} \subset \mathcal{H}_D$. In this case, the restriction $D_{\mathrm{app}}$ is automatically compact, since it is of finite rank. As a consequence, both types of ill-conditioning described above may occur simultaneously. This highlights once again why regularization is not merely convenient but unavoidable in numerical practice.
## G.2 Ridge-regression

Returning to the definition given in Section 2.5, recall that Ridge regression amounts to adding $\alpha^2 I_d$ (for some $\alpha > 0$) to the Gram matrix $G_\theta$ introduced in Equation (9):

$$ \theta_{t+1} \leftarrow \theta_t - \eta G_{\theta_t}^{\dagger} \nabla \ell(\theta_t); \qquad G_{\theta_t, p, q} := \langle \partial_p u_{\theta_t}, \partial_q u_{\theta_t} \rangle_{L^2(\Omega, \mu)}. \qquad (9) $$

We can reformulate this observation in the following way: given our model

$$ u : \mathbb{R}^P \to L^2(\Omega \to \mathbb{R}, \mu), \qquad (32) $$

consider the regularized model

$$ u^{\alpha} : \left\{ \begin{array}{ll} \mathbb{R}^P & \to L^2(\Omega, \mu) \times \mathbb{R}^P \\ \boldsymbol{\theta} & \mapsto (u_{\boldsymbol{\theta}}, \alpha \boldsymbol{\theta}). \end{array} \right. \qquad (33) $$

The Gram matrix of this regularized model is exactly $G_\theta + \alpha^2 I_d$. Suppose further that regression is performed with respect to some function $f \in L^2(\Omega, \mu)$. Then we must adapt the objective to the regularized model, replacing $f$ with the pair

$$ (f, \alpha \boldsymbol{\theta}) \in L^2(\Omega, \mu) \times \mathbb{R}^P. \qquad (34) $$

A straightforward computation shows that, for all $1 \le p \le \min(P, S)$,

$$ \begin{aligned} \langle \partial_p u_{\boldsymbol{\theta}}^{\alpha}, (f, \alpha \boldsymbol{\theta}) - u_{\boldsymbol{\theta}}^{\alpha} \rangle_{L^2(\Omega, \mu) \times \mathbb{R}^P} &= \langle \partial_p u_{\boldsymbol{\theta}}, f - u_{\boldsymbol{\theta}} \rangle_{L^2(\Omega, \mu)} + \alpha \langle e^{(p)}, \boldsymbol{\theta} - \boldsymbol{\theta} \rangle_{\mathbb{R}^P} \\ &= \langle \partial_p u_{\boldsymbol{\theta}}, f - u_{\boldsymbol{\theta}} \rangle_{L^2(\Omega, \mu)}. \end{aligned} \qquad (35) $$

Thus, regression of $(f, \alpha \boldsymbol{\theta})$ with the regularized model is exactly equivalent to Ridge regression. Equivalently, Ridge regression corresponds to replacing the original model $u$ by the regularized model $u^\alpha$, and replacing the objective $f$ by $(f, \alpha \boldsymbol{\theta})$. From this point of view, the choice of $\alpha \boldsymbol{\theta}$ as the secondary target may be interpreted as a default assumption in the absence of prior information on the parameters: one simply uses the current parameters as a reference target.

We can now extract several fundamental facts:

1. As $\alpha \to 0$, the regularized model $u^\alpha$ tends in operator norm to the unregularized model $(u, 0)$ (i.e. $u$ by abuse of notation). Indeed,

$$ \sup_{\|\boldsymbol{\theta}\|_{\mathbb{R}^P}=1} \|du_{\boldsymbol{\theta}}^{\alpha} - (du_{\boldsymbol{\theta}}, 0)\|_{L^2(\Omega, \mu) \times \mathbb{R}^P} = \alpha \sup_{\|\boldsymbol{\theta}\|_{\mathbb{R}^P}=1} \|\boldsymbol{\theta}\|_{\mathbb{R}^P} = \alpha. \qquad (36) $$

2. The model $u^\alpha$ is injective and continuous. Since $du_\theta$ is continuous (as $\mathbb{R}^P$ is finite-dimensional), the only possible source of non-injectivity is Ker $du_\theta^\alpha$. But

$$ \mathrm{Ker} \, du_{\boldsymbol{\theta}}^{\alpha} = \mathrm{Ker} \, du_{\boldsymbol{\theta}} \cap \mathrm{Ker}(\alpha I_{\mathbb{R}^P}) \subset \mathrm{Ker}(\alpha I_{\mathbb{R}^P}) = \{0\}, \qquad (37) $$

hence injectivity. Restricting $u^\alpha$ to its image makes it algebraically bijective, and the inverse is continuous since

$$ \alpha \|\boldsymbol{\theta}\|_{\mathbb{R}^P} \le \|du_{\boldsymbol{\theta}}^{\alpha}\|_{L^2(\Omega, \mu) \times \mathbb{R}^P}. \qquad (38) $$

By the equivalence stated in Equation (27), this implies that $(du_\theta^\alpha)^{-1}$ is continuous. Consequently, $\mathrm{Im} \, du_\theta^\alpha$ is closed in $L^2(\Omega, \mu) \times \mathbb{R}^P$, since it is the inverse image of a closed set under $(du_\theta^\alpha)^{-1}$. Therefore least-squares solution is well-defined.

3. The least-squares solution of $u^\alpha = (f, 0)$ is influenced by $\alpha$ as follows: $(f, 0)$ is projected onto

$$ \mathrm{Im} \, du_{\boldsymbol{\theta}}^{\alpha} = \mathrm{Span} \left( (\partial_p u_{\boldsymbol{\theta}}, \alpha e^{(p)}) : 1 \le p \le P \right). \qquad (39) $$

In particular, even if $f \in \mathrm{Im} \, du_\theta$ and $f \ne 0$, we still have $(f, 0) \notin \mathrm{Im} \, du_\theta^\alpha$ (since $du_\theta(0) = 0$). Consequently,

$$ (\Pi_{\mathrm{Im} \, du_{\boldsymbol{\theta}}^\alpha}^{\perp} (f, 0))_1 \ne f, \qquad (40) $$

where the subscript 1 denotes projection onto the first component in $L^2(\Omega, \mu) \times \mathbb{R}^P$.
We illustrate these phenomena in Figure 18a.

Building on the above analysis, we now show that Ridge regression can be extended to the functional setting. To this end, let us reconsider the operator $D: \mathcal{H}_D \to L^2(\Omega, \mu)$ introduced in Section G.1. Analogously to what we did for the parametric model $u$, we define the regularized operator at level $\alpha > 0$ as

$$D^{\alpha}: \left\{ \begin{array}{ll} \mathcal{H}_D & \to L^2(\Omega, \mu) \times \mathcal{H}_D \\ u & \mapsto (D[u], \alpha u) \end{array} \right. . \qquad (41)$$

The corresponding target becomes the regularized objective $(f, \alpha u)$.

At this level of generality, the equivalence with Gram-matrix regularization no longer holds, since we are dealing with infinite-dimensional operators for which no direct Gram-matrix representation exists. Nevertheless, the fundamental properties remain valid, namely:

1. When $\alpha \to 0$, the regularized operator $D^\alpha$ converges to $(D, 0)$ in the operator-norm sense, i.e. to $D$ by a mild abuse of notation. Indeed, we have

$$\sup_{\|u\|_{\mathcal{H}_D}=1} \|D^{\alpha}[u] - (D, 0)\|_{L^2(\Omega, \mu) \times \mathcal{H}_D} = \alpha \sup_{\|u\|_{\mathcal{H}_D}=1} \|u\|_{\mathcal{H}_D} = \alpha. \qquad (42)$$

2. The operator $D^\alpha$ is injective and continuous. Indeed, $D$ is continuous by the very construction of $\mathcal{H}_D$ (see Section G.1), and injectivity follows since

$$\mathrm{Ker} D^{\alpha} = \mathrm{Ker} D \cap \mathrm{Ker}(\alpha I_{\mathcal{H}_D}) \subseteq \mathrm{Ker}(\alpha I_{\mathcal{H}_D}) = \{0\}. \qquad (43)$$

Restricting $D^\alpha$ to its image makes it algebraically bijective, and the inverse is continuous: we have $\alpha \|u\|_{\mathcal{H}_D} \le \|D^\alpha[u]\|_{L^2(\Omega, \mu) \times \mathcal{H}_D}$, which by the equivalence in Equation (27) implies that $(D^\alpha)^{-1}$ is continuous. Consequently, $\mathrm{Im} D^\alpha$ is closed in $L^2(\Omega, \mu) \times \mathcal{H}_D$, since it is the inverse image of a closed set under $(D^\alpha)^{-1}$. Therefore least-squares solution is well-defined.

3. Least-squares solutions of the regularized problem $D^\alpha[u] = (f, 0)$ are impacted by $\alpha$ in the following way: we are projecting $(f, 0)$ onto

$$\mathrm{Im} D^{\alpha} = \mathrm{Span} \left( (D[h], \alpha h) : h \in \mathcal{H}_D \right). \qquad (44)$$

In particular, even if $f \in \mathrm{Im} D$ with $f \ne 0$, we have $(f, 0) \notin \mathrm{Im} D^\alpha$ (since $D[0] = 0$), and hence

$$(\Pi_{\mathrm{Im} D^{\alpha}}^{\perp}(f, 0))_1 \ne f, \qquad (45)$$

where the subscript 1 denotes the first coordinate in $L^2(\Omega, \mu) \times \mathcal{H}_D$.

We illustrate these phenomena in Figure 18b.

In summary, Ridge regression can be interpreted as a modification of the operator $D$, rendering it injective and continuously invertible on its image. However, this comes at a price: the regularized solutions are never exact solutions of the original equation $D[u] = f$, even when $\alpha$ is arbitrarily small, since we are in fact solving a different operator equation. This marks a fundamental distinction from cutoff regularization, which instead acts directly on the approximation space, as we shall see in the next section.
<figure><img src="image_25.png" /><figcaption>(a) Illustration of parametric Ridge regression.</figcaption></figure>

The green region represents the solution space, while the blue regions denote the target spaces. As $\alpha \to 0$, the regularized graph $\Gamma_{du_{\theta}^\alpha}$ of $du_{\theta}^\alpha$ approaches the graph $\Gamma_{du_{\theta}}$ of $du_{\theta}$, with the angle between them vanishing at rate $\arctan(\alpha)$. The key consequence is that the projection of the objective $f$ onto $\mathrm{Im} \, du_{\theta}^\alpha$ follows a non-linear path as $\alpha \to 0$, coinciding with $\Pi_{\mathrm{Im} \, du_{\theta}} f$ only asymptotically.

<figure><img src="image_26.png" /><figcaption>(b) Illustration of functional Ridge regression.</figcaption></figure>

The green region represents the solution space, while the blue regions denote the target spaces. As $\alpha \to 0$, the regularized graph $\Gamma_{D^\alpha}$ of $D^\alpha$ approaches the graph $\Gamma_{\mathcal{H}_D}$ of $D$, with the angle between them vanishing at rate $\arctan(\alpha)$. The key consequence is that the projection of the objective $f$ onto $\mathrm{Im} \, D^\alpha$ follows a non-linear path as $\alpha \to 0$, coinciding with $\Pi_{\mathrm{Im} \, D} f$ only asymptotically.

<caption>Figure 18: Illustrations of Ridge regression.</caption>

## G.3 Cutoff regression

As in Section G.2, let us return to the setting of Section 2.5. In Equation (11), we introduced cutoff regularization from the SVD perspective: given the differential $du_{\theta}$ of the model $u$, at the point $\theta$, and its singular value decomposition $du_{\theta} = V_{\theta} \Delta_{\theta} U_{\theta}^T$, the cutoff-regularized pseudo-inverse $du_{\theta}^{\dagger\alpha}$ at level $\alpha > 0$ is defined as

$$du_{\theta}^{\dagger\alpha} := U_{\theta} \Delta_{\theta}^{\dagger\alpha} V_{\theta}^{T}; \qquad \Delta_{\theta,p}^{\dagger\alpha} := \begin{cases} \Delta_{\theta,p}^{-1} & \text{if } \Delta_{\theta,p} \ge \alpha, \\ 0 & \text{otherwise} \end{cases}, \quad 1 \le p \le P. \qquad (46)$$

Let us reinterpret this construction. Denote by $N_{\alpha} \in \mathbb{N}$ the number of singular values larger than $\alpha$. Equivalently, assuming $(\Delta_{\theta,p})_{1 \le p \le P}$ is non-increasing,

$$N_{\alpha} := \arg \max_{p \in \mathbb{N}} \{ \Delta_{\theta,p} \ge \alpha \}. \qquad (47)$$

Define

$$\Theta_{\alpha} := \mathrm{Span}\{U_{\theta,p} : 1 \le p \le N_{\alpha}\}, \quad T_{N_{\alpha}}^{0} \mathcal{M} := \mathrm{Span}\{V_{t,p} : 1 \le p \le N_{\alpha}\}, \qquad (48)$$

so that $T_{N_{\alpha}}^{0} \mathcal{M} = du_{\theta}(\Theta_{\alpha})$. We then have

$$(du_{\theta|_{\Theta_{\alpha}}}^{T_{N_{\alpha}} \mathcal{M}})^{-1} = du_{\theta}^{\dagger\alpha}, \qquad (49)$$

meaning that the restriction $du_{\theta}^{\alpha} := du_{\theta|_{\Theta_{\alpha}}}$ of $du_{\theta}$ to the domain $\Theta_{\alpha}$ becomes invertible once its codomain is restricted to its image $T_{N_{\alpha}}^{0} \mathcal{M}$, with inverse given precisely by the cutoff pseudo-inverse $du_{\theta}^{\dagger\alpha}$. Moreover, for any $h \in \Theta_{\alpha}$,

$$\|du_{\theta}(h)\|_{L^2(\Omega, \mu)} = \|V_{\theta} \Delta_{\theta} U_{\theta}^T h\|_{L^2(\Omega, \mu)} \stackrel{\text{unitary}}{=} \|\Delta_{\theta} U_{\theta}^T h\|_{\mathbb{R}^P} \stackrel{h \in \Theta_{\alpha}}{\ge} \alpha \|U_{\theta}^T h\|_{\mathbb{R}^P} \stackrel{\text{unitary}}{=} \alpha \|h\|_{\mathbb{R}^P}. \qquad (50)$$

In other words, Equation (27) is satisfied by $du_{\theta}^{\alpha}$.
Thus, while ridge regularization modifies the model itself, cutoff regularization instead restricts the domain of the model so that, on this restricted domain, Equation (27) holds and the model becomes invertible. We summarize the fundamental properties:

1. We have
$$ \bigcap_{\alpha>0} (\mathbb{R}^P \setminus \Theta_\alpha) = \text{Ker } du_\theta, \qquad (51) $$

that is, $\lim_{\alpha \to 0} \mathbb{R}^P \setminus \Theta_\alpha = \text{Ker } du_\theta$, since for all $\alpha > \beta$ we have $\Theta_\alpha \subset \Theta_\beta$ and then $\mathbb{R}^P \setminus \Theta_\beta \subset \mathbb{R}^P \setminus \Theta_\alpha$. Similarly, $\lim_{\alpha \to 0} T_{N_\alpha}^0 \mathcal{M} = \text{Im } du_\theta$. Moreover, for each $\alpha > 0$, the restriction $du_\theta^\alpha$ coincides with $du_\theta$ on $\Theta_\alpha$.

2. By Equation (50), $du_\theta^\alpha$ is injective and continuous. Restricting it to its image $T_{N_\alpha}^0 \mathcal{M}$ makes it bijective and bicontinuous, with inverse exactly the cutoff pseudo-inverse $du_\theta^{\alpha'}$. In particular $du(\Theta_\alpha)$ is closed in $L^2(\Omega, \mu)$, since it is the inverse image of a closed set under $du_\theta^{\alpha'}$. Therefore least-squares solution is well-defined.

3. Solving the least-squares problem $du_\theta^\alpha = f$ is now altered in the following way: the target $f$ is first projected onto $T_{N_\alpha}^0 \mathcal{M} = \text{Im } du_\theta^\alpha$. In particular, if for some $\alpha > 0$ we already have $f \in \text{Im } du_\theta^\alpha$, then the regularized least-squares formulation recovers an *exact solution* to the problem. This stands in sharp contrast with Ridge regression, where such exact recovery can only occur *asymptotically* in the limit $\alpha \to 0$.

As in Section G.2, we now need to reinterpret the cutoff regularization in order to extend it to the functional setting. Let us return once more to the operator $D : \mathcal{H}_D \to L^2(\Omega, \mu)$ introduced in Section G.1. In general, one cannot define an SVD for such an operator (except when it is compact). We must therefore appeal to the spectral theorem for bounded self-adjoint operators, which relies on the notion of a *projection-valued measure* (also called a resolution of the identity). For our purposes, it will be sufficient to simply state the definition.

**Definition 1** (Projection-valued measure). Let $(X, \mathcal{A})$ be a measurable space, where $\mathcal{A}$ denotes its $\sigma$-algebra, and let $\mathcal{H}$ be a Hilbert space. A *projection-valued measure* (PVM) is a map

$$ \pi : \mathcal{A} \to \mathcal{L}_b(\mathcal{H} \to \mathcal{H}), $$

where $\mathcal{L}_b(\mathcal{H} \to \mathcal{H})$ denotes the set of bounded operators on $\mathcal{H}$, such that for every $A \in \mathcal{A}$, $\pi(A)$ is an orthogonal projection on $\mathcal{H}$, and the following properties hold:

1. $\pi(\emptyset) = 0$ and $\pi(X) = I_\mathcal{H}$, where $I_\mathcal{H}$ is the identity operator on $\mathcal{H}$;

2. $\pi(A \cap B) = \pi(A)\pi(B)$ for all $A, B \in \mathcal{A}$;

3. For every countable family $(A_i)_{i=1}^\infty$ of disjoint sets in $\mathcal{A}$,

$$ \pi\left(\bigcup_{i=1}^\infty A_i\right) = \sum_{i=1}^\infty \pi(A_i), $$

where the series converges in the strong operator topology.

Since projection-valued measures are measures, one can define integrals with respect to them. We refer to (Berezansky et al., 1996, Chapter 13) for details. We may now state the spectral theorem.

**Theorem 2.** Let $\mathcal{H}$ be a Hilbert space and let $A : \mathcal{H} \to \mathcal{H}$ be a self-adjoint operator. Then there exists a projection-valued measure $\pi$ on the Borel $\sigma$-algebra of $\mathbb{R}$ such that

$$ A = \int_{\mathbb{R}} \lambda \pi(d\lambda) = \int_{\sigma(A)} \lambda \pi(d\lambda), \qquad (52) $$

where $\sigma(A)$ denotes the spectrum of $A$.

A proof can be found in (Berezansky et al., 1996, Theorem 4.1, Section 4.1, Chapter 13). In particular, since $\pi$ is a projection-valued measure, we have by Definition 1:

$$ I_{\mathcal{H}} = \int_{\mathbb{R}} \pi(d\lambda). \qquad (53) $$

Since $D$ is continuous, we can define its adjoint $D^* : L^2(\Omega, \mu) \to \mathcal{H}_D$, and hence the self-adjoint operator $D^*D : \mathcal{H}_D \to \mathcal{H}_D$. Applying Theorem 2, we obtain a projection-valued measure $\pi_D$ on $\mathbb{R}$ endowed with its Borel $\sigma$-algebra,
such that

$$D^*D = \int_{\mathbb{R}_+} \lambda \pi_D(d\lambda), \qquad I_{\mathcal{H}_D} = \int_{\mathbb{R}_+} \pi_D(d\lambda), \qquad (54)$$

where the integration is restricted to $\mathbb{R}_+$ since $D^*D$ is a positive operator. We can then define

$$\Pi_D^\alpha := \int_{\alpha^2}^{+\infty} \pi_D(d\lambda), \qquad (55)$$

which is an orthogonal projection in $\mathcal{H}_D$ since $\pi_D$ is a projection-valued measure. We then define the regularized space $\mathcal{H}_D^\alpha$ at level $\alpha > 0$ by

$$\mathcal{H}_D^\alpha := \mathrm{Im} \Pi_D^\alpha \subset \mathcal{H}_D. \qquad (56)$$

For any $u \in \mathcal{H}_D^\alpha$, we compute

$$\begin{aligned} \|D[u]\|_{L^2(\Omega, \mu)}^2 &= \langle D[u], D[u] \rangle_{L^2(\Omega, \mu)} = \langle u, D^*D[u] \rangle_{\mathcal{H}_D} \\ &= \left\langle u, \int_{\mathbb{R}_+} \lambda \pi_D(d\lambda) u \right\rangle_{\mathcal{H}_D} \\ &\stackrel{u \in \mathcal{H}_D^\alpha}{=} \left\langle u, \int_{\mathbb{R}_+} \lambda_1 \pi_D(d\lambda_1) \int_{\alpha^2}^{+\infty} \pi_D(d\lambda_2) u \right\rangle_{\mathcal{H}_D} \\ &\stackrel{\pi_D \text{ PVM}}{=} \left\langle u, \int_{\alpha^2}^{+\infty} \lambda \pi_D(d\lambda) u \right\rangle_{\mathcal{H}_D} \\ &= \int_{\alpha^2}^{+\infty} \lambda \langle u, \pi_D(d\lambda) u \rangle_{\mathcal{H}_D} \\ &\geq \alpha^2 \int_{\alpha^2}^{+\infty} \langle u, \pi_D(d\lambda) u \rangle_{\mathcal{H}_D} \stackrel{u \in \mathcal{H}_D^\alpha}{=} \alpha^2 \langle u, u \rangle_{\mathcal{H}_D} = \alpha^2 \|u\|_{\mathcal{H}_D}^2. \end{aligned} \qquad (57)$$

That is,

$$\|D[u]\|_{L^2(\Omega, \mu)} \geq \alpha \|u\|_{\mathcal{H}_D}, \qquad (58)$$

so that Equation (27) is verified on $\mathcal{H}_D^\alpha$. We denote

$$D^\alpha := D|_{\mathcal{H}_D^\alpha}, \qquad (59)$$

the restriction of $D$ to the domain $\mathcal{H}_D^\alpha$. We can now list the fundamental properties:

1. We have

$$\bigcap_{\alpha>0} (\mathcal{H}_D \setminus \mathcal{H}_D^\alpha) = \mathrm{Ker} D \qquad (60)$$

that is, $\lim_{\alpha \to 0} \mathcal{H}_D \setminus \mathcal{H}_D^\alpha = \mathrm{Ker} D$, since for all $\alpha > \beta$, $\mathcal{H}_D^\alpha \subset \mathcal{H}_D^\beta$ by Property 3 of Definition 1. Moreover, by continuity of $D$, we also have $\lim_{\alpha \to 0} D[\mathcal{H}_D^\alpha] = \mathrm{Im} D$. Finally, for each $\alpha > 0$, $D^\alpha$ coincides with $D$ on $\mathcal{H}_D^\alpha$ by construction.

2. As established by Equation (27), $D^\alpha$ is injective and continuous. When restricted to its image, it is therefore bijective and bicontinuous, hence invertible. In particular $D[\mathcal{H}_D^\alpha]$ is closed in $L^2(\Omega, \mu)$, since it is the inverse image of a closed set under $(D^\alpha)^{-1}$. Therefore least-squares solution is well-defined.

3. The least-squares solution of $D^\alpha = f$ is now modified as follows: one projects $f$ onto $\mathrm{Im} D^\alpha$. In particular, if for some $\alpha > 0$ we already have $f \in \mathrm{Im} D^\alpha$, then the regularized least-squares formulation recovers an *exact solution* to the problem $D[u] = f$. This stands in sharp contrast with Ridge regression, where such exact recovery can only occur *asymptotically* in the limit $\alpha \to 0$.

## G.4 Connection to Green's Function

To further highlight the difference between the two regularization schemes, we now reinterpret them through the lens of Green's functions of the operator $D$. Schwencke & Furtlehner (2025, Theorem 2) established in the finite-dimensional case a connection between the natural gradient for PINNs and Green's functions. Their proof relies on Schwencke & Furtlehner (2025, Proposition 3), which will be our starting point. We restate the relevant definitions and results for completeness.
**Definition 2** (Schwencke & Furtlehner, 2025, Definition 9: generalized Green's function). Let $\mathcal{H}$ be an Hilbert space, $D: \mathcal{H} \to L^2(\Omega, \mu)$ be a linear differential operator, $\mathcal{H}_0 \subset \mathcal{H}$ a subspace isometrically embedded in $\mathcal{H}$ and $f \in L^2(\Omega, \mu)$. A generalized Green's function of $D$ on $\mathcal{H}_0$ is then any kernel function $g: \Omega \times \Omega \to \mathbb{R}$ such that the operator:

$$R_{\mathcal{H}_0}: \begin{cases} L^2(\Omega \to \mathbb{R}, \mu) & \to \mathcal{H} \\ f & \to (x \in \Omega \mapsto \int_{\Omega} g(x, s)f(s)\mu(ds)), \end{cases}$$

verifies the equation:

$$D \circ R_{\mathcal{H}_0} = \Pi_{D[\mathcal{H}_0]}^{\perp} \qquad (61)$$

**Proposition 2** (Schwencke & Furtlehner, 2025, Proposition 3). Let $D: \mathcal{H} \to L^2(\Omega, \mu)$ be a linear differential operator, and $\mathcal{H}_0 := \text{Span}(u_p: 1 \le p \le P) \subset \mathcal{H}$ a subspace isometrically embedded in $\mathcal{H}$. Then the generalized Green's function of $D$ on $\mathcal{H}_0$ is given by: for all $x, y \in \Omega$

$$g_{\mathcal{H}_0}(x, y) := \sum_{1 \le p, q \le P} u_p(x) G_{p,q}^{\dagger} D[u_q](y), \qquad (62)$$

with: for all $1 \le p, q \le P$,

$$G_{p,q} := \langle D[u_p], D[u_q] \rangle_{L^2(\Omega \to \mathbb{R}, \mu)}. \qquad (63)$$

**Our goal.** We aim to

(i) generalize Schwencke & Furtlehner (2025, Proposition 3) to arbitrary Reproducing Kernel Hilbert Spaces;

(ii) establish a direct connection to the regularization framework introduced earlier. This will provide a novel reinterpretation of the Green's function in the regularized operator setting.

**Operator framework.** Consider the operator $D: \mathcal{H}_D \to L^2(\Omega, \mu)$ from Section G.1, and assume that there exists an RKHS $\mathcal{H}_0$ isometrically embedded in $\mathcal{H}_D$ (for instance, any finite-dimensional RKHS, see Schwencke & Furtlehner, 2025, Corollary 1). For Schwencke & Furtlehner (2025, Definition 9) to be well-posed, the range $D[\mathcal{H}_0]$ must be a closed subspace of $L^2(\Omega, \mu)$. As argued earlier, this is guaranteed if $D$ is continuously invertible: indeed, in this case

$$D[\mathcal{H}_0] = (D^{-1})^{-1}[\mathcal{H}_0], \qquad (64)$$

and the inverse image of a closed subspace under a continuous operator is closed.

**Key observation.** Thus, to generalize Schwencke & Furtlehner (2025, Proposition 3), we require $D$ to be continuously invertible. Conveniently, this is precisely the property enforced by the regularization schemes we introduced earlier.

In what follows, we first focus on the cutoff regularization, which offers the clearest interpretation in terms of Green's functions. We then briefly revisit the case of Ridge regression. Before delving further into our main goal, let us first establish two general facts.

**Lemma 1.** Let $(\mathcal{H}_0, \|\cdot\|_{\mathcal{H}_0})$ be an RKHS on a set $X$ with reproducing kernel $k$. Suppose that $\|\cdot\|_{\text{bis}}$ is a norm equivalent to $\|\cdot\|_{\mathcal{H}_0}$. Then $(\mathcal{H}_0, \|\cdot\|_{\text{bis}})$ is also an RKHS.

*Proof.* The key point is to show that there exists a reproducing kernel for the inner product $\langle \cdot, \cdot \rangle_{\text{bis}}$ associated with $\|\cdot\|_{\text{bis}}$. Our argument follows the simple reasoning in Paulsen & Raghupathi (2016, Definitions 1–2).

Since, for every $x \in X$, the point evaluation functional

$$\delta_x: u \in \mathcal{H}_0 \mapsto u(x) \qquad (65)$$

is continuous with respect to $\|\cdot\|_{\mathcal{H}_0}$ by the definition of an RKHS, it is also continuous with respect to the equivalent norm $\|\cdot\|_{\text{bis}}$. Therefore, by the Riesz representation theorem, for each $x \in X$, there exists a unique element $k_x^{\text{bis}} \in \mathcal{H}_0$ such that for all $u \in \mathcal{H}_0$

$$\langle k_x^{\text{bis}}, u \rangle_{\text{bis}} = u(x). \qquad (66)$$

In particular, this defines a reproducing kernel for the norm $\|\cdot\|_{\text{bis}}$, given by

$$k_{\text{bis}}(x, y) = \langle k_x^{\text{bis}}, k_y^{\text{bis}} \rangle_{\text{bis}} = k_x^{\text{bis}}(y), \quad \forall x, y \in X. \qquad (67)$$

Hence $(\mathcal{H}_0, \|\cdot\|_{\text{bis}})$ is indeed an RKHS. $\square$
**Lemma 2.** Let $\mathcal{H}_A, \mathcal{H}_B$ be two Hilbert spaces. If $U: \mathcal{H}_A \to \mathcal{H}_B$ is an isometry, then

$$ U^*U = I_{\mathcal{H}_A}, \qquad UU^* = \Pi_{\mathrm{Im} U}. \qquad (68) $$

In particular $\mathrm{Im} U$ is closed in $\mathcal{H}_B$.

*Proof.* The first identity follows from the fact that for all $x, y \in \mathcal{H}_A$,

$$ \langle x, U^*U[y] \rangle_{\mathcal{H}_A} = \langle U[x], U[y] \rangle_{\mathcal{H}_B} = \langle x, y \rangle_{\mathcal{H}_A}. \qquad (69) $$

Thus $(U^*U(y) - y) \in \mathcal{H}_A^\perp$, i.e. $U^*U = I_{\mathcal{H}_A}$. For the second, the key point is to show that $\mathrm{Im} U$ is closed, i.e. $\mathrm{Im} U = \overline{\mathrm{Im} U}$.

Let $y \in \overline{\mathrm{Im} U}$, and $(y_n) \in \mathrm{Im} U^\mathbb{N}$ with $y_n \to y$. Since $(y_n)$ is Cauchy, and $y_n = U(x_n)$ for some $(x_n) \in \mathcal{H}_A^\perp$, we have

$$ \|U(x_n) - U(x_m)\|_{\mathcal{H}_B} = \|x_n - x_m\|_{\mathcal{H}_A}, \qquad (70) $$

so $(x_n)$ is also Cauchy and converges to $x \in \mathcal{H}_A$, since $\mathcal{H}_A$ is complete. Since $U$ is an isometry, we have for all $x \in \mathcal{H}_A$

$$ \|U(x)\|_{\mathcal{H}_B} = \|x\|_{\mathcal{H}_A}. \qquad (71) $$

In particular, $U$ is bounded with operator norm $\|U\| = 1$, and hence continuous. Thus $U(x) = y$, hence $y \in \mathrm{Im} U$. We conclude that $\mathrm{Im} U$ is closed in $\mathcal{H}_B$. Finally:

* For $y \in \mathrm{Im} U$, say $y = U(x)$, we have

$$ UU^*(y) = U(U^*U)(x) = U(x) = y. \qquad (72) $$

* For $y \in (\mathrm{Im} U)^\perp$, we check that $UU^*(y) = 0$. Indeed, for any $z \in \mathcal{H}_B$,

$$ \langle z, UU^*(y) \rangle_{\mathcal{H}_B} = \langle UU^*(z), y \rangle_{\mathcal{H}_B} = 0, \qquad (73) $$

since $UU^*(z) \in \mathrm{Im} U$. Thus $UU^*(y) \in \mathcal{H}_B^\perp$, i.e. $UU^*(y) = 0$. $\square$

We are interested in the restriction of $D$ to the domain $\mathcal{H}_0$. Since the restriction $D^*D: \mathcal{H}_D \to \mathcal{H}_D$ does not, a priori, map $\mathcal{H}_0$ into itself, we first need to adapt the setting in order to apply the spectral theorem of Theorem 2.

Because $\mathcal{H}_0 \subset \mathcal{H}_D$ isometrically, we have for all $u, v \in \mathcal{H}_0$:

$$ \begin{aligned} \langle D[u], D[v] \rangle_{L^2(\Omega, \mu)} &= \langle D[\Pi_{\mathcal{H}_0} u], D[\Pi_{\mathcal{H}_0} v] \rangle_{L^2(\Omega, \mu)} \\ &= \langle \Pi_{\mathcal{H}_0} u, D^*D[\Pi_{\mathcal{H}_0} v] \rangle_{\mathcal{H}_D} \\ &= \langle u, (\Pi_{\mathcal{H}_0} D^*D\Pi_{\mathcal{H}_0})[v] \rangle_{\mathcal{H}_D}, \end{aligned} \qquad (74) $$

where we used in the last step that $\Pi_{\mathcal{H}_0}$ is self-adjoint.

We can therefore apply the spectral theorem Theorem 2 to the bounded self-adjoint operator $\Pi_{\mathcal{H}_0} D^*D\Pi_{\mathcal{H}_0}: \mathcal{H}_0 \to \mathcal{H}_0$, obtaining the analogue of the decomposition in Equation (54):

$$ \Pi_{\mathcal{H}_0} D^* D \Pi_{\mathcal{H}_0} = \int_{\mathbb{R}_+} \lambda \pi_D^{\mathcal{H}_0}(d\lambda), \qquad I_{\mathcal{H}_0} = \int_{\mathbb{R}_+} \pi_D^{\mathcal{H}_0}(d\lambda). \qquad (75) $$

**Regularized spaces.** Fixing $\alpha > 0$, and analogously to Equations (55) and (56), we define the regularized projection and subspace:

$$ \Pi_{D, \mathcal{H}_0}^{\alpha} := \int_{\alpha^2}^{+\infty} \pi_D^{\mathcal{H}_0}(d\lambda), \qquad \mathcal{H}_{D, \mathcal{H}_0}^{\alpha} := \mathrm{Im} \Pi_{D, \mathcal{H}_0}^{\alpha} \subset \mathcal{H}_0 \subset \mathcal{H}_D. \qquad (76) $$

Let $k: \Omega \times \Omega \to \mathbb{R}$ be the reproducing kernel of $\mathcal{H}_0$. Then, by Paulsen & Raghupathi (2016, Theorem 2.5), $\mathcal{H}_{D, \mathcal{H}_0}^\alpha$ is an RKHS with reproducing kernel

$$ k_{\alpha}(x, y) := \Pi_{D, \mathcal{H}_0}^{\alpha}[k(x, \cdot)](y), \quad \forall x, y \in \Omega. \qquad (77) $$
**Norm equivalence.** Since $\mathcal{H}_{D, \mathcal{H}_0}^{\alpha} \subset \mathcal{H}_0 \subset \mathcal{H}_D$, inequality in Equation (26) remains valid, i.e. for all $u \in \mathcal{H}_{D, \mathcal{H}_0}^{\alpha}$:

$$ \|D[u]\|_{L^2(\Omega, \mu)} \le \|u\|_{\mathcal{H}_D}. \qquad (78) $$

Furthermore, by an argument entirely analogous to Equation (57), we also have

$$ \|D[u]\|_{L^2(\Omega, \mu)} \ge \alpha \|u\|_{\mathcal{H}_D}, \quad \forall u \in \mathcal{H}_{D, \mathcal{H}_0}^{\alpha}. \qquad (79) $$

In particular, the functional

$$ \|\cdot\|_D : \left\{ \begin{array}{ll} \mathcal{H}_{D, \mathcal{H}_0}^{\alpha} & \rightarrow \mathbb{R} \\ u & \mapsto \|D[u]\|_{L^2(\Omega, \mu)} \end{array} \right. \qquad (80) $$

defines a norm equivalent to $\|\cdot\|_{\mathcal{H}_D}$ on $\mathcal{H}_{D, \mathcal{H}_0}^{\alpha}$. By Lemma 1, the pair $(\mathcal{H}_{D, \mathcal{H}_0}^{\alpha}, \|\cdot\|_D)$ is itself an RKHS with a reproducing kernel $k_D$.

**Isometry property.** The crucial observation is that $D$ is an isometry with respect to this norm. Indeed, for all $u, v \in \mathcal{H}_{D, \mathcal{H}_0}^{\alpha}$,

$$ \langle u, v \rangle_D = \langle D[u], D[v] \rangle_{L^2(\Omega, \mu)}. \qquad (81) $$

This allows us to characterize the associated Green's function.

**Theorem 1.** The generalized Green's function of the operator $D$ in the regularized space $\mathcal{H}_{D, \mathcal{H}_0}^{\alpha}$ is given, for all $x, y \in \Omega$, by

$$ g_D(x, y) := D[k_D(x, \cdot)](y), \qquad (12) $$

*Proof.* For all $f \in L^2(\Omega, \mu)$ and $x \in \Omega$,

$$ \begin{aligned} \int_{\Omega} g_D(x, s)f(s)\mu(ds) &= \langle g_D(x, \cdot), f \rangle_{L^2(\Omega, \mu)} \\ &= \langle D[k_D(x, \cdot)], f \rangle_{L^2(\Omega, \mu)} \\ &= \langle k_D(x, \cdot), D^*f \rangle_D \\ &= (D^*f)(x). \end{aligned} \qquad (82) $$

Since $D$ is an isometry, Lemma 2 gives $DD^* = \Pi_D[\mathcal{H}_{D, \mathcal{H}_0}^{\alpha}]$. Therefore,

$$ D[x \mapsto \int_{\Omega} g_D(x, s)f(s)\mu(ds)] = D[D^*f] = \Pi_D[\mathcal{H}_{D, \mathcal{H}_0}^{\alpha}]f, \qquad (83) $$

which precisely shows that $g_D$ is a generalized Green's function. $\square$

The key insight of Theorem 1 is that, in the PINNs setting—and most notably in our algorithm—we implicitly construct the reproducing kernel $k_D$ associated with the norm $\|\cdot\|_D$ on the regularized tangent space $T_{\theta}^{\alpha} \mathcal{M}$ of the neural network manifold $\mathcal{M}$, at cutoff level $\alpha$. This kernel is precisely the PINNs NNTK introduced by Schwencke & Furtlehner (2025).

A crucial consequence is that the regularization of the Gram matrix is not merely a “numerical trick” to guarantee stability: it is the very mechanism that ensures the Green’s function is well defined.

**Conceptual interpretation.** This perspective also offers a profound interpretation of the procedure: rather than attempting to invert the operator $D$ directly, we build a kernel $k_D$ whose associated metric makes $D$ an isometry, and thus ensures that $D^*$ acts as the generalized left-inverse of $D$. The magic of the kernel lies in the following facts:

(i) We never need to compute $D^*$ explicitly, since it is implicitly encoded in the relation

$$ \langle D[k_D(x, \cdot)], f \rangle_{L^2(\Omega, \mu)} = \langle k_D(x, \cdot), D^*f \rangle_D. \qquad (84) $$

(ii) The same formula allows us to directly evaluate the generalized solution $D^*f$: indeed, for all $x \in \Omega$, the reproducing property gives

$$ D^*f(x) = \langle k_D(x, \cdot), D^*f \rangle_D. \qquad (85) $$

**Comparison with Ridge regression.** An analogous analysis holds for Ridge regression. However, instead of inverting $D$ “via isometry,” we invert the augmented operator $(D, \alpha I_{\mathcal{H}_D})$.
H Proofs

H.1 Statement and proof of Proposition 3

We start by recalling the following statements from Schwencke & Furtlehner (2025).

**Definition** (Schwencke & Furtlehner, 2025, Definition 4). A linear operator $A: \mathcal{H} \to \mathcal{H}$ is an integral operator given that there is $k: \Omega \times \Omega \to \mathbb{K}, \mathbb{K} \in \{\mathbb{R}, \mathbb{C}\}$, such that: for all $f \in \mathcal{H}$, for all $x \in \Omega$

$$A(f)(x) = \langle k(x, \cdot), f \rangle_{\mathcal{H}}. \qquad (86)$$

**Lemma** (Schwencke & Furtlehner, 2025, Lemma 1). Let us be $\mathcal{H}_0 := \text{Span}(u_p: 1 \le p \le P) \subset \mathcal{H}$ and consider the Gram matrix $G_{pq} := \langle u_p, u_q \rangle_{\mathcal{H}}$ of $(u_p)$ and its eigen-decomposition $G = U\Delta^2U^t$. Then:

$$L_p := \sum_{1 \le q \le P} u_q U_{q,p} \Delta_p^\dagger, \qquad (87)$$

is an orthonormal basis of $\mathcal{H}_0$. In particular, $\Pi_{\mathcal{H}_0}$ is an integral operator whose kernel is:

$$k(x, y) = \sum_{1 \le p, q \le P} u_p(x) G_{p,q}^{\dagger} u_q(y). \qquad (88)$$

Furthermore $L_p$ are the left-singular vector of the so-called synthesis operator<fnref n="3" />:

$$\mathcal{T} := \left\{ \begin{array}{ll} \mathbb{R}^P & \to \mathcal{H}_0 \\ \alpha & \mapsto \sum_{1 \le p \le P} \alpha_p u_p \end{array} \right. . \qquad (89)$$

**Proposition 3.** Given the scalar loss

$$\ell(\boldsymbol{\theta}) := \mathcal{L}(u_{\boldsymbol{\theta}}) \stackrel{(6)}{=} \frac{1}{2} \|u_{\boldsymbol{\theta}} - f\|_{L^2(\Omega, \mu)}^2, \qquad (90)$$

the Natural Gradient update of Equation (8)

$$u_{\boldsymbol{\theta}_{t+1}} \leftarrow u_{\boldsymbol{\theta}_t} - \eta \Pi_{\boldsymbol{\theta}_t \mathcal{M}} (\nabla \mathcal{L}_{u_{\boldsymbol{\theta}_t}}); \qquad \boldsymbol{\theta}_{t+1} \leftarrow \boldsymbol{\theta}_t - \eta \mathrm{d}u_{\boldsymbol{\theta}_t}^{\dagger} (\Pi_{\boldsymbol{\theta}_t \mathcal{M}} (\nabla \mathcal{L}_{u_{\boldsymbol{\theta}_t}})) \qquad (8)$$

can be equivalently written as

$$\boldsymbol{\theta}_{t+1} \leftarrow \boldsymbol{\theta}_t - \eta G_{\boldsymbol{\theta}_t}^{\dagger} \nabla \ell(\boldsymbol{\theta}_t); \qquad G_{\boldsymbol{\theta}_t, p, q} := \langle \partial_p u_{\boldsymbol{\theta}_t}, \partial_q u_{\boldsymbol{\theta}_t} \rangle_{L^2(\Omega, \mu)}. \qquad (9)$$

*Proof.* Since the tangent space $T_{\boldsymbol{\theta}}\mathcal{M}$ of Equation (7):

$$T_{\boldsymbol{\theta}}\mathcal{M} := \text{Im}(du_{\boldsymbol{\theta}}) = \text{Span}(\partial_p u_{\boldsymbol{\theta}}: 1 \le p \le P) \subset \mathcal{H}, \qquad (7)$$

is finite-dimensional, we may invoke Schwencke & Furtlehner (2025, Lemma 1). This result shows that the Natural Neural Tangent Kernel (NNTK), given by

$$NNTK_{\boldsymbol{\theta}}(x, y) := \sum_{1 \le p, q \le P} (\partial_p u_{\boldsymbol{\theta}}(x)) G_{\boldsymbol{\theta}pq}^{\dagger} (\partial_q u_{\boldsymbol{\theta}}(y))^t, \qquad G_{\boldsymbol{\theta}p,q} := \langle \partial_p u_{\boldsymbol{\theta}}, \partial_q u_{\boldsymbol{\theta}} \rangle_{\mathcal{H}}, \qquad (91)$$

is the kernel of the orthogonal projection $\Pi_{\Pi_{\boldsymbol{\theta}}\mathcal{M}}^{\perp}$ onto $T_{\boldsymbol{\theta}}\mathcal{M}$. Therefore, by Equation (86), for all $x \in \Omega$,

$$\Pi_{\Pi_{\boldsymbol{\theta}}\mathcal{M}}^{\perp}(\nabla \mathcal{L}_{|u_{\boldsymbol{\theta}}}) (x) = \langle NNTK_{\boldsymbol{\theta}}(x, \cdot), \nabla \mathcal{L}_{|u_{\boldsymbol{\theta}}} \rangle_{\mathcal{H}} \\ \stackrel{(91)}{=} \sum_{1 \le p, q \le P} \partial_p u_{\boldsymbol{\theta}}(x) G_{\boldsymbol{\theta}pq}^{\dagger} \langle \partial_q u_{\boldsymbol{\theta}}, \nabla \mathcal{L}_{|u_{\boldsymbol{\theta}}} \rangle_{\mathcal{H}}. \qquad (92)$$

Next, note that

$$\langle \partial_q u_{\boldsymbol{\theta}}, \nabla \mathcal{L}_{|u_{\boldsymbol{\theta}}} \rangle_{\mathcal{H}} = d\mathcal{L}_{|u_{\boldsymbol{\theta}}} (\partial_q u_{\boldsymbol{\theta}}) \stackrel{\text{chain rule}}{=} \partial_q \mathcal{L}(u_{\boldsymbol{\theta}}) \stackrel{(90)}{=} \partial_q \ell(\boldsymbol{\theta}). \qquad (93)$$

Therefore, by linearity of $du_{\boldsymbol{\theta}}^{\dagger}$,

$$du_{\boldsymbol{\theta}}^{\dagger} (\Pi_{\Pi_{\boldsymbol{\theta}}\mathcal{M}}^{\perp} (\nabla \mathcal{L}_{|u_{\boldsymbol{\theta}}})) \stackrel{(92),(93)}{=} \sum_{1 \le p, q \le P} du_{\boldsymbol{\theta}}^{\dagger} (\partial_p u_{\boldsymbol{\theta}}) G_{\boldsymbol{\theta}pq}^{\dagger} \partial_q \ell(\boldsymbol{\theta}). \qquad (94)$$

<footnote n="3">³Name and notation are taken from Adcock & Huybrechts (2019).</footnote>
Finally, observe that $\partial_p u_\theta = du_\theta(e^{(p)}),$ where $e^{(p)}$ is the p-th canonical basis vector of $\mathbb{R}^p$. If $du_\theta$ were invertible, we would directly obtain

$$du_{\theta}^{\dagger}(\partial_{p} u_{\theta}) = e^{(p)}, \qquad (95)$$

which would complete the argument. However, this invertibility does not hold in general.

To address this, recall that $du_\theta$ can be identified with the synthesis operator $\mathcal{T}$ introduced in Equation (89) of Schwencke & Furtlehner (2025, Lemma 1). From the final part of that lemma, we know that $\mathrm{Im} \, du_\theta^\dagger = \mathrm{Im} \, G_\theta^\dagger$. Consequently,

$$G_{\theta}^{\dagger} e^{(p)} = G_{\theta}^{\dagger} du_{\theta}(\partial_{p} u_{\theta}). \qquad (96)$$

Putting all pieces together yields the desired update rule, thereby completing the proof. □

## H.2 Ridge-regression implementation ANaGRAM

In the following, we show that a Ridge-regression can be implemented in ANaGRAM's update rule given by Equation (10).

**Proposition 4.** A Ridge-regression can be implemented in the SVD-based update Equation (10) by replacing the pseudo-inverse $\hat{\Delta}^\dagger$ with

$$\left( \frac{\hat{\Delta}_{t,i}}{\hat{\Delta}_{t,i}^2 + S\alpha} \right)_{1 \le i \le r_{\mathrm{svd}}} . \qquad (97)$$

*Proof.* As shown in (Schwencke & Furtlehner, 2025, Section E), the ANaGRAM's update of Equation (10):

$$\theta_{t+1} \leftarrow \theta_t - \eta \hat{\phi}^{\dagger} \nabla \widehat{\mathcal{L}}_{\theta_t}; \qquad \hat{\phi}_{i,p} := \partial_p u_{\theta}(x_i); \qquad (\widehat{\nabla \mathcal{L}}_{\theta})_i := u_{\theta}(x_i) - f(x_i), \qquad (10)$$

is equivalent to the update with the empirical matrix $\hat{G}_\theta$:

$$\theta_{t+1} \leftarrow \theta_t - \eta \hat{G}_{\theta_t}^{\dagger} \nabla \ell(\theta_t); \qquad \hat{G}_{\theta_t} := \frac{1}{S} \hat{\phi}_{\theta_t}^t \hat{\phi}_{\theta_t}, \qquad (98)$$

where $\ell$ is defined in Equation (5):

$$\ell(\theta) := \frac{1}{2S} \sum_{i=1}^{S} (u_{\theta}(x_i) - f(x_i))^2. \qquad (5)$$

Thus, we get immediately

$$\nabla \ell(\theta_t) = \frac{1}{S} \hat{\phi}^t \nabla \widehat{\mathcal{L}}_{\theta} = \frac{1}{S} \hat{U} \hat{\Delta} \hat{V}_{\theta}^t \nabla \widehat{\mathcal{L}}_{\theta}, \qquad (99)$$

where we used the SVD decomposition of $\hat{\phi}$:

Using Equation (100) again, we have

$$\hat{G}_{\theta} = \frac{1}{S} \hat{U} \hat{\Delta}_{\theta}^{2} \hat{U}_{\theta}^{t}, \qquad (101)$$

thus for all $\alpha > 0$

$$\hat{G}_{\theta} + \alpha I_d = \frac{1}{S} \hat{U} \hat{\Delta}_{\theta}^2 \hat{U}_{\theta}^t + \alpha \hat{U} \hat{U}^t = \hat{U} \left( \mathrm{diag} \left( \frac{\hat{\Delta}_{\theta_i}^2}{S} + \alpha \right)_{1 \le i \le r_{\mathrm{svd}}} \right) \hat{U}_{\theta}^t, \qquad (102)$$

which implies

$$(\hat{G}_{\theta} + \alpha I_d)^{-1} = \hat{U} \left( \mathrm{diag} \left( \frac{S}{\hat{\Delta}_{\theta_i}^2 + S\alpha} \right)_{1 \le i \le r_{\mathrm{svd}}} \right) \hat{U}_{\theta}^t. \qquad (103)$$

This finally yields

$$(\hat{G}_{\theta} + \alpha I_d)^{-1} \nabla \ell(\theta_t) \stackrel{(99)}{=} \hat{U} \left( \mathrm{diag} \left( \frac{S}{\hat{\Delta}_{\theta_i}^2 + S\alpha} \right)_{1 \le i \le r_{\mathrm{svd}}} \right) \hat{U}_{\theta}^t \frac{1}{S} \hat{U} \hat{\Delta} \hat{V}_{\theta}^t \nabla \widehat{\mathcal{L}}_{\theta} \qquad (104)$$

$$= \hat{U} \left( \mathrm{diag} \left( \frac{\hat{\Delta}_{t,i}}{\hat{\Delta}_{\theta_i}^2 + S\alpha} \right)_{1 \le i \le r_{\mathrm{svd}}} \right) \hat{V}_{\theta}^t \nabla \widehat{\mathcal{L}}_{\theta}, \qquad (105)$$

which concludes the proof. □
## H.3 Proof of Proposition 1

To prove Proposition 1, we need the following lemma:

**Lemma 3.** For $1 \le M \le N \le r_{svd}$:

$$ (RCE_M^S)^2 - (RCE_N^S)^2 = \frac{1}{S} \left\| \Pi_N^M \hat{V}^T \hat{\nabla} \mathcal{L} \right\|_{\mathbb{R}^S}^2. \quad (106) $$

*Proof.* Let us first recall the definition of the $RCE_N^S$ in Equation (13), namely

$$ RCE_N^S := \frac{1}{\sqrt{S}} \left\| \hat{V} \Pi_N^0 \hat{V}^T \hat{\nabla} \mathcal{L} - \hat{\nabla} \mathcal{L} \right\|_{\mathbb{R}^S}. \quad (13) $$

Fixing $1 \le N \le M \le r_{svd}$ and applying the same reasoning as in Equation (116) to $RCE_M^S$ and $RCE_N^S$ (see the proof of Proposition 1 in Section H.3), we get

$$ S(RCE_M^S)^2 = \widehat{\nabla}_{\mathcal{L}_{\theta}}^t \widehat{\nabla}_{\mathcal{L}_{\theta}} - \sum_{p=1}^{M} (\hat{V}_{\theta_p}^t \widehat{\nabla}_{\mathcal{L}_{\theta}})^2; \quad S(RCE_N^S)^2 = \widehat{\nabla}_{\mathcal{L}_{\theta}}^t \widehat{\nabla}_{\mathcal{L}_{\theta}} - \sum_{p=1}^{N} (\hat{V}_{\theta_p}^t \widehat{\nabla}_{\mathcal{L}_{\theta}})^2, \quad (107) $$

and therefore

$$ \begin{aligned} S\left((RCE_N^S)^2 - (RCE_M^S)^2\right) &= \sum_{p=1}^{N} (\hat{V}_{\theta_p}^t \widehat{\nabla}_{\mathcal{L}_{\theta}})^2 - \sum_{p=1}^{M} (\hat{V}_{\theta_p}^t \widehat{\nabla}_{\mathcal{L}_{\theta}})^2 \\ &\stackrel{M \le N}{=} \sum_{p=M+1}^{N} (\hat{V}_{\theta_p}^t \widehat{\nabla}_{\mathcal{L}_{\theta}})^2 = \sum_{p=M+1}^{N} (e^{(p)^T} \hat{V}^T \widehat{\nabla}_{\mathcal{L}_{\theta}})^2 \\ &= \sum_{p=M+1}^{N} (\widehat{\nabla}_{\mathcal{L}_{\theta}}^t \hat{V} e^{(p)}) (e^{(p)^T} \hat{V}^T \widehat{\nabla}_{\mathcal{L}_{\theta}}) \\ &= \widehat{\nabla}_{\mathcal{L}_{\theta}}^t \hat{V} \left( \sum_{p=M+1}^{N} e^{(p)^T} e^{(p)^T} \right) \hat{V}^T \widehat{\nabla}_{\mathcal{L}_{\theta}} \quad (108) \\ &= \langle \hat{V}^T \widehat{\nabla}_{\mathcal{L}_{\theta}}, \Pi_N^M \hat{V}^T \widehat{\nabla}_{\mathcal{L}_{\theta}} \rangle_{\mathbb{R}^S} \\ &\stackrel{\Pi_N^M = \Pi_N}{=} \langle \Pi_N^M \hat{V}^T \widehat{\nabla}_{\mathcal{L}_{\theta}}, \Pi_N^M \hat{V}^T \widehat{\nabla}_{\mathcal{L}_{\theta}} \rangle_{\mathbb{R}^S} \\ &= \left\| \Pi_N^M \hat{V}^T \widehat{\nabla}_{\mathcal{L}_{\theta}} \right\|_{\mathbb{R}^S}^2, \end{aligned} $$

where we use in the penultimate equality, the fact that $\Pi_N^M$ is an orthogonal projection. $\square$

**Remark 5.** The above lemma provides an interesting property that gives a further understanding of what is happening during the flattening, i.e. $RCE_M^S - RCE_N^S \approx 0$. In particular, as $(RCE_M^S)^2 - (RCE_N^S)^2 = (RCE_M^S - RCE_N^S)(RCE_M^S + RCE_N^S)$, therefore flattening for the components in the range $[N_{flat}, r_{cutoff}]$ means that $\frac{1}{S} \left\| \Pi_N^M \hat{V}^T \widehat{\nabla} \mathcal{L} \right\|_{\mathbb{R}^S}^2 \approx 0$. In other words, the problem is "learned" for those components, as the projection of the functional gradient (which is propotional to the error) on their corresponding span is null. The proof of this lemma is provided in Section H.3.

**Proposition 1.** $RCE_N^S$ is a non-increasing function of $N$, i.e. for all $1 \le M, N \le r_{svd}$:

$$ M \le N \implies RCE_M^S \ge RCE_N^S. \quad (15) $$

Furthermore, assuming that $(x_i)_{i=1}^S$ are i.i.d sampled from $\mu$, we have $\mu$-almost surely

$$ \lim_{S \to \infty} RCE_N^S = \left\| \nabla_{\mathcal{L}_{u\theta}} - \Pi_{T_N^0 M}^\perp \nabla_{\mathcal{L}_{u\theta}} \right\|_{L^2(\Omega, \mu)} = \left\| \Pi_{[T_N^0 M]^\perp} \nabla_{\mathcal{L}_{u\theta}} \right\|_{L^2(\Omega, \mu)}, \quad (16) $$

where $T_N^M M = \text{Span}(V_{t,i} : M \le i \le N)$, while $(V_{t,i})_{1 \le i \le r_{svd}}$ are the right singular-vectors of the differential $du_\theta$ ordered in a decreasing order according to their associated singular values.
Proof. The first statement is a direct consequence of Lemma 3 proven above.

Let us now show that the second statement takes place. Since $\nabla \mathcal{L}_{u_{\theta}} \in L^2(\Omega, \mu)$ and $\mathrm{Im} \, du_{\theta} \subset L^2(\Omega, \mu)$, the law of large numbers yields that for all $1 \le p, q \le P$

$$ \lim_{S \to \infty} \frac{1}{S} \sum_{i=1}^{S} [\nabla \mathcal{L}_{u_{\theta}}(x_i)]^2 = \lim_{S \to \infty} \frac{1}{S} \widehat{\nabla \mathcal{L}_{\theta}^t \nabla \mathcal{L}_{\theta}} = \int_{\Omega} [\nabla \mathcal{L}_{u_{\theta}}(x)]^2 \mu(dx) \quad \text{a.s.} \quad (109) $$

$$ \lim_{S \to \infty} \frac{1}{S} \sum_{i=1}^{S} \partial_p u_{\theta}(x_i) \nabla \mathcal{L}_{u_{\theta}}(x_i) = \lim_{S \to \infty} \frac{1}{S} \hat{\phi}_{\theta_p}^t \widehat{\nabla \mathcal{L}_{\theta}} = \int_{\Omega} \partial_p u_{\theta}(x) \nabla \mathcal{L}_{u_{\theta}}(x) \mu(dx) \quad \text{a.s.} \quad (110) $$

$$ \lim_{S \to \infty} \frac{1}{S} \sum_{i=1}^{S} \partial_p u_{\theta}(x_i) \partial_q u_{\theta}(x_i) = \lim_{S \to \infty} \frac{1}{S} \hat{\phi}_{\theta_p}^t \hat{\phi}_{\theta_q} = \int_{\Omega} \partial_p u_{\theta}(x) \partial_q u_{\theta}(x) \mu(dx) \quad \text{a.s.} \quad (111) $$

In particular, this implies

$$ \lim_{S \to \infty} \frac{1}{S} \hat{\phi}^t \hat{\phi} = G_{\theta} = U_{\theta} \Delta_{\theta}^2 U_{\theta}^t \quad \text{a.s.} \quad (112) $$

Since the eigenvectors (and eigenvalues) are continuous functions of the matrix coefficients (by polynomial dependence through the characteristic polynomial) and taking into account that $\frac{1}{S}\hat{\phi}^t\hat{\phi} = \frac{1}{S}\hat{U}\Delta_{\theta}^2\hat{U}^t$, this yields

$$ \lim_{S \to \infty} \hat{U} = U_{\theta} \quad \text{a.s.}; \qquad \lim_{S \to \infty} \frac{1}{S} \hat{\Delta}^2 = \Delta_{\theta}^2 \quad \text{a.s.} \quad (113) $$

By continuity of the square root and the inverse on $\mathbb{R}^*_+$, we get that for all $1 \le p \le P$ such that $\Delta_{\theta_p} > 0$

$$ \lim_{S \to \infty} \sqrt{S} \hat{\Delta}_{\theta_p}^{-1} = \Delta_{\theta_p}^{-1} \quad \text{a.s.} \quad (114) $$

and thus for all $1 \le p \le P$ such that $\Delta_{\theta_p} > 0$, we have almost surely

$$ \begin{aligned} \lim_{S \to \infty} \frac{1}{\sqrt{S}} \hat{V}_{\theta_p}^T \widehat{\nabla \mathcal{L}_{\theta}} &= \lim_{S \to \infty} \sqrt{S} \hat{\Delta}_{\theta_p}^{-1} \hat{U}_{\theta_p}^t \left( \sum_{q=1}^{P} e^{(q)} e^{(q)^t} \right) \frac{1}{S} \hat{\phi}^t \widehat{\nabla \mathcal{L}_{\theta}} \\ &= \sum_{q=1}^{P} \left( \lim_{S \to \infty} \sqrt{S} \hat{\Delta}_{\theta_p}^{-1} \right) \left( \lim_{S \to \infty} \hat{U}_{\theta_p}^t e^{(q)} \right) \left( \lim_{S \to \infty} \frac{1}{S} \hat{\phi}_{\theta_q}^t \widehat{\nabla \mathcal{L}_{\theta}} \right) \\ &= \sum_{q=1}^{P} \Delta_{\theta_p}^{-1} U_{\theta_p}^t e^{(q)} \int_{\Omega} \partial_q u_{\theta}(x) \nabla \mathcal{L}_{u_{\theta}}(x) \mu(dx) \\ &= \int_{\Omega} du_{\theta} (U_{\theta_p} \Delta_{\theta_p}^{-1}) (x) \nabla \mathcal{L}_{u_{\theta}}(x) \mu(dx) \\ &= \int_{\Omega} V_{\theta_p}(x) \nabla \mathcal{L}_{u_{\theta}}(x) \mu(dx), \end{aligned} \quad (115) $$

where we used in the last equality, the identification of the singular vectors of $du_{\theta}$ in (Schwencke & Furtlehner, 2025, Lemma 1 p. 24, section C.2). Returning to the definition of the $\mathrm{RCE}_N^S$ in Equation (13), namely

$$ \mathrm{RCE}_N^S := \frac{1}{\sqrt{S}} \left\| \hat{V} \Pi_N^0 \hat{V}^T \widehat{\nabla \mathcal{L}} - \widehat{\nabla \mathcal{L}} \right\|_{\mathbb{R}^S}, \quad (13) $$
we get

$$
\begin{align*}
S (\mathrm{RCE}_N^S)^2 &= \langle \hat{V} \Pi_N^0 \hat{V}^t \hat{\nabla} \mathcal{L}_\theta - \hat{\nabla} \mathcal{L}_\theta, \hat{V} \Pi_N^0 \hat{V}^t \hat{\nabla} \mathcal{L}_\theta - \hat{\nabla} \mathcal{L}_\theta \rangle_{\mathbb{R}^S} \\
&= \hat{\nabla} \mathcal{L}_\theta^t \hat{\nabla} \mathcal{L}_\theta + \hat{\nabla} \mathcal{L}_\theta^t \hat{V} \underbrace{\Pi_N^0 \hat{V}^t \hat{V}}_{=I_d} \Pi_N^0 \hat{V}^t \hat{\nabla} \mathcal{L}_\theta - 2 \hat{\nabla} \mathcal{L}_\theta^t \hat{V} \Pi_N^0 \hat{V}^t \hat{\nabla} \mathcal{L}_\theta \\
&= \hat{\nabla} \mathcal{L}_\theta^t \hat{\nabla} \mathcal{L}_\theta - \hat{\nabla} \mathcal{L}_\theta^t \hat{V} \Pi_N^0 \hat{V}^t \hat{\nabla} \mathcal{L}_\theta \\
&= \hat{\nabla} \mathcal{L}_\theta^t \hat{\nabla} \mathcal{L}_\theta - \hat{\nabla} \mathcal{L}_\theta^t \hat{V} \left( \sum_{p=1}^N e^{(p)} e^{(p)^t} \right) \hat{V}^t \hat{\nabla} \mathcal{L}_\theta \\
&= \hat{\nabla} \mathcal{L}_\theta^t \hat{\nabla} \mathcal{L}_\theta - \sum_{p=1}^N (\hat{V}_{\theta_p}^t \hat{\nabla} \mathcal{L}_\theta)^2,
\end{align*}
$$

where in the second equality, we use the fact that $\hat{V}$ is orthogonal and $\Pi_N^0$ is a projection. Combining Equations (109) and (115), this yields

$$
\lim_{S \to \infty} (\mathrm{RCE}_N^S)^2 = \int_{\Omega} \nabla \mathcal{L}_{u_\theta}(x)^2 \mu(dx) - \sum_{p=1}^{N} \left( \int_{\Omega} V_{\theta_p}(x) \nabla \mathcal{L}_{u_\theta}(x) \mu(dx) \right)^2 \quad \text{a.s.} \quad (117)
$$

By Fubini's theorem, we have almost surely

$$
\begin{align*}
\sum_{p=1}^{N} \left( \int_{\Omega} V_{\theta_p}(x) \nabla \mathcal{L}_{u_{\theta}}(x) \mu(dx) \right)^2 &= \int_{\Omega^2} \nabla \mathcal{L}_{u_{\theta}}(x) \left( \sum_{p=1}^{N} V_{\theta_p}(x) V_{\theta_p}(y) \right) \nabla \mathcal{L}_{u_{\theta}}(y) \mu(dx) \mu(dy) \\
&= \int_{\Omega} \nabla \mathcal{L}_{u_{\theta}}(x) \Pi_{\text{Span}(V_{\theta_i}: 1 \le i \le N)}^{\perp} \nabla \mathcal{L}_{u_{\theta}}(x) \mu(dx) \tag{118} \\
&= \left\| \Pi_{\text{Span}(V_{\theta_i}: 1 \le i \le N)}^{\perp} \nabla \mathcal{L}_{u_{\theta}} \right\|_{L^2(\Omega, \mu)}^2,
\end{align*}
$$

where in the second equality, we used (Schwencke & Furtlehner, 2025, Theorem 4 p. 23, section C.2) and the fact that
$(\Pi_{\text{Span}(V_{\theta_i}: 1 \le i \le N)}^{\perp})^2 = \Pi_{\text{Span}(V_{\theta_i}: 1 \le i \le N)}^{\perp}$ in the third. Therefore, from Equation (117) and Equation (118)

$$
\lim_{S \to \infty} (\mathrm{RCE}_N^S)^2 = \| \nabla \mathcal{L}_{u_\theta} \|_{L^2(\Omega, \mu)}^2 - \| \Pi_{\mathrm{Span}(V_{\theta_i}: 1 \le i \le N)}^\perp \nabla \mathcal{L}_{u_\theta} \|_{L^2(\Omega, \mu)}^2 \quad \text{a.s.,} \quad (119)
$$

$$
= \| \nabla \mathcal{L}_{u_\theta} - \Pi_{\text{Span}(V_{\theta_i}: 1 \le i \le N)}^{\perp} \nabla \mathcal{L}_{u_\theta} \|_{L^2(\Omega, \mu)}^2 \quad \text{a.s.,} \quad (120)
$$

where in the second equality, we use in the reverse order a reasoning similar to Equation (116). Finally, we obtain

$$
\| \nabla \mathcal{L}_{u_\theta} - \Pi_{\text{Span}(V_{\theta_i}: 1 \le i \le N)}^{\perp} \nabla \mathcal{L}_{u_\theta} \|_{L^2(\Omega, \mu)}^2 = \| \Pi_{\text{Span}(V_{\theta_i}: 1 \le i \le N)}^{\perp} \nabla \mathcal{L}_{u_\theta} \|_{L^2(\Omega, \mu)}^2 \quad (121)
$$

which comes from the canonical decomposition in Hilbert spaces, i.e. using that Span($V_{\theta_i}$ : $1 \le i \le N$) is a closed
subspace and

$$
\nabla \mathcal{L}_{u_{\theta}} = \Pi_{\text{Span}(V_{\theta_i}: 1 \le i \le N)}^{\perp} \nabla \mathcal{L}_{u_{\theta}} + \Pi_{\text{Span}(V_{\theta_i}: 1 \le i \le N)}^{\perp} \nabla \mathcal{L}_{u_{\theta}}. \quad (122)
$$

This completes the proof.

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

**Corollary 1.** For $1 \le M \le N \le r_{svd}$:

$$
\lim_{S \to \infty} (RCE_M^S)^2 - (RCE_N^S)^2 = \| \Pi_{T_N^M}^\perp \nabla \mathcal{L}_{u_\theta} \|_{L^2(\Omega)}^2 \quad (123)
$$

*Proof.* Apply Proposition 1 to Equation (106) of Lemma 3.

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