# Learning and extrapolating scale-invariant processes

Anaclara Alvez-Canepa, Cyril Furtlehner, François Landes

► To cite this version:

Anaclara Alvez-Canepa, Cyril Furtlehner, François Landes. Learning and extrapolating scale-invariant processes. 2025. hal-05470825

HAL Id: hal-05470825

https://inria.hal.science/hal-05470825v1

Preprint submitted on 22 Jan 2026

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Learning and extrapolating scale-invariant processes

Anaclara Alvez$^{1,2}$, Cyril Furtlehner$^{1,3}$, and François Landes$^{1,2}$

$^{1}$Laboratoire Interdisciplinaire des Sciences du Numérique (LISN)

$^{2}$Université Paris-Saclay and

$^{3}$INRIA-Saclay

(Dated: January 22, 2026)

Machine Learning (ML) has deeply changed some fields recently, like Language and Vision. In the case of complex systems, spectacular breakthroughs happened too (e.g. for protein folding) and more are expected to come. Our question is: how and to which extent can one regress scale-free processes, i.e. processes displaying power law behavior, like earthquakes or avalanches? The events one is interested in predicting are the large ones, i.e. events that are typically rare in the training set, so we are basically in the extrapolation regime. While some recent works also tackle scale-free systems by proposing generative models closely aligned with the renormalization group framework, here instead we explore the problem of prediction in the extrapolation regime. We consider two paradigmatic problems that are statistically self-similar. The first one is a 2-dimensional fractional Gaussian field obeying linear dynamics, self-similar by construction and amenable to exact analysis. The second one is the Abelian sandpile model, exhibiting self-organized criticality.

The emerging paradigm of Geometric Deep Learning shows that including known symmetries into the model's architecture is key to success (as translation invariance for images). Here one may hope to extrapolate only by leveraging scale invariance, which is however a peculiar symmetry, as it involves possibly non-trivial coarse-graining operations and anomalous scaling. We perform experiments on various existing architectures like U-net, Riesz network (scale invariant by construction), or our own proposals: a wavelet-decomposition based Graph Neural Network (with discrete scale symmetry), a Fourier embedding layer and a Fourier-Mellin Neural Operator. Based on these experiments and a complete characterization of the linear case, we identify the main issues relative to spectral biases and coarse-grained representations, and discuss how to alleviate them with the relevant inductive biases.

# I. INTRODUCTION

In the last decade machine learning (ML) and more specifically deep neural networks (DNN) have thoroughly renewed the research perspectives in many fields like Natural Language Processing and Computer Vision. Despite indisputable successes, the introduction of ML approaches to describe or predict complex physical systems remains a challenge [1, 2], driven by the specificity of physics data. Unlike natural image or language datasets, physics data are governed by symmetries, conservation laws, causal relations, and complex systems frequently involves rare, highly non-linear events that are critical to their dynamics. These characteristics impose specific statistical constraints that traditional, general purpose ML models are not inherently designed to handle, requiring the introduction of innovative approaches to ensure both accuracy and interpretability. By complex physics we mean systems displaying qualitatively similar phenomena over a broad range of scales, where all or many scales are relevant and thus cannot be decoupled. Setting aside multifractality—a hallmark of complexity—we observe that even monofractal systems, which by definition exhibit self-similarity under scale transformations, already pose significant challenges to standard machine learning models.

The purpose of this work is to discuss two important aspects of this problem. The first aspect concerns spectral bias [3], that most ML models suffer from, and which becomes a critical issue in this context. Generically

most ML models learn only the low frequencies, as these dominate the loss, making them unable or very slow to learn also the high frequency features. This problem will appear as a crucial point in the examples we present. The second aspect deals with the inability of standard ML models to extrapolate. While conventional machine learning models excel at interpolation within the training distribution, they struggle to extrapolate beyond it or into its tails. For complex systems the rare events are usually the extreme ones, that play a crucial role in the physical process: they cannot be neglected.

Formally we can define the ability of an ML model to extrapolate as follows. Given some physical process in the form of a function $f_{GT}(\varphi) = \psi$ (deterministic or not) connecting a pair $(\varphi, \psi)$ of input-output quantities (not specified here, but think for instance of a field configuration before and after an avalanche), and a corresponding real observable $s = g(\varphi, \psi) \in \mathbb{R}$ (e.g. the size of the avalanche), we want to be able to train a neural network $f_{\theta}$ on a subset $S_{\text{train}} \subset \{(\varphi, \psi), \mathbb{1}_{\{g(\varphi, \psi)<s^{\text{th}}\}}\}$ corresponding to observation below some threshold $s^{\text{th}} \in \mathbb{R}$ and still be able to generalize well to any pair $(\varphi, \psi)$, in particular contained in the subset $\mathbb{1}_{\{g(\varphi, \psi)\ge s^{\text{th}}\}}$ which are typically those which are too rare to be observed during training. Intuitively this seems possible if the true process $f_{GT}$ obeys a form of scale-symmetry, such that large-scale events ($s > s^{\text{th}}$) are similar and share the same properties as the small scales ones, up to some well adapted rescaling. The question of how to inform the ML model with the appropriate scale symmetry is the
subject of the present work.

The practical importance of symmetries in Neural Networks has been recognized since the beginning, with the early introduction by Y. Le Cun of convolutional networks[4, 5], that ensure (approximate) translation symmetry. Another more recent example is Graph Neural Networks [6], that are built so as to ensure node permutation equivariance, since node labeling should not impact modeling or predictions. For general groups, the idea of building a neural network that is exactly equivariant by construction, has then been formalized in a series of papers [7]. Let us recall the definition of equivariance. An operator (e.g. a Neural Network) $f: X \to Y$ is equivariant with respect to a group $G$ if and only if:

$$ \forall g \in G: \quad \rho_Y(g) \circ f = f \circ \rho_X(g) \quad (1) $$

where $\rho_X, \rho_Y$ are the representations of the group $G$ in the vector spaces $X, Y$. Intuitively, equivariance means that if the input is transformed according to $g$, the output also transforms according to the same transformation $g$. Invariance is a special case of equivariance where one chooses a trivial representation $\forall g \in G, \rho_Y(g) = I$, thus restricting the expressivity of the operator $f$.

Depending on the complexity of the symmetry group, the number of generators can be very large and possibly infinite. Since standard implementations are based on the regular representation, i.e. one latent variable per generator, the size of the equivariant layer can be potentially infinite, hence a tradeoff has to be found between the computational or memory burdens and the faithfulness of the representations (projection of the input onto the basis corresponding to the generators). For the rotation group, useful in many scientific applications, there is a large body of work, from the early theoretical papers [8] up to the diffusion of a dedicated torch package [9]. In this steerable network approach [8], the point is to have representations built in the appropriate basis, i.e. that corresponding to the irreducible representation of the group of interest, and perform products in the proper space, so as to maintain equivariance.

In the context of Computer Vision, the notions of scale invariance or equivariance have recently been introduced [10, 11]. In appendix C we recall these definitions and show to which extent they apply to statistically self-similar processes. It is worth noticing that many works focus on the problem of correctly classifying images regardless of the size of the object of interest in the picture, i.e. deal with an invariant task. In this sense, they deal with a markedly different problem from the one we want to address. Spatial transformers [12] is a general purpose method to perform arbitrary affine transformation of arbitrary input feature map to facilitate classification for instance, but is not as such able to generalize to unseen scales. Some other methods are based on multi-scale channels where either the input [13] or the (steerable) filters are replicated at different scales [11]. Some methods use continuous filters to generate multiscale representation of the input [14, 15] like Gaussian derivative

filters [16] or Riesz ones [17], the latter yielding scale invariant networks by construction. Fourier input is also considered in [18] in order to correct for antialiasing when input feature maps are downscaled with continuous filters.

Among these methods only a few of them may readily generalize to unseen scales, notably the Riesz networks [17] can be successfully used to segment cracks of all sizes on tomographic images of concrete, even when trained solely on a subset of crack sizes. A key element to this success is of course that cracks are self-similar objects: the shapes of small cracks is distributed similarly to that of larger cracks.

Independently from this line of work, for complex systems, connections between the renormalization group (RG) flow and neural networks have emerged as a compelling alternative approach to conventional supervised learning models. While such methods are typically limited in scope, notable exceptions exist – such as proposals to construct generative models directly from RG equations [19, 20]. Beside image processing, the development of scale invariant/equivariant architectures in the context of analyzing and forecasting critical systems with non-trivial scaling behavior appears to us as a virgin ground.

Our context is distinct from the usual omputer vision setup, where the output, e.g. a classification label, is invariant w.r.t. to the symmetry group. In our context we instead seek an equivariance in a statistical sense where Eq. (1) is now replaced by

$$ \forall g \in G: \quad \rho_Y(g) \circ f \stackrel{d}{=} f \circ \rho_X(g) \quad (2) $$

meaning that a change of the input by a group element does change the output distribution. This stems from the requirement that for a self-similar field, scaling operations—such as zooming in or out—necessitate generating missing components of the data and accounting for anomalous scaling behavior. Specifically, when the input is rescaled, either its fine (small-scale) and coarse (large-scale) components, as well as those of the output, must be synthesized to preserve the invariance of their joint distribution. The simplest example is zooming in on a Brownian motion: to 'look closer,' one must sample small scale Brownian paths, as the process is inherently stochastic. For the type of data we study (e.g., fractional Gaussian fields), scaling introduces non-trivial, exponent-dependent behavior (e.g., tied to the Hurst exponent) that must explicitly be accounted for. Our goal is to design a neural network that inherently incorporates these scaling laws.

We investigate these questions under the following functional regression framework where we have a pair of input-output fields $(\varphi, \psi)$ both defined on a finite square lattice $\Omega \subset \mathbb{Z}^2$ of size $|\Omega| = L^2$. The fields are related by

$$ \psi(\mathbf{x}) = F[\varphi](\mathbf{x}) + \epsilon(\mathbf{x}) \quad (3) $$

where $F$ is an unknown operator, $\epsilon$ is a noise field also defined on the square lattice, and all $\varphi, \psi$ and $\epsilon$ belong to some family of scale-free fields with long range
spatial correlations and self-similarity characterized by critical exponents. F and the various fields are assumed to be such that equation (3) is equivariant w.r.t. scale symmetry, equation (2). Given a certain number of observations $\{(\varphi^{(s)}, \psi^{(s)}), s = 1, \dots, M\}$, the goal is then to predict $\psi$ from $\varphi$. Note that this regression problem has recently been framed using the neural operator approach in the literature [21–24].

Concretely we consider examples with different levels of difficulty, out of which lessons can be learned. The first self-similar toy-problem of functional regression we consider involves Fractional Gaussian fields (FGF) [25] subject to a linear functional map $F$, where the extrapolation task is defined by training on data with truncated spectrum and testing on data with full spectrum. In the cases of primary interest—specifically, when masking either the high-frequency or low-frequency ends of the spectrum—this corresponds, in direct space, to super-resolution or extrapolation to larger scales, respectively. In practice, this means testing on examples at scales not encountered during training.

Our second example concerns non-linear regression in the context of self-organized criticality: the sandpile model [26, 27] is an avalanche process paradigmatic of scale free-processes, where the size distribution of events $S$ follows a $1/S$ scaling. The extrapolation is then defined as the ability to predict large avalanches when training solely on small and intermediate ones.

Equipped with these two quite different test cases, we investigate how to inform the network with the scale symmetry, in order to first fit the training set in its full (broad) spectrum, and possibly extrapolate to an extra piece of the spectrum, under the assumption of a self-similar task. Inspired by the mathematical framework of steerable networks that some of us worked with [28, 29], we are looking to work in the appropriate bases, such that operations are by construction scale-invariant or equivariant.

## II. DATASETS AND TASKS:

In this section we present the two physical systems to be studied – the Fractional Gaussian field (FGF) and the Abelian Sandpile Model (ASM) – and specify the tasks we want to solve. The FGF serves as a relatively basic, well-understood example of a scale invariant system, to which we add simple dynamics. The ASM is a toy model for self-organized criticality and avalanche phenomena, more complex due to the fact that, despite having avalanche processes with long-range correlations (large susceptibility), the spatial correlations decay rapidly.

## A. Fractional Gaussian field

We first consider a simple multi-dimensional Fractional Gaussian field (FGF). Formally it is defined in the continuum in arbitrary dimension as [25]

$$ \varphi = (-\Delta)^{-\beta/4} W $$

where $\Delta$ is the Laplace operator, $W$ a white noise on $\mathbb{R}^d$ and $\beta \in \mathbb{R}^{+*}$ the fractional exponent. It is fully characterized by the covariance structure of its Fourier modes given by

$$ \mathbb{E}[\tilde{\varphi}(\mathbf{k})\tilde{\varphi}(\mathbf{k}')] \propto S(\mathbf{k})\delta(\mathbf{k} - \mathbf{k}') $$

with what is usually called power spectrum:

$$ S(\mathbf{k}) \propto |\mathbf{k}|^{-\beta}. \qquad (4) $$

Formally its Boltzmann distribution reads

$$ P(\varphi) = \frac{1}{Z} e^{-\frac{1}{2} \int d^2 \mathbf{k} E(\mathbf{k})} \qquad (5) $$

with the energy given by

$$ E(\mathbf{k}) \stackrel{\text{def}}{=} k^{\beta} |\tilde{\varphi}(\mathbf{k})|^2 \qquad (6) $$

With this energy, the Boltzmann distribution is invariant under rescaling. Note that the standard Gaussian free field corresponds to $\beta = 2$. The FGF is also commonly parametrized by the Hurst exponent $H = \frac{\beta-d}{2} \in [0, 1]$ to quantify the roughness and self-similarity of the field. The real space covariance structure in this range of parameters is then given by

$$ \mathbb{E}[\varphi(\mathbf{x})\varphi(\mathbf{x}')] \propto |\mathbf{x} - \mathbf{x}'|^{2H} $$

and is typically long-range, tunable by $H$. According to this scaling, we see that the FGF displays self-similarity, in the sense that

$$ \varphi(s\mathbf{x}) \stackrel{d}{=} s^{\frac{\beta-d}{2}} \varphi(\mathbf{x}), \qquad (7) $$

$$ \tilde{\varphi}(\mathbf{k}/s) \stackrel{d}{=} s^{\frac{\beta+d}{2}} \tilde{\varphi}(\mathbf{k}). \qquad (8) $$

These equations can be seen as the renormalization equation for the FGF, which involves no other exponent than the anomalous dimension of the field.

To be concrete, in the following we will discretize our FGF $\varphi: \Omega \subset \mathbb{Z}^2 \to \mathbb{C}$ on an $L \times L$ lattice, which amounts to generate a collection of random Fourier coefficients $\tilde{\varphi}(\mathbf{k})$:

$$ \varphi(\mathbf{x}) = \frac{1}{|\Omega|} \sum_{\mathbf{k}} e^{i2\pi \mathbf{k} \cdot \mathbf{x}} \tilde{\varphi}(\mathbf{k}), \qquad (9) $$

with $\mathbf{k} = (k_x, k_y)$ parameterized by $L^2$ pairs of in integers as $\mathbf{k}_{i,j} = (-\frac{1}{2} + \frac{i}{L}, -\frac{1}{2} + \frac{j}{L})$ with $i, j = 0, \dots, L-1$. Each
coefficient’s real and imaginary part is chosen to follow a normal distribution:

$$Re[\tilde{\varphi}(\mathbf{k})] \stackrel{d}{=} Im[\tilde{\varphi}(\mathbf{k})] \stackrel{d}{=} \mathcal{N}(0, |\mathbf{k}|^{-\beta}), \quad (10)$$

except for the choice $\tilde{\varphi}(0) = 0$, so that field average is zero by construction: $\sum_{\mathbf{x}} \varphi(\mathbf{x}) = 0$.

Using this field, we now consider two different simple linear processes preserving the FGF distribution, hence displaying a form of equivariance. The first one corresponds to a phase-mixing dynamics and the second one to a local spectral flow associated to scale transform and rotations:

$$\frac{\partial \tilde{\varphi}(\mathbf{k}, t)}{\partial t} = -i\nu(\mathbf{k})\tilde{\varphi}(\mathbf{k}, t) \quad (\text{Phase Mixing}) \quad (11)$$

$$\frac{\partial \tilde{\varphi}(\mathbf{k}, t)}{\partial t} = i(\eta T_{\text{scale}}(\beta) + \omega T_{\text{rot}}) \tilde{\varphi}(\mathbf{k}, t) \quad (\text{Spectral flow}) \quad (12)$$

where $\nu(\mathbf{k})$ is an arbitrary phase function, while $T_{\text{scale}}$ and $T_{\text{rot}}$ are two (commuting) self-adjoint operators

$$T_{\text{scale}}(\beta) \stackrel{\text{def}}{=} -i \left( \mathbf{k} \cdot \nabla_{\mathbf{k}} + 1 + \frac{\beta}{2} \right) \quad (13)$$

$$T_{\text{rot}} \stackrel{\text{def}}{=} i(k_x \partial_{k_y} - k_y \partial_{k_x}) \quad (14)$$

corresponding respectively to scale and rotation transformations with speed controlled respectively by parameters $\eta, \omega \in \mathbb{R}$. These two processes represent the only local linear and stationary transformations that we can think of that preserve the FGF distribution; we consider them separately because they do not commute and want to see how to treat them specifically. The corresponding operators (in particular by definition of $T_{\text{scale}}(\beta)$) are indeed self-adjoint with respect to the inner product

$$\langle \varphi_1, \varphi_2 \rangle_{\beta} = \int d^2 \mathbf{k} k^{\beta} \varphi_1(\mathbf{k}) \varphi_2^*(\mathbf{k})$$

induced by the norm defining the energy (6). The solutions to these processes are found in closed form,

$$\tilde{\varphi}(\mathbf{k}, t) = \tilde{\varphi}(\mathbf{k}, 0)e^{-i\nu(\mathbf{k})t}, \quad (\text{Phase Mixing})$$

$$\tilde{\varphi}(\mathbf{k}, t) = \tilde{\varphi}(\mathcal{R}_{\omega t} \mathbf{k} s(t), 0) s(t)^{\beta/2+1}, \quad (\text{Spectral flow})$$

with $\mathcal{R}_{\alpha \mathbf{k}} = (k_x \cos(\alpha) + k_y \sin(\alpha), -k_x \sin(\alpha) + k_y \cos(\alpha))$ and $s(t) = e^{\eta t}$. This means the evolution is simply a change of phase at every $\mathbf{k}$ with respect to the initial configuration $\tilde{\varphi}(\mathbf{k}, 0)$ in the first case, while in the second there is a combined radial and orthoradial current $J(\mathbf{k}, t) = k^{\beta}|\phi(\mathbf{k}, t)|^2(\eta \mathbf{k} - \omega \mathbf{k}^{\perp})$ of energy flow spiraling toward large [resp. small] scales when $\eta > 0$ [resp. $\eta < 0$]. The inference task then consists in regressing $\psi(\mathbf{k}) = \varphi(\mathbf{k}, 1)$ from $\varphi(\mathbf{k}, 0)$. The phase mixing illustrates the case of a purely local transformation (in $\mathbf{k}$ space) while for the spectral flow we have an example of a non-local transformation (at least in $\mathbf{k}$ space).

<figure><img src="image_2.png" /><figcaption>FIG. 1. Sample input and output with lattice size $L^2 = 64 \times 64$ for the FGF spectral flow inference task. The transformation applied corresponds to a rescaling of factor $s = 0.9$ and a rotation of $\pi/4$.</figcaption></figure>

Both conserve the energy (6). In addition in both cases the scaling properties (7,8) are preserved by the dynamics, since a phase shift doesn't affect the power spectrum (which in turns, controls the scaling property) in the first case while in the second case the energy flow conserves also the energy except possibly at the edges of the domain, which means that some energy has to be injected into the system at small scales to compensate for the energy leaving the system at large scale.

In Fig. 1 we show a pair of input and output samples for the FGF spectral flow process. Note that the corners of the image are distorted due to the rotation which will induce some limitation in the performance of models corresponding to idealistic conditions.

Our goal is to evaluate the model's ability to extrapolate across scales, ensuring that the test set includes scales not present in the training data. To achieve this, we set the Fourier coefficients within a band $k_1 < k < k_2$ to zero for all samples in the training set. At test time, the model is evaluated on samples containing the full range of frequencies. This thereby includes both super-resolution experiment and large scale extrapolation when hiding respectively high and low frequencies. In order to measure the level of extrapolation we define an extrapolation factor for both the super-resolution and large scale extrapolation as follows. If $[1/k_{\text{max}}, 1/k_{\text{min}}]$ is the range of considered scales in the test set and if $1/k_{\text{occ}}$ designate the scale above [resp. below] which the spectrum is occulted in the train set of the large scale extrapolation [resp. super-resolution] experiment we define, for a scaling factor $s < 1$ of the transformation the extrapolation factor as

$$f_{\text{extrapol}} \stackrel{\text{def}}{=} \begin{cases} \frac{sk_{\text{occ}}}{k_{\text{min}}} & (\text{large scale extrapolation}) \\ \frac{sk_{\text{max}}}{k_{\text{occ}}} & (\text{super-resolution}) \end{cases} \quad (15)$$

which corresponds to the ratio between the largest scale $1/k_{\text{min}}$ that can be predicted in the test and largest one $1/s k_{\text{occ}}$ in the train [resp. the smallest in the train $1/k_{\text{occ}}$]
and the smallest one $1/k_{\text{occ}}$ in the train] in the large scale extrapolation [resp. super-resolution] experiment.

## B. Abelian Sandpile Model

The Abelian sandpile model (ASM) is a paradigmatic example of self-organized criticality, a concept introduced by Bak, Tang, and Wisenfeld in 1987 [30] to group together physical systems that, without needing to fine-tune any control parameter, exhibit critical-point like behavior. In practice, this often relates to the dissipation parameter being vanishing. These systems spontaneously evolve from any configuration towards a set of states that display scale invariance and power law distributions of physical quantities [31, 32]. There exists a certain number of variations of the ASM, here we focus on one of them.

The ASM is defined on a square lattice with $N$ sites. To each site $i$ we associate a variable $z_i$ representing the height of the sandpile at each site, and a threshold value $z_{ci}$. The dynamics consists in two steps:

1. Select a site at random with uniform probability, and add a grain of sand to that site

$$z_i \mapsto z_i + 1 \qquad (16)$$

2. If at any site $z_i > z_{ci}$, that site becomes unstable. Topple any unstable sites

$$z_j \mapsto z_j - \Delta_{ij} \quad \text{for all } j = 1, \dots, N \qquad (17)$$

where we take the Laplacian matrix $\Delta$ to be an $N \times N$ matrix containing ones in the positions corresponding to nearest neighbours, and zeros everywhere else.

Repeat until all sites are stable. The sequence of topplings that occur until stability is reached is called an avalanche.

In practice we will take the critical height to be $z_{ci} = 4$. Given this value of the Laplacian, during the avalanche the energy is conserved in the bulk of the lattice and is dissipated at the boundary. This dissipation ensures that avalanches remain finite.

After a possible transient regime, the dynamics reaches a set of typical configurations (recurrent under deterministic driving) and the ASM displays critical behavior, in the form of power law distributions. There are only two independent critical exponents that control these power laws. The first one corresponds to the distribution of avalanche sizes $P(s)$, where $s$ is defined as the number of sites that topple during an avalanche. The power law is written as

$$P(s) \sim s^{-\tau}. \qquad (18)$$

The second one, often called dynamical exponent, is related to the duration of the avalanches, defined as the

number of consecutive updates that need to be done to reach stability. It has the scaling law:

$$t \sim r^z, \qquad (19)$$

where $r$ is the linear size (radius) of the avalanche. Other critical exponents, such as those corresponding to the avalanche duration distribution and linear size distribution, can be written in terms of $\tau$ and $z$ [27, 33].

Unlike the FGF, the ASM does not have an energy function, and it is not clear a priori how the renormalization group can be applied to extract the critical behavior. An approximate spatial renormalization group scheme for the ASM was proposed in the 90's by Vespignani et. al. [34], in which the stochastic process that defines the dynamics is renormalized, rather than the field itself. This scheme yields values of the critical exponents that are in overall agreement with numerical simulations.

In rough terms, the idea of this renormalization scheme is to define coarse-grained variables. Blocks at a given scale are said to be stable if, when energy is transferred to it, no energy is transferred to its neighbors (topplings inside the block at smaller scales may still occur). Similarly, a block is said to be critical if it will topple when energy is transferred to it, and will transfer energy to its neighbours at the same scale. The dynamics can be characterized by the density $\rho$ of critical cells, and the probabilities $p = (p_1, p_2, p_3, p_4)$, where $p_i$ is the probability that energy will be transferred to $i$ neighbours when a cell topples. Imposing certain rules on how to go from one scale to the next one, RG equations can be obtained for $(\rho, p)$. However, the approximations used in this procedure are not fully controlled [27] and some discrepancies were found between the predicted scaling behaviour and what is observed in simulations.

The task we want to solve is the following: starting from the unstable configuration right before an avalanche occurs, we want to predict which sites will topple (the shape of the avalanche). In Fig. 2 we show a sample input-label pair. The input field $\varphi(\mathbf{x})$, shown on the left, is an $L \times L$ lattice of discrete variables corresponding to the height of the sandpile at each site $\mathbf{x}$. There is one single unstable site at height $\varphi(\mathbf{x}_0) = 5$, marking the starting point of the avalanche. The target field $\psi(\mathbf{x})$ is a binary map, where the sites with value $\psi(\mathbf{x}) = 1$ are the ones that toppled during the avalanche that results from the initial configuration $\varphi(\mathbf{x})$. Our model's output $\hat{\psi}(\mathbf{x})$ can be interpreted as the probability of site $\mathbf{x}$ to topple during the avalanche. We can threshold this map to make it binary and obtain a predicted avalanche shape.

Although the process given by eq. (17) to go from the input to the output is deterministic and relatively simple, the rapidly decaying correlations in the input map makes it impossible to guess from the naked eye, and difficult to learn for a neural network.

Extrapolation to large scale avalanches is of particular interest, as it provides a realistic baseline for applications to real systems in which large events are rare. An example is earthquake prediction, for which large-scale events
<figure><img src="image_3.png" /><figcaption>FIG. 2. Sample input and label with lattice size 32 × 32 for the regression task in the ASM. In the input, the unstable site corresponding to the epicenter of the avalanche has a height of 5. The output is a binary map, the ones marking the sites that toppled during the avalanche.</figcaption></figure>

are scarce [35]. In order to test extrapolation in the ASM, the natural setup is to train on samples with a maximum avalanche size $S$, leaving out the largest avalanches, and then test events of size $s > S$. Since there are few large avalanches however, this training is not markedly different from the standard one. Presently we train on all samples and test on all. Note that while in the case of the FGF the input and target fields share the same self-similarity properties (same exponent $\beta$ for input and target fields), here instead the scaling properties differ between input and output, which are different types of fields, with different critical exponents. This greatly increases the difficulty of the extrapolation problem.

### III. MACHINE LEARNING MODELS

Here we describe the neural network architectures we designed to solve the proposed tasks, along with known models that serve as baselines. Considering our problem setting (3), models can be described in the framework of neural operators (NO) [21]:

$$ \hat{\psi}[\varphi](\mathbf{x}) = \sigma\left(\sum_{\mathbf{x}' \in \Omega} F_{\theta}(\mathbf{x}, \mathbf{x}')\varphi(\mathbf{x}') + b_{\theta}(\mathbf{x})\right) \quad (20) $$

where $\sigma$ is an activation function and $F_{\theta}$ and $b_{\theta}$ are respectively a linear operator and a functional bias, to be learned. As we shall see, except from the U-Net, all models considered in the following can be cast in this framework or its iterated version (i.e. many layers of NO (20)). In order to realize extrapolation across scales, we explore different ways of inserting scale invariance in this framework. The Fourier embedding network is a model designed specifically for the FGF problem, that operates in Fourier space, and shares weights across all frequencies. Then we propose the Fourier-Mellin neural operator (FMNO), as a way to obtain scale invariance via convolutions in scales. The Riesz network [17] makes use of scale-equivariant operations too. Finally, we propose a GNN architecture that operates on a wavelet representation of the data, somewhat inspired by the idea of

[19]. Additionally, we use a U-net without built-in scale invariance as a baseline for all tasks.

#### A. U-net

The U-net [36] is a widely used architecture designed originally for segmentation tasks. It follows an encoder-decoder structure: the encoder operates by performing consecutive convolutions and coarse-grainings on the input, then the decoder upsamples the representation obtained to the original input size. The information at each encoding layer is recorded and then used by the decoder at the corresponding level, thus preserving the fine-grain information and mixing it with higher-level features. For a more detailed explanation, see Appendix A 1.

The hierarchical scale structure of this network makes it a reasonable choice to treat the problems we are interested in. However, as the weights learned at different scales are independent from each other, it does not have any type of scale-invariance, so it is not expected to extrapolate to scales unseen during training. In this work we use it as a baseline to compare with the architectures that do have some type of implementation of scale invariance.

#### B. Riesz network

One way of constructing a scale-invariant neural network is to impose the constraint that all operations within it must commute with rescalings. This approach is used in [17] to build the Riesz network, an architecture designed for segmenting pictures of multiscale cracks in concrete. This architecture is based on the Riesz transform, an operator that commutes with continuous scaling transformations.

Given a continuous map $\varphi: \mathbb{R}^2 \to \mathbb{R}$, the Riesz transform $\mathcal{R}_j: L(\mathbb{R}^2) \to L(\mathbb{R}^2)$ is defined as

$$ \mathcal{R}_j[\varphi](\mathbf{x}) = C \lim_{\varepsilon \to 0} \int_{\mathbb{R}^2 \setminus \mathcal{B}_\varepsilon} d\mathbf{x}' \frac{x_j' \varphi(\mathbf{x} - \mathbf{x}')}{|\mathbf{x}'|^3}, \quad (21) $$

where $j = 1, 2$ indicates the 2 spatial directions and $\mathcal{B}_\varepsilon$ is the ball of radius $\varepsilon$ centered at the origin. In Fourier space the Riesz transform takes the simpler form

$$ \tilde{\mathcal{R}}_j[\varphi](\mathbf{k}) = -i \frac{k_j}{|\mathbf{k}|} \tilde{\varphi}(\mathbf{k}). \quad (22) $$

This operator commutes with naive rescaling transformations, i.e.

$$ \mathcal{R}_i \left[ \varphi \left( \frac{\cdot}{s} \right) \right] = \left[ \mathcal{R}_i \varphi \right] \left( \frac{\cdot}{s} \right). \quad (23) $$

For proof of this property, refer to [17]. One can define higher order Riesz transforms by successively applying $\mathcal{R}_j$. In particular we have three second-order Riesz
transform given by

$$ \mathcal{R}^{(m,n)}[\varphi](\mathbf{x}) = \mathcal{R}_1^m \circ \mathcal{R}_2^n[\varphi](\mathbf{x}), \quad m+n=2; m, n \in \mathbb{N}. \quad (24) $$

The basic component of the Riesz network is the Riesz layer. We consider for simplicity linear combinations of the two first-order and the three second-order Riesz transforms

$$ J[\varphi] = C_0 + \sum_{k=1}^{2} C_k \mathcal{R}_k[\varphi] + \sum_{m+n=2} C_{mn} \mathcal{R}^{(m,n)}[\varphi]. \quad (25) $$

The coefficients ($C_0, C_k, C_{mn}$) of the linear combination will be the learnable parameters of the network. Higher order Riesz transforms could in principle be included, but they increase computational cost and are not needed for the scope of this work.

The Riesz layer number $\ell$ will take an input with a number $C(\ell)$ of channels $\varphi^{(\ell)} = (\varphi_1^{(\ell)}, \dots, \varphi_{C(\ell)}^{(\ell)})$, and output a number of $C(\ell + 1)$ of channels via

$$ \varphi_c^{(\ell+1)} = \sum_{c'=1}^{C(\ell)} J_{\ell}^{(c,c')}[\varphi_{c'}^{(\ell)}], \quad (26) $$

meaning that a given channel $c$ at layer $\ell + 1$ is connected to the channel $c'$ at layer $\ell$ via the application and subsequent linear combination of the five Riesz filters considered, and the result is summed over all input channels. $J_{\ell}^{(c,c')}$ is defined as

$$ J_{\ell}^{(c,c')}[\varphi] = C_0 + \sum_{k=1}^{2} C_k^{(\ell,c,c')} \mathcal{R}_k[\varphi] + \sum_{m+n=2} C_{mn}^{(\ell,c,c')} \mathcal{R}^{(m,n)}[\varphi], \quad (27) $$

allowing for different learnable parameters for each connection between an input channel $c'$ and an output channel $c$, and for each layer $\ell$. Layer $\ell$ thus has $6C(\ell)C(\ell+1)$ real parameters to be learned.

The Riesz network is then built by successive applications of Riesz layers alternated with non-linearities. Note that the non-linearities should also preserve the scale-equivariance property – this is the case for rectified linear units activations (ReLU) [17], which we will use in this work. Any function that is piecewise linear in $]-\infty, 0[\text{ and in } ]0, \infty[\text{ fulfills this criterion, which means in practice that we can use only ReLU and leaky ReLU. In practice, the Riesz layer is implemented in Fourier space, so the network is very similar to a Fourier neural operator [22] in the sense that, for instance, for a 2-layer network we have

$$ \hat{\psi}[\varphi](\mathbf{x}) = \sigma(\mathcal{F}^{-1}\tilde{J}_2\mathcal{F}\sigma(\mathcal{F}^{-1}\tilde{J}_1\mathcal{F}(\varphi(\mathbf{x})))) \quad (28) $$

where $\varphi^{(\ell+1)} = \tilde{J}_\ell \varphi^{(\ell)}$ is the operation done in (26), and $\mathcal{F}, \mathcal{F}^{-1}$ are the Fourier transform and its inverse, applied channel-wise. The number of layers and channels is chosen depending on the task, see appendix for details.

### C. Fourier embedding network

In the FGF phase mixing task, the main difficulty lies in the spectral bias. The high frequency components make it hard for standard architectures to perform the regression. However, the time evolution we want to learn is local in Fourier space, i.e. each Fourier coefficient $\varphi(\mathbf{k})$ is multiplied by a factor depending only on its own mode $\mathbf{k}$, not on the other $\varphi(\mathbf{k}')$, $\mathbf{k}' \neq \mathbf{k}$. Taking this into consideration we propose a model that operates in Fourier space, with the form

$$ F_{\theta}[\tilde{\varphi}] = g_{\theta} \odot \tilde{\varphi} \quad (29) $$

where $g_{\theta}(\mathbf{k})$, a function of the frequencies of the lattice, is learned and multiplied element wise with $\tilde{\varphi}$. The ground truth imposed by the dynamics corresponds to

$$ g(\mathbf{k}) = e^{-i\nu(\mathbf{k})}, \quad (30) $$

with the concrete choice

$$ \nu(\mathbf{k}) = \nu|\mathbf{k}|^2 \quad (31) $$

made for the experiments, $\nu$ being a real parameter. In practice we use an MLP for $g_{\theta}(\mathbf{k})$, but rather than applying it directly to the frequency modes $|\mathbf{k}|$, inspired by [37], we first embed these modes themselves in Fourier. The embedding is $f(\mathbf{k}) = (f_1(\mathbf{k}), \dots, f_P(\mathbf{k}))$ with

$$ f_p(\mathbf{k}) = e^{i2\pi|\mathbf{k}|p/P}, \quad (32) $$

where we often set $P = 100$. This embedding is designed to facilitate the training of the MLP. It depends only on the modulus of $\mathbf{k}$, not its direction, so as to impose rotational symmetry. The network then reads:

$$ F_{\theta}(\tilde{\varphi}(\mathbf{k})) = \text{MLP}_{\theta}(f(|\mathbf{k}|)) \odot \tilde{\varphi}(\mathbf{k}), \quad (33) $$

where the MLP's output is a single scalar (complex) value. Note that the same MLP is applied at every value of $\mathbf{k}$, meaning that the weights are shared across all modes. This network cannot perform well on the rescaling nor rotation tasks because it is local in $\mathbf{k}$ space by construction and specialized to the phase shifts task. Hence we do not perform experiments in these directions.

### D. Fourier-Mellin network

In addition to the issue of spectral bias, the spectral flow of the FGF presents an additional problem, since the process is now non-local in Fourier space. As explained in more detail in section V, the appropriate basis for this task is the Fourier-Mellin (FM) basis, in which the kernel that describes the rescaling and rotating processes is diagonal.

We propose a model that operates in the Fourier-Mellin basis:

$$ F_{\theta}[\varphi_{(\lambda, \mu)}^{*}] = \theta_{(\lambda, \mu)} \odot \varphi_{(\lambda, \mu)}^{*}, \quad (34) $$
where $\varphi^*_{(\lambda, \mu)}$ denotes the Fourier-Mellin transform of the field $\varphi^* = \mathcal{F}[\tilde{\varphi}^*]$, where

$$ \tilde{\varphi}^*(u, \alpha) \stackrel{\text{def}}{=} \tilde{\varphi}\left[\mathbf{k} = \left(e^u \cos(\alpha), e^u \sin(\alpha)\right)\right] $$

is the Fourier component as function of $u = \log|\mathbf{k}|$ and $\alpha = \text{angle}(\mathbf{k}/|\mathbf{k}|)$. Then $\lambda$ and $\mu$ denote respectively the Fourier-Mellin radial and angular momentum that are conjugate to $u$ and $\alpha$. Each Fourier-Mellin mode being multiplied by a weight, this is a linear model in FM space. The number of independent parameters is thus the number of components used when projecting on the FM basis.

Thus, the challenge is not on the ML side, but lies in the implementation of the Fourier-Mellin transform for discrete data. As in [38], we do it by first performing a change of basis from cartesian $(\mathbf{k})$ to logpolar coordinates $(u, \alpha) = (\log|\mathbf{k}|, \text{angle}(\mathbf{k}/|\mathbf{k}|))$, followed by the application of a Fourier transform. Our input field $\varphi(\mathbf{x})$ is homogeneous and the origin does not play any special role, apart from being the center of the spectral flow task we consider. On the contrary, the origin of $\tilde{\varphi}(\mathbf{k})$ is not arbitrary, as $\mathbf{k} \sim 0$ corresponds to the low frequencies. The logpolar setting is thus well motivated. This coordinate change implies interpolating between the input's cartesian grid and a logpolar grid. As this interpolation is not invertible, some numerical error is introduced when inverting the transformation.

To minimize the reconstruction error of the coordinate change, we introduce the Dirichlet kernel, which gives the interpolation between the two grids:

$$ \tilde{\varphi}_{\mathbf{k}^*} = \sum_{\mathbf{k} \in \Omega} D(\mathbf{k}^* - \mathbf{k}) \tilde{\varphi}_{\mathbf{k}}, \quad (35) $$

with $\{\mathbf{k}^* \in \Omega^*\}$ denoting the points of the logpolar grid, $\{\mathbf{k} \in \Omega\}$ the points of the cartesian grid, and

$$ D(\mathbf{q}) = D_{1d}(q_x) D_{1d}(q_y), \quad (36) $$

with the 1d Dirichlet kernel $D_{1d}$ defined as

$$ D_{1d}(q) = \begin{cases} \frac{1}{L} \frac{\sin(\pi L q)}{\sin(\pi q)} & q \neq 0 \\ 1 & q = 0 \end{cases} \quad (37) $$

with $L$ the linear size of the Cartesian grid. Note that the field being interpolated must be periodic, which is not the case for the radial component in logpolar coordinates, due to the power spectrum. Therefore, we multiply by the inverse spectrum in order to have a statistically periodic gaussian field, and take the Fourier-Mellin transform of the field with constant variance. We restore the power spectrum in the decoder step, after going back to cartesian coordinates. The details of this method are explained further in the Appendix A.

The Fourier-Mellin transform is then numerically implemented by computing

$$ \varphi^*_{(\lambda, \mu)} = \frac{1}{M} \sum_{\mathbf{k}, \mathbf{k}^* \in \Omega, \Omega^*} D(\mathbf{k}^* - \mathbf{k}) |\mathbf{k}|^{\beta/2} \varphi_{\mathbf{k}} e^{i(2\pi \log |\mathbf{k}^*|\lambda + \alpha\mu)}, \quad (38) $$

where $M$ is the number of Fourier-Mellin modes chosen and $\alpha = \text{angle}(\mathbf{k}^*)$. In particular, we can choose the resolution of the logpolar grid, which can be different from the original grid. The inverse transformation is numerically done by means of the pseudo-inverse $D^\dagger$ of the Dirichlet kernel, using its SVD, so we have

$$ \tilde{\varphi}_{\mathbf{k}} = \frac{1}{|\mathbf{k}|^{\beta/2}} \sum_{\mathbf{k}^* \in \Omega^*} D^{\dagger}(\mathbf{k}^* - \mathbf{k}) \times \left( \frac{1}{M} \sum_{(\lambda, \mu)} \varphi^*_{(\lambda, \mu)} e^{-i(2\pi \log |\mathbf{k}^*|\lambda + \alpha\mu)} \right). \quad (39) $$

The inversion of the kernel is the main source of numerical error in this method. We find empirically that the linear size of the logpolar grid needs to be at least around 3 times larger than that of the cartesian grid in order to have a good reconstruction of the field.

## E. Wavelet-GNN

We propose an original architecture consisting of a Graph Neural Network (GNN) operating on a graph-structured representation of the input in wavelet space. This takes inspiration from the Wavelet Conditional Renormalization Group approach developed in [19], where the wavelet representation is used to build a hierarchical generative model that is applied to physical systems, such as the $\varphi^4$ model and gravitational lensing data.

The wavelet transform is a decomposition of a signal into a basis of functions called wavelets. There are different families of wavelet bases. In this work, we use the Daubechies wavelets, which are functions that are both localized in direct and Fourier space [39], making them well-suited to represent data at different scales. In practice we use the library made available in [40], which allows for backpropagation through the wavelet transform in PyTorch.

A 2-dimensional discrete wavelet transform applied to an input $\varphi_0 = \varphi$ of size $L \times L$ gives a set of wavelet coefficients $\varphi_1$ of size $L/2 \times L/2$, with 3 channels each roughly representing the horizontal, vertical, and diagonal [41] details, i.e. the fast degrees of freedom. The transform also yields a coarse-grained map $\varphi_1$ of the input, of size $L/2 \times L/2$, that contains the slow degrees of freedom. This totals to $(L/2)^2(3+1) = L^2$ coefficients.

The wavelet transform is invertible, so the original map $\varphi_0$ can be perfectly recovered from $(\tilde{\varphi}_1, \varphi_1)$. Applying the wavelet transform again to $\varphi_1$, we can again separate the fast degrees of freedom at the next scale, obtaining $(\tilde{\varphi}_2, \varphi_2)$. Iterating this procedure, we build an explicit scale representation of the input, as illustrated in Fig. 3. Strictly speaking, reconstruction only needs the set of coefficients $(\varphi_{\log_2(L)}, \tilde{\varphi}_{\log_2(L)}, \tilde{\varphi}_{\log_2(L)-1}, \dots, \tilde{\varphi}_1)$, where we assume that $L$ is a power of 2. In practice, it is convenient to also use the coarse-grained maps
($\varphi_{\log_2(L)}, \varphi_{\log_2(L)-1}, \dots, \varphi_1$) as input features.

<figure><img src="image_4.png" /><figcaption>FIG. 3. Wavelet decomposition of the input. $\varphi_0$ can be perfectly recovered from $\varphi_2, \bar{\varphi}_2, \bar{\varphi}_1$. At each level, the number of “pixels” is divided by 4.</figcaption></figure>

For our purposes, given an input sample of size $L \times L$, with $L = 2^J$ a power of 2, we arrange the wavelet coefficients and coarse-grained maps into a graph with a tree structure. There are $2^{J-1} \times 2^{J-1}$ leave nodes in the tree (level $\ell = 1$), and to each such node $i$ is assigned a 4-component feature vector, $(\bar{\varphi}_1(i), \varphi_1(i))$. Subsequent levels of the tree are built in the same way, each level $\ell$ containing $2^{J-\ell} \times 2^{J-\ell}$ nodes and corresponding to a certain scale of the wavelet representation, with the root node (level $\ell = J$) containing the 4 features $(\bar{\varphi}_J, \varphi_J)$.

Directed edges are drawn between nearest neighbors over the current $2^\ell \times 2^\ell$ lattice. Directed edges are also drawn between levels, in such a way that nodes connected together correspond to the same spatial region. In the bulk of the tree, a node at level $\ell$ will have four children nodes at level $\ell - 1$, and one parent node at level $\ell + 1$. We assign a tag to each (directed) edge, indicating the type of connection it represents. There are 9 different tags: one for the parent node, 4 possible nearest neighbours and 4 possible children. Having these 9 possible independent relations enables the model to distinguish spatial orientation and scale direction.

This representation of the data allows to use GNN architectures in which the basic operations are given by a relational graph convolutional operator. We pick one of the simplest instantiations of such a convolution [42]:

$$ \mathbf{h}'_i = \sum_{r \in \mathcal{R}} \frac{1}{|\mathcal{N}_r(i)|} \sum_{j \in \mathcal{N}_r(i)} \Theta_r \cdot \mathbf{h}_j, \quad (40) $$

where $i$ is now the node index (not directly related to spatial location), $\mathbf{h}_i$ is the feature vector at node $i$, $r$ is the tag index, and $\mathcal{N}_r(i)$ is the set of $r$-type neighbours of $i$. The learnable weight matrices $\Theta_r$ are distinct for each tag type but shared across all similar edges.

One convolutional GNN layer consists in the nodes update defined by (40), followed by a non-linearity applied

<figure><img src="image_5.png" /><figcaption>FIG. 4. Graph representation of a $4 \times 4$ input map, where we only show the edges outgoing from two nodes (highlighted in red). Each level in the tree corresponds to a scale of the wavelet decomposition (3 levels here). The 9 edge types are colored differently, in green for the child to parent relation, in shades of blue for the parent-children relations, and shades of orange/yellow for the neighbors relations. The node features $h_i$ contain the three wavelet coefficients and the value of the coarse-grained field at each position at the given scale. In dashed lines we show the edges that would continue if the tree continued downward.</figcaption></figure>

node-wise and feature-wise. This design allows for a simple empirical implementation of scale invariance: by sharing the weights $\Theta_r$ across the levels of the tree, the information transfer is done identically at each scale. In practice, we allow for different weights at the top $\Theta_r^{J-1}$ and at the bottom level $\Theta_r^1$, to account for the boundary effects related to these 2 special scales.

To obtain an output map with the same size as the input, we choose to have 4 channels in the nodes of the last representation (after the last convolution layer), so as to be able to interpret them as wavelet coefficients and coarse-grained features. We then perform an inverse wavelet transform, using the implementation in [40]. We can count the parameters in a layer: the bulk nodes of the tree share 9 matrices $\Theta_r$ of shape $c_{out} \times c_{in}$, where $c_{in}, c_{out}$ are the number of channels in input and output, respectively. The weight matrices at the top and bottom of the tree are independent from the bulk one: $\Theta_r^{L-1}$ has only 1 tag, parent, and $\Theta_r^1$ has 5 tags (neighbors, parent). The total is thus $15c_{out}c_{in}$ parameters.

## F. Loss function

For the FGF, we perform a pixel-wise regression of a complex-valued map and thus use by default the mean-squared error (MSE) loss function:

$$ \mathcal{L}(\psi, \hat{\psi}) \stackrel{\text{def}}{=} \frac{1}{NL^2} \sum_{s=1}^{N} \|\hat{\psi}^{(s)} - \psi^{(s)}\|^2 \quad (41) $$
which can be indifferently written in the direct pixel or in its Fourier representation

$$ \|\hat{\psi}^{(s)} - \psi^{(s)}\|^2 = \sum_{\mathbf{x} \in \Omega} |\hat{\psi}_{\mathbf{x}}^{(s)} - \psi_{\mathbf{x}}^{(s)}|^2 \quad (42) $$

$$ = \sum_{\mathbf{k} \in \Omega} |\hat{\psi}_{\mathbf{k}}^{(s)} - \tilde{\psi}_{\mathbf{k}}^{(s)}|^2. \quad (43) $$

In practice, for the phase mixing task we operate in Fourier space to save a gradient retro-propagation through the fast Fourier transform (FFT).

In the phase mixing task, spectral bias is strong, but our Fourier embedding network is able to tackle it. For the spectral flow task also, the gradient descent over this loss suffers from a strong spectral bias induced by the power law behavior of the spectrum. For the Fourier-Mellin network, in order to avoid that, we can redefine the loss by changing the norm: instead of the $L^2$ norm, we use that induced by the energy function (6). In discrete form, the loss then corresponds to:

$$ \mathcal{L}(\psi, \hat{\psi}) = \frac{1}{NL^2} \sum_{\mathbf{k} \in \Omega} k^{\beta} |\hat{\psi}_{\mathbf{k}}^{(s)} - \tilde{\psi}_{\mathbf{k}}^{(s)}|^2 \quad (44) $$

$$ \approx \frac{1}{NL^2} \sum_{(\lambda, \mu) \in \Omega^*} |\hat{\psi}_{(\lambda, \mu)}^{(s)} - \psi_{(\lambda, \mu)}^{(s)}|^2 \quad (45) $$

where the last equation involves the Fourier-Mellin coefficients, see equation (38). In the continuous limit this last equation should be an equality, but in practice, due to finite size effects and discretization problems discussed in the appendix, the relation is only approximate. Still it is more efficient to use the Fourier-Mellin representation of the loss, again to save the retro-propagation through several FFTs and an interpolation layer. Similarly, the loss for the U-Net and Riesz network in the spectral flow task is taken in direct space.

In the case of the ASM, the task can be seen as a pixel-wise binary classification, so we use the Binary Cross Entropy (BCE) Loss. The network outputs a map giving the “probability” of each site to belong to the avalanche about to unfold. Here, BCE reads:

$$ \mathcal{L}(\psi, \hat{\psi}) = -\frac{1}{NL^2} \sum_{\mathbf{x},s} \left[ w_{\mathbf{x}} \psi_{\mathbf{x}}^{(s)} \log \sigma(\hat{\psi}_{\mathbf{x}}^{(s)}) + (1 - \psi_{\mathbf{x}}^{(s)}) \log (1 - \sigma(\hat{\psi}_{\mathbf{x}}^{(s)})) \right], \quad (46) $$

where we average over all pixels $\mathbf{x}$ in the lattice and over samples $s$, $\sigma$ is the sigmoid function (applied to force output into the range [0, 1]), and $w_{\mathbf{x}}$ is a reweighting factor that boosts the term corresponding to positive pixels. Indeed, the ASM dataset is highly imbalanced: since most avalanches are small, the binary target maps contain mostly zeros, so the reweighting factor is needed to prevent the network from learning the trivial solution of

predicting zeros everywhere. In practice, we use a constant factor $w_{\mathbf{x}} = w$ given by the ratio between the total number of positive and negative pixels in the target fields $\psi_{\mathbf{x}}$ from the training dataset

$$ w = \frac{\sum_{s, \mathbf{x}} (1 - \psi_{\mathbf{x}}^{(s)})}{\sum_{s, \mathbf{x}} \psi_{\mathbf{x}}^{(s)}}. \quad (47) $$

On our finite dataset, it evaluates to $w \simeq 12.25$.

# IV. RESULTS

## A. FGF phase mixing task

We present the results obtained with the U-Net, the Riesz network, and the Fourier embedding network on the FGF task. The dataset consists of 50k samples with lattice size $L = 128$, with the target being the time evolution of the field $\varphi(\mathbf{k}, t = 1)$, for a phase function (31) with parameter $\nu = 50$.

We use a 4-level U-Net with a total of 485,673 parameters. For the Riesz network, we use 4 layers, with a total of 30,057 parameters. For the Fourier embedding network, we use an embedding of $P = 100$ features, and a 2-layer MLP with both hidden layers of size 64, giving a total of 10,689 parameters.

We evaluate the performance of each model by computing the average mean squared error as a function of $|\mathbf{k}|$, normalized by the power spectrum

$$ E(|\mathbf{k}|) = \left( \frac{1}{N} \sum_s |\psi_{\mathbf{k}}^s - \hat{\psi}_{\mathbf{k}}^s|^2 \right) |\mathbf{k}|^\beta, \quad (48) $$

as shown in Fig. 5. Both U-Net and Riesz network only learn the low frequencies, which illustrates the problem of spectral bias. In contrast, using an equivalent (respectively the same) loss function, the Fourier embedding network manages to appropriately deal with spectral bias.

For the Fourier embedding network, which is the only network that manages to perform well in the standard task, we additionally perform two extrapolation experiments, by hiding a range of frequencies in the training set. In the first one, we hide (set to 0) the lowest part of the spectrum $0 < |\mathbf{k}| < 0.1$, which corresponds to extrapolation to large scales. In the second one, we hide the range $0.1 < |\mathbf{k}| < 0.2$, corresponding to interpolation to a set of middle scales. The results are shown in Fig. 6.

We observe here that the network learns for the three cases all the frequencies present in the training set. However, in both extrapolation experiments the test error is high for the hidden range of frequencies, which indicates that the model is struggling to extrapolate across scales. It still performs significantly better than a random predictor, while both the U-Net and the Riesz Network are unable (by design for the U-Net) to extrapolate.

For these three experiments, we can compare the learned function $g_\theta(\mathbf{k})$, as given by eq. 29 to the ground
<figure><img src="image_6.png" /><figcaption>FIG. 5. Normalized Mean Squared Error $E(|\mathbf{k}|)$ (equation (48)) evaluated on the test set for the Riesz network, the U-net, and the Fourier embedding network. The last one perfectly recovers the signal. The dashed gray line corresponds to the random predictor.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIG. 6. Normalized Mean Squared Error $E(|\mathbf{k}|)$ (equation (48)) evaluated on the test set, for the Fourier embedding network in the extrapolation setup. We compare the baseline training where all modes are present in the training set (full spectrum) to two extrapolation cases where ranges of frequencies are hidden. Hidden frequency ranges are highlighted with shaded areas. The dashed gray line corresponds to the random predictor.</figcaption></figure>

truth $g(\mathbf{k}) = e^{-i\nu|\mathbf{k}|^2}$. In Fig. 7 we see that the output of the model perfectly matches the target in the baseline experiment without extrapolation, but deviates from it in the hidden range for the extrapolation experiments.

In the FGF task, the energy distribution over scales is broad (see Eq. 4), which makes the learning difficult, as evidenced by the failure of U-Net and Riesz networks to succeed. Designing an architecture that explicitly operates in a Fourier representation, we induce the proper

<figure><img src="image_8.png" /><figcaption>FIG. 7. Function $g_{\theta}(|\mathbf{k}|)$ (equation 29) learned by the Fourier embedding network compared to the ground truth $g(\mathbf{k}) = e^{-i\nu|\mathbf{k}|^2}$. In the first panel (standard experiment), the true and predicted values overlap so perfectly that they cannot be distinguished. In the other panels (extrapolation experiments), the mismatch are contained within the range of frequencies that was hidden during training.</figcaption></figure>

architectural bias, which proves to be key to correct fitting.

In summary, we have introduced a testbed for extrapolation to unseen scales together with a baseline model, the Fourier Embedding Network (FEN), which shows the difficulty of this out-of-domain generalization task. We believe that the key to solving this issue is to encode the proper scaling property inherent to the task into the architecture.
## B. FGF spectral flow regression task

For this task we use a dataset of 250 samples with lattice size $L = 64$. To perform scale transforms and rotations, data are prepared using the analytic continuation in $k$ space provided by the Dirichlet kernel (equation A1 in Appendix A). The rescaling is done with a factor of $s = 0.9$ and we experiment on two cases: a dataset with no rotation, and one with a rotation of angle $\alpha = \frac{\pi}{4}$. We test the U-Net and the Riesz network, as well as the FM network designed for this task. For the U-Net and Riesz network we use the same architecture as in the phase mixing task, while for the FM network we use $M = 256$ Fourier-Mellin modes in both the radial and angular components to build the latent representation of the data, which gives $M^2 = 256^2$ weights to be learned.

In Fig. 8 we compare the three models on the baseline task without extrapolation, where all frequencies are present in the training data. As before, we measure the relative error in Fourier space, as a function of $|k|$. In the case of pure scaling (no rotation), the U-net and the Riesz network manage to learn on the low end of the spectrum, but perform significantly worse than the FM network, which despite also showing some spectral bias, is able to learn the intermediate frequencies up to $|k| = 0.45$, which is the scale at which noise is inserted in the data in the spectral flow process. In the case of scaling combined with rotation, we observe that both the U-net and Riesz network fail completely, while the FM network maintains its performance. This is likely due to the process being less local in Fourier space.

Additionally, in Fig. 9 we plot the weights learned by the FM network. Note that the process corresponds to a translation in logpolar coordinates, the rotation being a translation in the angular coordinate, and the rescaling a translation in the radial coordinate. Since the Fourier-Mellin transform is a Fourier transform applied in logpolar coordinates, this means that the transformation can be seen as a phase shift in Fourier-Mellin space. We observe that indeed, the weights learned match the theoretical phase shift that defines the process, at least around the center, where the FM modes are stronger. The difference between the two is due to residual spectral bias in Fourier Mellin caused by finite size effect when doing the cartesian to logpolar interpolation. The fact is that the loss of the learned model is lower than the loss obtained with the theoretical phase shift. The theoretical phase shift is actually not optimal because the scale transform and rotation done on the input field to generate the output data is done as already said by interpolation on the square lattice, not on the log-polar lattice. As a result the phase-shift in FM corresponding to these transformation is only approximate and the loss that we obtain is actually the lowest one that can be achieved with a purely local model in FM.

We consider two extrapolation settings: extrapolation to high frequencies i.e. super-resolution, where we hide the high end of the spectrum in the training data; and

<figure><img src="image_9.png" /><figcaption>FIG. 8. Normalized Mean Squared Error $E(|k|)$ (equation (48)) evaluated on the test set, for the Riesz network, U-net and FM network on the spectral flow FGF task, with a rescaling factor of 0.9. The first panel corresponds to scaling without rotation, and the second one to scaling combined with a rotation of $\pi/4$. The dashed gray line corresponds to the random predictor. The white line at $E(|k|) = 1 = 10^0$ corresponds to a constant predictor $\hat{\psi} = 0$: error above this line can be considered a failure.</figcaption></figure>

extrapolation to low frequencies, where we hide the low end. For the first case, frequencies above $|k| = 0.45$ correspond to the noise inserted of the rescaling, so the cut should be made below that range. In figure Fig. 10 we compare the baseline experiment and the two extrapolation cases. The FM Network is able to extrapolate both to high and low frequencies: although the error is slightly higher than the baseline in the hidden range, it is still significantly below the order of constant or random predictors. The model is especially successful at extrapolating to low frequencies, where even though the hidden range is larger, the curve is very close to the baseline.

We compare the error curves as a function of the width of the hidden range, for both extrapolation settings, in Fig. 11. Again, the performance is remarkably good, especially for extrapolation to low frequencies, even when a large portion of the spectrum is hidden during training. We remark that the number of points in the input corresponding to a given value of $|k|$ increases with $|k|$, with only a few input points accounting for the lower frequencies. This explains why we are able to hide a larger
<figure><img src="image_10.png" /><figcaption>FIG. 9. Learned weights of the FM network on an FGF transformed dataset with scaling factor $s = 0.9$ and rotation angle $\alpha = \pi/4$ (left), compared to the theoretical phase shift in FM space corresponding to this transformation (right). The training data consist of samples of size $L^2 = 64 \times 64$, and we use $M^2 = 256 \times 256$ Fourier-Mellin features.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIG. 10. Normalized Mean Squared Error $E(|\mathbf{k}|)$ (equation (48)) evaluated on the test set, for the FM network on the spectral flow FGF task, with a rescaling factor of 0.9, and a rotation angle of $\pi/4$. We compare the baseline training where all modes are present in the training data with two extrapolation settings, where we hide the high frequencies (from $k_{occ} = 0.325$) and the low frequencies (up to $k_{occ} = 0.3$) respectively with corresponding extrapolation factors up to $f_{extrapol} = 17.3$ and $f_{extrapol} = 1.29$. The dashed gray line corresponds to the random predictor. The white line at $E(|\mathbf{k}|) = 1 = 10^0$ corresponds to a constant predictor $\hat{\psi} = 0$.</figcaption></figure>

range of low frequencies.

We have shown that the FM model manages to alleviate the problem of spectral bias that affects the U-net and Riesz network at the spectral flow task, dealing well with the non-locality induced by the rotation, which causes both the other models to fail completely. Additionally and most importantly, this model performs remarkably well at extrapolation when we hide a range of frequencies from the training data, especially for large scale extrapolation with extrapolation factors of the order of $f_{extrapol} > 10$ are easily obtained, while super-resolution also work but less neatly with extrapolation factors of the order of $f_{extrapol} < 10$. This is expected

<figure><img src="image_12.png" /><figcaption>FIG. 11. Normalized Mean Squared Error $E(|\mathbf{k}|)$ (equation (48)) evaluated on the test set, for the FM network on the spectral flow FGF task, with a rescaling factor of $s = 0.9$, and a rotation angle of $\pi/4$. We compare the error curves for the two extrapolation settings (extrapolation to high (top) and low (bottom) frequencies respectively) as we increase the hidden frequency range. The range of extrapolation are therefore between $f_{extrapol} \in [11.5, 23]$ and $f_{extrapol} \in [1.29, 2.35]$ respectively for large scale extrapolation and the super-resolution settings.</figcaption></figure>

as the model is designed to learn the FM phases which correspond to fully delocalized features in $\mathbf{k}$ space, hence appropriate for extrapolation. The residual spectral bias in FM make it however more difficult for super-resolution than for large scale extrapolation, also simply because the density of modes hence of information is much larger at small than at large scale.

### C. Abelian Sandpile Model

For the ASM we consider a dataset consisting of $10^6$ samples with lattice size 32. The data is generated by running the dynamics and registering one event every hundred avalanches to provide some decorrelation. Note that the number of training samples is much smaller than the size of the configuration space of the ASM, which scales exponentially with the lattice size [27]. We present the results obtained with the U-Net, the Riesz network, and the Wavelet GNN.

In Fig. 12 we show sample outputs obtained from the
<figure><img src="image_13.png" /><figcaption>FIG. 12. Sample outputs of the Riesz network for the ASM task, on lattice size 32 × 32. The left panel corresponds to the target outputs and the right panel to the predicted ones.</figcaption></figure>

Riesz network compared to their respective ground truth. The output map is interpreted as representing the probability for each site to topple during the avalanche. As expected, the model is able to capture the starting point of the avalanche as well as some of its shape. There is a visible bias towards larger avalanches, which is due to the boosting of the positive term in the loss function. To obtain a prediction, we determine an optimal threshold by means of the receiver operating characteristic (ROC) curve, which is the curve obtained by varying the threshold of the binary classifier and plotting the true positive rate (TPR) against the false positive rate (FPR). The area under the curve (AUC) gives a (threshold-independent) metric of the quality of the classifier, with a perfect classifier having an AUC of 1, and a random one, 0.5. In Fig. 13 we show the ROC curves obtained for each of the three models. The Riesz network gives the best performance out of the three models, but only by a small margin compared to the U-Net. The Wavelet GNN has the lowest AUC of the three, but the performance is still comparable to that of the other two.

The optimal threshold is obtained by maximizing the balanced accuracy, defined in terms of the TPR and the true negative rate (TNR) as $A = \frac{1}{2}(TPR + TNR)$. The avalanche prediction obtained by applying this threshold is marked as a contour on the samples in Fig. 12. We finally compare the size distribution of the predicted avalanches to the empirical distribution obtained from the training samples. In Fig. 14 we show the size distribution obtained from each of the three models by applying the optimal threshold. Note that the domain of sizes above $S \gtrsim L^2/2$ corresponds to events occuring with a probability below $P(S) \lesssim 10^{-3}$, making them so

<figure><img src="image_14.png" /><figcaption>FIG. 13. ROC curve comparison for the U-Net, Riesz network and Wavelet GNN.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>FIG. 14. Size distribution of the predictions from the U-Net, Riesz network and Wavelet GNN, compared to the empirical distribution in the training data.</figcaption></figure>

rare that inference in this regime effectively constitutes extrapolation, given the scarcity of such samples in the training data.

We can observe that the predicted distributions have in all cases a bias towards larger avalanches, confirming the qualitative observation of the output samples. In the case of the U-Net and the Riesz network, the distribution is accurate for the small events.

To conclude with this experiment, we first observe that the Riesz network—despite its significantly smaller number of parameters—slightly outperforms the U-Net, which is heavily parameterized. This suggests that its built-in scale invariance helps for this task. On the other hand, the GNN, which also incorporates built-in scale invariance, albeit in a discrete form, exhibits lower performance. One possible explanation is that the GNN starts
from a wavelet representation, which may complicate the construction of effective filters for avalanche prediction. Indeed, other experiments with U-Nets (not shown here) using wavelet representations also yielded poor performance. In addition, edge effects in scale were likely not properly dealt with in our GNN architecture. This could deserve further development.

# V. THEORETICAL ANALYSIS

In this section we want first to characterize linear models which are scale equivariant and then study asymptotically the learning dynamics of such models in a relevant scaling limit.

## A. Scale Equivariant Linear models

We denote by FGF($\beta$) the fractional field distribution corresponding to exponent $\beta$. Let us first state the property that should obey an equivariant functional map w.r.t. scale symmetry and rotation, leaving aside translation for the moment. We assume this functional $F$ to take a 2-d field $\varphi \sim \text{FGF}(\beta)$ as input and to output another 2-field $\psi \sim \text{FGF}(\gamma)$ with exponents $\beta$ and $\gamma$, possibly distinct. Using the Fourier representation this corresponds to have (from now on we drop the tilda in $\tilde{\varphi}$ for denoting Fourier components of the field)

$$ \psi[\varphi](\mathbf{k}) = \int d^2\mathbf{k}' F(\mathbf{k}, \mathbf{k}')\varphi(\mathbf{k}') \quad (49) $$

with $\mathbf{k} = (k_x, k_y)$ a 2-d vector which we also represent in polar coordinates as $(k = |\mathbf{k}|, \alpha)$. Assuming rotation symmetry we have

$$ F(\mathbf{k}, \mathbf{k}') = F(k, k', \alpha - \alpha') $$

The equivariant property corresponds to having both $\varphi$ and $\psi$ satisfying the scale symmetry (8) but with possibly distinct exponents. In the Fourier representation this reads

$$ \varphi_s(\mathbf{k}) = \frac{1}{s^{1+\beta/2}} \varphi\left(\frac{\mathbf{k}}{s}\right) \stackrel{d}{=} \varphi(\mathbf{k}) \quad (50) $$

$$ \psi_s(\mathbf{k}) = \frac{1}{s^{1+\gamma/2}} \psi\left(\frac{\mathbf{k}}{s}\right) \stackrel{d}{=} \psi(\mathbf{k}) \quad (51) $$

where these equality hold both in distribution, meaning that the rescaled fields $\varphi_s \sim \text{FGF}(\beta)$ and $\psi_s \sim \text{FGF}(\gamma)$ follows the same distributions as the original ones as can be checked from the definition (5,6). Now imposing equivariance on $F$ means that $\varphi_s$ and $\psi_s$ should also be related by (49). This then leads to (see Appendix C 2) impose the following form for $F$:

$$ F(\mathbf{k}, \mathbf{k}') = \frac{1}{k^{\gamma/2}} \frac{f\left(\frac{\mathbf{k}'}{k}, \alpha - \alpha'\right)}{k k'} k'^{\beta/2}. \quad (52) $$

A convention is chosen here such that the kernel $f(k', \alpha - \alpha')/(kk')$ is a map from the Hilbert space $\mathcal{H}_{L_2}$ of square-integrable functions onto itself, with respect to the standard $L_2$ norm. In this form we have $F_\theta = D_F^{-\gamma/2} D_{FM} D_F^{\beta/2}$ where $D_F$ is a diagonal operator in the Fourier representation with $D_F(k) = k$, while $D_{FM}$ is strictly scale invariant and diagonal in the “Fourier-Mellin” representation. This representation uses the following orthonormal basis w.r.t. to the standard $L_2$ inner product

$$ \phi_{(\lambda, \mu)}(k, \alpha) \stackrel{\text{def}}{=} \frac{1}{2\sqrt{\pi \log(L)}} \frac{e^{i(2\pi \lambda \log(k) + \mu \alpha)}}{k} \quad (53) $$

indexed by $\lambda \in \mathbb{R}$ and $\mu \in \mathbb{Z}$ being respectively the eigenvalue of the 2-d scale and rotation operator (13,14), with $\beta = 0$ as far as the scale operator is concerned (see Appendix C 3). Here the normalization corresponds to assume an IR cutoff for $k$ equal to $1/L$. Then the relation between $f$ and $D_{FM}$ basically corresponds to a Fourier transform in scale.

To conclude this section let us remark that the form of $F$ insuring equivariance can be relaxed by allowing the functional to be equivariant up to some arbitrary phase shift i.e. gauge transformation. This does still preserve the scaling properties of the output. In this respect we extend our family of operators to

$$ F(\mathbf{k}, \mathbf{k}') = \frac{1}{k^{\gamma/2}} \frac{\left| f\left(\frac{\mathbf{k}'}{k}, \alpha - \alpha'\right) \right|}{k k'} k'^{\beta/2} e^{i(\omega_{\text{out}}(\mathbf{k}) - \omega_{\text{in}}(\mathbf{k}'))} $$

with $\omega_{\text{in,out}}$ some real functions. In the case where $F$ is a one-to-one map, corresponding for instance the end points of trajectories which do not cross each others like for an inviscid flow, then the phases can be transported along trajectories and strict equivariance is satisfied with respect to the whole group of transformations, including phase shifts. The phase mixing example considered in Section (II A) actually corresponds to a special case of such equivariant dynamical process.

## B. Projection on the feature space

Since we are doing a linear regression, the learning task amounts to finding a good projection of the underlying data-generating functional map onto an operator space defined by our set of pre-defined features, in a similar way as previous works [43]. Here we want to specify the nature of this projection. Let us consider an arbitrary functional $F$ corresponding to equation (49) and a set of $M$ features $\{f_p, p = 1, \dots, M\}$. We look for a decomposition of this kernel as

$$ F(\mathbf{k}, \mathbf{k}') = F^{\|}(\mathbf{k}, \mathbf{k}') + F^{\perp}(\mathbf{k}, \mathbf{k}') $$

such that $F^{\|}$ corresponds to the part that can be recovered with the features at hand, while $F^{\perp}$ is the orthogonal part in a sense to be specified. The set of features
defines a space $\mathcal{E}_f$ of operators onto which we expect a decomposition of the form

$$F^{\|}(\mathbf{k}, \mathbf{k}') = \sum_{p=1}^{M} w_p f_p(\mathbf{k}, \mathbf{k}'), \quad (54)$$

with coefficients $w_p \in \mathbb{R}$. These are obtained by orthonormal projection onto $\mathcal{E}_f$, w.r.t. an inner product defined on a larger operator space containing both $F$ and $\mathcal{E}_f$. For this it is required that two operators $F_1$ and $F_2$ are orthogonal whenever the corresponding $\psi_1[\varphi]$ and $\psi_2[\varphi]$ obtained from (49) are orthogonal in average, i.e. that:

$$\sum_{\mathbf{k} \in \Omega} |\mathbf{k}|^{\gamma} \mathbb{E}_{\varphi \sim FGF(\beta)} [\psi_1[\varphi](\mathbf{k}) \psi_2^{*}[\varphi](\mathbf{k})] = 0$$

where $*$ denotes complex conjugation. Here we implicitly use the inner product associated to the energy (6) of the $FGF(\gamma)$ distribution. This then corresponds to have

$$\sum_{\mathbf{k}, \mathbf{k}' \in \Omega} \frac{k^{\gamma}}{k'^{\beta}} F_1(\mathbf{k}, \mathbf{k}') F_2^*(\mathbf{k}, \mathbf{k}') = 0.$$

Consequently we define an inner product between operators adapted to input fields $\varphi \sim FGF(\beta)$ and output fields $\psi \sim FGF(\gamma)$ as

$$\langle F_1, F_2 \rangle_{\beta, \gamma} \stackrel{\text{def}}{=} \frac{1}{L^2} \sum_{\mathbf{k}, \mathbf{k}' \in \Omega} \frac{k^{\gamma}}{k'^{\beta}} F_1(\mathbf{k}, \mathbf{k}') F_2^*(\mathbf{k}, \mathbf{k}'), \quad (55)$$

where the $L^2$ normalization stands here for later convenience. Equipped with this definition, the weights of the projection (54) are given by

$$w_p = \sum_{q=1}^{M} [C^{-1}]_{pq} \langle F, f_q \rangle_{\beta, \gamma}, \quad p = 1, \dots, M. \quad (56)$$

with $C$ the covariance of the features, those coefficients are given by

$$C_{pq} = \langle f_p, f_q \rangle_{\beta, \gamma}, \quad p, q = 1, \dots, M. \quad (57)$$

This corresponds to perform the regression with an infinite amount of data. In practice of course we have a finite number of data points and the regression is the result of a learning process which is described in the next section.

### C. Linearized learning dynamics

Let us now state the learning dynamics in full generality. From now on we take into account the fact that we are on a finite lattice $\Omega$, so that all integrals are replaced by summation of $\mathbf{k}$ over $\Omega$. The data consists in a set of $N$ pairs of fields $\{ (\varphi^{(s)}, \psi^{(s)}), s = 1, \dots, N \}$ related by

$$\psi(\mathbf{k}) = \sum_{\mathbf{k}' \in \Omega} K(\mathbf{k}, \mathbf{k}') \varphi(\mathbf{k}') + \sigma(\mathbf{k}) \epsilon(\mathbf{k}) \quad (58)$$

with $\varphi \sim FGF(\beta)$ and $\psi \sim FGF(\gamma)$ the joint distribution being denoted $\sim FGF(\beta, \gamma)$, $\Omega \subset \mathbb{Z}^2$, with $|\Omega| = L^2$ and where $K$ a kernel of the form (52) and $\epsilon(\mathbf{k}) = \mathcal{N}(0, 1)$ is iid, with both $K$ and $\sigma(\mathbf{k})$ satisfying the necessary constraints to ensure that $\psi \sim FGF(\gamma)$. In particular we have $\sigma(\mathbf{k}) \propto 1/k^{\gamma/2}$. The model $\psi_{\theta}[\varphi]$ to be learned is defined by

$$\psi_{\theta}[\varphi](\mathbf{k}) = \sum_{\mathbf{k}' \in \Omega} F_{\theta}(\mathbf{k}, \mathbf{k}') \varphi(\mathbf{k}') \quad (59)$$

$$\text{with } F_{\theta}(\mathbf{k}, \mathbf{k}') = \sum_{p}^{M} f_p(\mathbf{k}, \mathbf{k}') \theta_p \quad (60)$$

where $\theta$ is the parameter vector, assumed to be of dimension $M$. In full generality, $F_{\theta}(\mathbf{k}, \mathbf{k}')$ could in principle be a more complex architecture, with more parameters. Here we are interested by the dynamics of the linear regime, where $f_p(\mathbf{k}, \mathbf{k}')$'s are independent from $\theta$. We minimize the $L_2$ loss w.r.t. $\theta$:

$$\mathcal{L}(\theta) = \frac{1}{2} \mathbb{E}_{(\varphi, \psi) \sim \widehat{FGF}(\beta, \gamma)} \left[ \sum_{\mathbf{k} \in \Omega} |\psi_{\theta}[\varphi](\mathbf{k}) - \psi(\mathbf{k})|^2 \right], \quad (61)$$

where $\widehat{FGF}(\beta, \gamma)$ represents the empirical joint distribution of $\varphi$ and $\psi$ related by (58).

Let us turn now to learning dynamics. For reasons which will become more apparent when considering the scaling limits, we introduce a learning rate scaling like $\eta/L^2$ with $\eta = \mathcal{O}(1)$ non-negative. The gradient dynamics then reads

$$\dot{\theta} = -\frac{\eta}{L^2} \frac{\partial \mathcal{L}(\theta)}{\partial \theta} = -\frac{\eta}{L^2} \mathbb{E}_{(\varphi, \psi) \sim \widehat{FGF}(\beta, \gamma)} \left[ \sum_{(\mathbf{k}, \mathbf{k}') \in \Omega^2} \mathbf{f}(\mathbf{k}, \mathbf{k}') \varphi(\mathbf{k}') (\psi_{\theta}[\varphi](\mathbf{k}) - \psi(\mathbf{k})) \right],$$

with $\mathbf{f}(\mathbf{k}, \mathbf{k}')$ the vector of features $f_p(\mathbf{k}, \mathbf{k}')$.

The linear dynamics is of the form

$$\eta^{-1} \dot{\theta} = -\hat{C}\theta + \hat{Z}$$

with $\hat{C}$ an $M \times M$ empirical matrix and $\hat{Z}$ an $M$-dimensional empirical vector given by

$$\hat{C}_{pq} = \frac{1}{L^2} \sum_{\mathbf{k}, \mathbf{k}', \mathbf{k}'' \in \Omega^3} f_p(\mathbf{k}, \mathbf{k}') f_q(\mathbf{k}, \mathbf{k}'') \times \mathbb{E}_{\varphi \sim FGF(\beta)} [\varphi(\mathbf{k}') \varphi(\mathbf{k}'')] \quad (62)$$

$$\hat{Z}_p = \frac{1}{L^2} \sum_{\mathbf{k}, \mathbf{k}' \in \Omega^2} f_p(\mathbf{k}, \mathbf{k}') \mathbb{E}_{\varphi, \psi \sim \widehat{FGF}(\beta, \gamma)} [\psi(\mathbf{k}) \varphi(\mathbf{k}')] \quad (63)$$

leading to

$$\theta(t) = \hat{G}(t) \hat{Z},$$
with

$$ \hat{G}(t) \stackrel{\text{def}}{=} \hat{C}^{-1} (1 - e^{-\eta \hat{C} t}) \qquad (64) $$

the propagator in time of the learning process. The difference with the usual linear regression setting, discussed in [44] is that we have now $|\Omega| = L^2$ regressions in parallel indexed by $\mathbf{k}$ instead of a single scalar one, and the feature matrix is averaged over all these distinct regressions by summing over $\mathbf{k}$ in (62,63). Also as discussed in [44], we see in (64) that early stopping acts as a regularization, since at time $t$, only modes with eigenvalues larger than $1/(\eta t)$ are learned. In our case this can lead to a strong spectral bias issue, when the modes of $\hat{C}$ decay as a power law in $k$. Just notice that the population matrix $C$ of $\hat{C}$ corresponding to taking the limit $N \to \infty$ at fixed $M$ is actually given by (57). This indicates that in order to efficiently address the spectral bias issue, it may be sufficient to use a preconditioner that would orthonormalize the features w.r.t. the inner product introduced above in (55), leading to $C = \mathbb{I}$. A second issue is the extrapolation one. As previously stated we have $L^2$ regressions problems indexed by $\mathbf{k}$ and for this we need to share information between those regression problems. From the previous discussion on equivariance, the obvious way to proceed is to look for features in the form of smooth function of $\mathbf{k}$ that diagonalize the kernel $K$, i.e. functions having the functional form (52). Hence our use of the Fourier-Mellin basis in the building of our FM Network to address the spectral flow task.

### D. Asymptotic limits of learning trajectories

We are now interested to obtain learning trajectories in some meaningful asymptotic limits. From the previous discussion, it appears that the meaningful scaling parameter is given by

$$ \rho = \frac{N N_s}{P} \qquad (65) $$

where $N$ is the number of independent samples, $N_s$ the number of independent scalar observable per sample and $P$ the number of parameters of the model. In the context of FGF regression, the Fourier components are independent so in principle we should have $N_s = L^2$. Actually we only have $N_s = s^2 L^2$ for the spectral flow experiment with $s < 1$ since only the output modes $\psi(\mathbf{k})$ of the output with $\max(|k_x|, |k_y|) < s/2$ contain information. In the context of extrapolation scenario considered for the FGF experiments, some of the components are hidden during the train and the effective number of observations is further reduced as will be detailed in Section V E. In absence of inductive bias for the inference model, the feature space is supposed to span an operator space of dimension $L^4$, which means that $P = O(L^4)$ in this case. If instead some prior knowledge is given, like scale and rotation (equ-) invariance, the operator could

be considered in the diagonal form of some specific representation, like Fourier or Fourier-Mellin in our study cases, and then the number of parameters is reduced to $P = O(L^2)$. This allows us to consider a meaningful asymptotic regime where $N = O(1)$, enabling the analysis of global quantities such as training and test errors. In such case we recover a dynamical behavior corresponding to the population limit when $N$ becomes large. When considering instead local quantities like mode dependent train and test error, these should become deterministic when $N, P \to \infty$ in the proportional scaling regime $N/P$ fixed. This will become more obvious in the following.

Let us look into this in more details. Call $\psi_{\theta_t}[\varphi]$ the model learned at learning time $t$. First, the evolution of the training and test errors as functions of time can be given in a general form in terms of the time propagator (64). This leads to equations (D1,D2) given in Appendix D 2. when considering the mean squared error (MSE) as the default metric, consistent with the loss function (61). These equations actually simplify significantly when considering instead the metric based on the FGF energy (6), leading to the following definition for the train and test errors as a function of $t$

$$ E_{\text{train}}(t) = \mathbb{E}_{(\varphi, \psi) \sim \widehat{\text{FGF}}(\beta, \gamma)} \left[ \sum_{\mathbf{k} \in \Omega} |\mathbf{k}|^{\gamma} |\psi_{\theta_t}[\varphi](\mathbf{k}) - \psi(\mathbf{k})|^2 \right], $$

$$ E_{\text{test}}(t) = \mathbb{E}_{(\varphi, \psi) \sim \text{FGF}(\beta, \gamma)} \left[ \sum_{\mathbf{k} \in \Omega} |\mathbf{k}|^{\gamma} |\psi_{\theta_t}[\varphi](\mathbf{k}) - \psi(\mathbf{k})|^2 \right]. $$

while modifying the loss (61) accordingly. The corresponding gradient descent has also the advantage to compensate for the spectral bias which originates from the power-law behavior of the spectrum and to simplify the learning equations as explained in Appendix D 2, so we will concentrate only on this formulation which corresponds to equations (D3,D4). While these equations can be dealt with in principle in full generality, once we are given the spectrum of the population matrix, for sake of simplicity (in particular to deal with the first term of (D4)) we focus on the case where $C = \mathbb{I}$. In that case, using random matrix theory we eventually arrive at closed form expressions for the train and test errors as function of time (see Appendix D). The way to proceed is quite standard (see e.g. [43, 45] and refs. herein respectively in the context of scalar regression and energy based models). For this we need first to introduce the mean spectral density $\bar{\nu}$ of $\hat{C}$, which for $C = \mathbb{I}$ is the ordinary Marchenko-Pastur distribution:

$$ \bar{\nu}(x) = (1 - \rho) \mathbb{I}_{\{\rho < 1\}} \delta(x) + \frac{\rho}{2\pi x} \sqrt{(x - x_{-})(x_{+} - x)}, $$

with

$$ x_{\pm} = \left(1 \pm \frac{1}{\sqrt{\rho}}\right)^2. $$

Then we need to introduce two quantities depending on
the decomposition of

$$K = K^{\parallel} + K^{\perp}$$

involved in (59), corresponding respectively to the component $K^{\parallel}$ on the feature space and its orthogonal space. Using the inner product (55) defined above for operators we call

$$\kappa \stackrel{\text{def}}{=} \langle K^{\parallel}, K^{\perp} \rangle_{\beta, \gamma}$$

and $\sigma_{\text{eff}}^2$ the effective variance of the noise which emerges from the fact that the regression is misspecified. By definition $K^{\perp}$ cannot be recovered from the features and adds an additional contribution to the original noise entering in the definition of the data model (58)

$$\sigma_{\text{eff}}^2 \stackrel{\text{def}}{=} \sum_{k \in \Omega} k^{\gamma} \sigma^2(k) + \langle K^{\perp}, K^{\perp} \rangle_{\beta, \gamma}. \quad (66)$$

This noise variance can actually be decomposed component-wise onto the modes as

$$\sigma_{\text{eff}}^2(k) = k^{\gamma} \sigma^2(k) + \frac{1}{L^2} \sum_{\mathbf{k}' \in \Omega} \frac{k^{\gamma}}{k'^{\gamma}} K^{\perp}(\mathbf{k}, \mathbf{k}')^2.$$

Assuming without loss of generality $\langle K, K \rangle_{\beta, \gamma} = 1$, we have $\langle K^{\perp}, K^{\perp} \rangle_{\beta, \gamma} = 1 - \kappa$.

These two quantities fully determine the learning dynamics at large scale for the train, and test errors. Thanks to the preconditioning corresponding to the input and output field scaling behavior, the $k$ dependency in equations (D1,D2) disappears and we are eventually left with the following simplified expression for the train and test error in the asymptotic limit (see Appendix D):

$$E_{\text{train}}(t) = \kappa \int_{0}^{\infty} \bar{\nu}(dy)y[1 - yj(y,t)]^2 + \sigma_{\text{eff}}^2 \left[1 - \frac{1}{\rho} + \frac{1}{\rho} \int \bar{\nu}(dy)[1 - yj(y,t)]^2\right] \quad (67)$$

$$E_{\text{test}}(t) = \kappa \int_{0}^{\infty} \bar{\nu}(dy)[1 - yj(y,t)]^2 + \sigma_{\text{eff}}^2 \left[1 + \frac{1}{\rho} \int dy \bar{\nu}(dy)yj(y,t)^2\right] \quad (68)$$

where

$$j(x, t) = \frac{1 - e^{-xt}}{x}$$

with learning rate $\eta = 1$ by convention, is the function controlling the dynamics of $\theta$ in equation (64). As shown in the Appendix D 1, this simplified case actually corresponds to make use of a preconditioner adapted to the spectrum of the output.

### E. The case of FGF spectral flow in the extrapolation regime

Let us consider now what happens in the extrapolation setting of the spectral flow experiment considered in Section IV B. The range of scales on the square lattice is $k_{x,y} \in [k_{\min}, k_{\max}]$ with $k_{\min} = \frac{1}{L}$ and $k_{\max} = 0.5$. When considering the rescaling, which maps the component $\mathbf{k}$ to $s\mathbf{k}$ with $s < 1$, all the components $\psi(\mathbf{k})$ of the output, situated in the domain $\max(|k_x|, |k_y|) > sk_{\max}$ are randomly generated, which corresponds to inserting the amount of noise (66)

$$\sigma_{\text{eff}}^2 = 1 - s^2 \quad (s < 1) \quad (69)$$

into the output field, when this one is assumed to be normalized to one. In addition, in the extrapolation regime, let us call $k_{\text{occ}}$ the occultation scale, which for the super-resolution experiment corresponds to hide the components $\mathbf{k}$ in the domain $\max(|k_x|, |k_y|) > k_{\text{occ}}$ in the fields, while the large scale extrapolation experiments corresponds to hide components in the domain $\max(|k_x|, |k_y|) < k_{\text{occ}}$. The set of meaningful observation (in the training set) of the output field in the super-resolution case then corresponds to $\max(|k_x|, |k_y|) \in [k_{\text{occ}}, sk_{\max}]$, others components either corresponding to noise or being hidden. Instead in the large scale extrapolation the meaningful observations correspond to $\max(|k_x|, |k_y|) > sk_{\text{occ}}$ which justifies the definition (15) of the extrapolation factor given in Section II A. To see the effect of extrapolation on the performance of the prediction we introduce also the ratio

$$r \stackrel{\text{def}}{=} \frac{k_{\text{occ}}}{k_{\max}}$$

which is more convenient. As a function of $r$ and $s$ the equations (67,68) remain unchanged except for the coefficient $\sigma_{\text{eff}}(s)$ given by (69) and the samples/parameters ratio (65) of the global feature matrix

$$\rho(r, s) = \rho_0 \times \begin{cases} (s^2 - r^2) & \text{(large scale extrapolation)} \\ (sr)^2 & \text{(super-resolution)} \end{cases}$$

where $\rho_0$ is a default ratio. As we see, for the large scale extrapolation there is no information to process when $r \ge s$ with a linear decay close to this point while for super-resolution the information vanishes quadratically with $r$.

The expected behavior of the test error in absence of any spectral bias of the features, i.e. corresponding to $C = \mathbb{I}$ is displayed in Figure 15 for the large scale extrapolation setting (the super-resolution case being similar). Note that this curves are generic and basically assume that there is some information, not local in scale, which couples the various regressions of $\psi(\mathbf{k})$ to be done at different scales. Remarkably this tells us that in principle for sufficiently large input and output dimension
<figure><img src="image_16.png" /><figcaption>FIG. 15. Behavior of the test error measured when varying the level of extrapolation $r$ for $s = 0.8$ and $\kappa = 1$ in the large scale extrapolation setting. The inset shows how the optimal stopping point and the corresponding test error behave with $r$. The spike of $E_{opt}$ corresponds to the value of $r$ corresponding $\rho(r, s) = 1$ i.e. the interpolation point separating the under ($\rho > 1$) from the over-parameterized regime.</figcaption></figure>

($L \gg 1$), inference in the extreme extrapolation regime ($r$ close to $s$ in the large scale extrapolation task or close to zero in the super-resolution task) should be achievable with only a few training instances ($N$ of the order of one).

In practice, in our experiments we do not find us exactly in this idealized situation. For the phase mixing problem, the phase shift function $g_{\theta}(\mathbf{k})$ to be learned in (29) is not smooth enough to be decomposable into a few basic (Fourier modes in our case) delocalized functions of $\mathbf{k}$. Concerning the FGF spectral flow process, which in principle should be closer to the ideal case, it is still displaying spectral biases issues resulting from finite size effect on the square lattice to the log-polar lattice interpolation. Our current implementation of the Fourier-Mellin transform is computationally prohibitive, primarily due to the Dirichlet kernel and its pseudo-inverse, which represent fully connected layers. Developing a sparse implementation of these layers would enable us to scale the model sufficiently to observe the behaviors described on Figure 15.

## VI. CONCLUSION AND PERSPECTIVES

In this work we have investigated how scale invariance can be leveraged to perform extrapolation across scales on two study cases of scale invariant processes: linear FGF processes and the ASM. At odds with most of the existing literature on this subject, we consider a setting in which scale invariance is realized at a probabilistic level, meaning that we look for equivariance properties (2) of the architecture at the level of the probability distributions of the input and output fields.

The two examples we consider represent two extreme cases in terms of difficulty. On the one hand, for the FGF

we consider the two possible types of linear elementary processes leaving the FGF distribution invariant: a phase shift, and a spectral flow corresponding to rescaling and rotation of the input in momentum space. On the other hand the ASM is a non-linear map with input and output fields (height map and avalanche cluster) pertaining to very different kinds of scale-free fields. The idea of considering FGF processes is to identify basic constraints that the neural operator should obey. For each FGF process we develop new architectures – the Fourier embedding and Fourier-Mellin network respectively – which are adapted to each of these elementary transformations. This opens the possibility to design more complex architectures by combining these building blocks, in order to address more complex inference tasks. Already we see that these elementary tasks confound sophisticated models like the U-Net or the Riesz network. They suffer either from very strong spectral bias in the phase mixing task, either fail entirely for the spectral flow task, not to mention the corresponding extrapolation tasks. For the FGF, while the basic tasks look straightforward, they are not in practice and our work provides a proof of concept demonstrating that the proposed elementary building blocks are effective for extrapolation. These building blocks could be used as basic layers of more complex deep models, in combination with non-linear activations like e.g. ReLU which would preserve the scale symmetry, in some way generalizing the Riesz network. This is left for further investigations.

From the theoretical analysis of Section V, we reach two important points. First it is clearly seen that a pre-conditioner for the gradient descent is mandatory to avoid filtering the small scale information. The relevant preconditioner can be directly deduced from the knowledge of the statistical properties of the input and output fields, which can be either considered to be known in advance, as we do for the FGF, either estimated from the data. Second, imposing equivariance with respect to scale lets us reduce the number of features by a factor $L^2$, as discussed in Section V D, and ensures the possibility of extrapolation, in principle with only a few examples as stated in Section D 3.

For the ASM instead the question of extrapolation remains open. Naive scale invariance alone (either Riesz or our GNN) is clearly not the solution and we did not yet attempt to test scale equivariant models, as we neither expect these to be sufficient to work. Some additional sophistication obviously needs to be added. Indeed, we believe that the model needs to learn an equivariant latent representation analogous to what the spatial RG might deliver manually, in order to be able to decode across scales the output field in a top-down auto-regressive manner, somewhat similar in spirit to what has been done successfully in [19] in the context of generative models. In this respect the ASM is particularly complex as the coarse graining operation is highly non-trivial, as shown in [34]. By contrast, in the FGF tasks considered here, the coarse graining step is directly given by the anoma-
lous scaling of the field and does not involve any latent representation. Future investigations should concentrate on problems of intermediate difficulty between these two, like for instance percolation problems where the coarse graining needs to be learned but is less involved than for ASM. Additionally a pragmatic strategy should consider first to learn the coarse grained latent representation in a supervised manner as a starting point, by learning directly pre-defined spatial RG flows.

## ACKNOWLEDGMENTS

Authors thank Misaki Ozawa and Sergio Chibbaro for insightful discussions and acknowledge financial support by the French ANR grant Scalp (ANR-24-CE23-1320).

## REFERENCES

[1] G. Carleo, I. Cirac, K. Cranmer, L. Daudet, M. Schuld, N. Tishby, L. Vogt-Maranto, and L. Zdeborová, Machine learning and the physical sciences, Reviews of Modern Physics **91**, 045002 (2019).

[2] C. Meng, S. Griesemer, D. Cao, S. Seo, and Y. Liu, When physics meets machine learning: A survey of physics-informed machine learning, Machine Learning for Computational Science and Engineering **1**, 20 (2025).

[3] N. Rahaman, A. Baratin, D. Arpit, F. Draxler, M. Lin, F. Hamprecht, Y. Bengio, and A. Courville, On the spectral bias of neural networks, in ICML (PMLR, 2019) pp. 5301–5310.

[4] Y. LeCun, B. Boser, J. Denker, D. Henderson, R. Howard, W. Hubbard, and L. Jackel, Backpropagation applied to handwritten zip code recognition, Neural computation **1**, 541 (1989).

[5] Z. Li, F. Liu, W. Yang, S. Peng, and J. Zhou, A survey of convolutional neural networks: analysis, applications, and prospects, IEEE transactions on neural networks and learning systems **33**, 6999 (2021).

[6] F. Scarselli, M. Gori, A. Tsoi, M. Hagenbuchner, and G. Monfardini, The graph neural network model, IEEE transactions on neural networks **20**, 61 (2008).

[7] M. M. Bronstein, J. Bruna, T. Cohen, and P. Veličković, Geometric deep learning: Grids, groups, graphs, geodesics, and gauges, arXiv preprint arXiv:2104.13478 (2021).

[8] N. Thomas et al., Tensor field networks: Rotation- and translation-equivariant neural networks for 3D point clouds, arXiv.

[9] M. Geiger, T. Smidt, A. M., B. K. Miller, W. Boomsma, B. Dice, K. Lapchevskyi, M. Weiler, M. Tyszkiewicz, S. Batzner, D. Madisetti, M. Uhrin, J. Frellsen, N. Jung, S. Sanborn, M. Wen, J. Rackers, M. Rød, and M. Bailey, Euclidean neural networks: e3nn (2022).

[10] D. Worrall and M. Welling, Deep scale-spaces: Equivariance over scale, Advances in Neural Information Processing Systems **32** (2019).

[11] I. Sosnovik, M. Szmała, and A. Smeulders, Scale-equivariant steerable networks, arXiv preprint arXiv:1910.11093 (2019).

[12] M. Jaderberg, K. Simonyan, A. Zisserman, et al., Spatiotemporal transformer networks, Advances in neural information processing systems **28** (2015).

[13] A. Kanazawa, A. Sharma, and D. Jacobs, Locally scale-invariant convolutional neural networks, arXiv preprint arXiv:1412.5104 (2014).

[14] T.-Y. Lin, P. Dollár, R. Girshick, K. He, B. Hariharan, and S. Belongie, Feature pyramid networks for object detection, in *Proceedings of the IEEE conference on computer vision and pattern recognition* (2017) pp. 2117–2125.

[15] E. H. Adelson, C. H. Anderson, J. R. Bergen, P. J. Burt, and J. M. Ogden, Pyramid methods in image processing, RCA engineer **29**, 33 (1984).

[16] T. Lindeberg, Scale-covariant and scale-invariant Gaussian derivative networks, Journal of Mathematical Imaging and Vision **64**, 223 (2022).

[17] T. Barisin, K. Schladitz, and C. Redenbach, Riesz networks: Scale-invariant neural networks in a single forward pass, Journal of Mathematical Imaging and Vision **1** (2024).

[18] M. A. Rahman and R. A. Yeh, Truly scale-equivariant deep nets with Fourier layers, Advances in Neural Information Processing Systems **36** (2024).

[19] T. Marchand, M. Ozawa, G. Biroli, and S. Mallat, Multiscale data-driven energy estimation and generation, Physical Review X **13**, 041038 (2023).

[20] F. Guth, S. Coste, V. De Bortoli, and S. Mallat, Wavelet score-based generative modeling, Advances in neural information processing systems **35**, 478 (2022).

[21] A. Anandkumar, K. Azizzadenesheli, K. Bhattacharya, N. Kovachki, Z. Li, B. Liu, and A. Stuart, Neural operator: Graph kernel network for partial differential equations, in *ICLR 2020 workshop on integration of deep neural models and differential equations* (2020).

[22] Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandkumar, Fourier neural operator for parametric partial differential equations, in *International Conference on Learning Representations* (2021).

[23] L. Lu, P. Jin, G. Pang, Z. Zhang, and G. Karniadakis, Learning nonlinear operators via deeponet based on the universal approximation theorem of operators, Nature machine intelligence **3**, 218 (2021).

[24] N. Kovachki, Z. Li, B. Liu, K. Azizzadenesheli, K. Bhattacharya, A. Stuart, and A. Anandkumar, Neural operator: Learning maps between function spaces with applications to pdes, Journal of Machine Learning Research **24**, 1 (2023).

[25] A. Lodhia, S. Sheffield, X. Sun, and S. Watson, Fractional Gaussian fields: A survey, Probability Surveys **13** (2016).

[26] D. Dhar, Self-organized critical state of sandpile automaton models, Physical Review Letters **64**, 1613 (1990).

[27] D. Dhar, The abelian sandpile and related models, Physica A: Statistical Mechanics and its applications **263**, 4 (1999).

[28] F. S. Pezzicoli, G. Charpiat, and F. P. Landes, Rotation-
equivariant graph neural networks for learning glassy liquids representations, SciPost Physics 16, 136 (2024).

[29] F. Pezzicoli, *Statistical Physics-Machine Learning Interplay: from Addressing Class Imbalance with Replica Theory to Predicting Dynamical Heterogeneities with SE (3)-equivariant Graph Neural Networks*, Ph.D. thesis, Université Paris-Saclay (2024).

[30] P. Bak, C. Tang, and K. Wiesenfeld, Self-organized criticality: An explanation of the 1/f noise, Physical review letters 59, 381 (1987).

[31] D. L. Turcotte, Self-organized criticality, Reports on progress in physics 62, 1377 (1999).

[32] N. W. Watkins, G. Pruessner, S. C. Chapman, N. B. Crosby, and H. J. Jensen, 25 years of self-organized criticality: Concepts and controversies, Space Science Reviews 198, 3 (2016).

[33] K. Christensen and Z. Olami, Sandpile models with and without an underlying spatial structure, Physical Review E 48, 3361 (1993).

[34] A. Vespignani, S. Zapperi, and L. Pietronero, Renormalization approach to the self-organized critical behavior of sandpile models, Physical Review E 51, 1711 (1995).

[35] B. Gutenberg and C. F. Richter, Earthquake magnitude, intensity, energy, and acceleration: (second paper), Bulletin of the seismological society of America 46, 105 (1956).

[36] O. Ronneberger, P. Fischer, and T. Brox, U-net: Convolutional networks for biomedical image segmentation, in *International Conference on Medical image computing and computer-assisted intervention* (Springer, 2015) pp. 234–241.

[37] M. Ma, Q. Zhu, Y. Zhan, Z. Yin, H. Wang, and Y. Zheng, Robustifying Fourier features embeddings for implicit neural representations, arXiv preprint arXiv:2502.05482 (2025).

[38] S. Derrode and F. Ghorbel, Robust and efficient Fourier–Mellin transform approximations for gray-level image reconstruction and complete invariant description, Computer vision and image understanding 83, 57 (2001).

[39] S. Mallat, *A wavelet tour of signal processing* (Elsevier, 1999).

[40] I. W. Selesnick, R. G. Baraniuk, and N. C. Kingsbury, The dual-tree complex wavelet transform, IEEE signal processing magazine 22, 123 (2005).

[41] What we call diagonal is actually what comes from combining horizontal and vertical basis elements.

[42] M. Schlichtkrull, T. N. Kipf, P. Bloem, R. Van Den Berg, I. Titov, and M. Welling, Modeling relational data with graph convolutional networks, in *The semantic web: 15th international conference, ESWC 2018, Heraklion, Crete, Greece, June 3–7, 2018, proceedings* 15 (Springer, 2018) pp. 593–607.

[43] C. Furtlehner, Free dynamics of feature learning processes, J.Stat.Phys. 190, 51 (2023).

[44] M. S. Advani, A. M. Saxe, and H. Sompolinsky, High-dimensional dynamics of generalization error in neural networks, Neural Networks 132, 428 (2020).

[45] G. Catania, A. Decelle, C. Furtlehner, and B. Seoane, A theoretical framework for overfitting in energy-based modeling, in *International Conference on Machine Learning* (PMLR, 2025).

[46] If we exclude the singular scaling that reduces space to a point, homoteties centered on the origin form a group.

[47] V. A. Marchenko and L. A. Pastur, Distribution of eigen-

values for some sets of random matrices, Mathematics of the USSR-Sbornik 1, 457 (1967).

[48] O. Ledoit and S. Péché, Eigenvectors of some large sample covariance matrix ensembles, Probability Theory and Related Fields 151, 233 (2011).

[49] R. Couillet, M. Debbah, and J. W. Silverstein, A deterministic equivalent for the analysis of correlated mimo multiple access channels, IEEE Transactions on Information Theory 57, 3493 (2011).

[50] W. Hachem, P. Loubaton, and J. Najim, Deterministic equivalents for certain functionals of large random matrices, The Annals of Applied Probability 17, 875 (2007).

[51] C. Chouard, *Sample covariance random matrices arising in artificial neural networks*, Ph.D. thesis, Université Paul Sabatier-Toulouse III (2023).

[52] Z. Burda, A. Görlich, A. Jarosz, and J. Jurkiewicz, Signal and noise in correlation matrix, Physica A: Statistical Mechanics and its Applications 343, 295 (2004).

[alpha]
## Appendix A: Model details

In this appendix we give additional details on some of the various architectures used throughout this work. In particular, we include the derivation of the Dirichlet kernel method employed by the FM network to perform the Fourier-Mellin transform.

### 1. U-Net

U-nets are a type of CNN architecture created for image segmentation, which is the task of highlighting or extracting an object of interest from an image.

These networks consist of a series of “downwards” layers followed by an equal number of “upwards” layers, forming in this way a hierarchical structure in which a representation at each scale is constructed from the previous one. This multiscale structure makes the U-Net a good baseline for the datasets treated in this work. They do not employ however any sort of scale invariance, as the weights learned at each scale are independent from each other. A diagram of our instantiation of this family of architectures is shown in Fig. 16.

The downwards layers apply a convolution to the input, then they coarse-grain it via a pooling operation, reducing the spatial dimensions by half. One of the key differences between the U-net and usual CNNs is that U-net records the features learned at each downwards step, to re-use that information later, in the upscaling steps.

The upwards layers perform an upscaling of the input, where a transposed convolution is applied as shown in figure 16 in order to double the spatial size, and this is convoluted with the output of the downwards layer of the same size, to aggregate the information corresponding to that scale.

### 2. Riesz network

The Riesz network is designed to be scale invariant due to its use of scale invariant filters. In Fig. 17 we show the five Riesz filters (up to order 2) that we use in practice to build the network, applied to a sample of the FGF.

The shortcoming of this network is the limited expressivity, since the filters are fixed, and the network learns only their linear combinations. Two simple ways of increasing the expressivity are to add more layers and more channels at each layer, or to include higher order filters. However, both these options are computationally costly, due to the need to employ Fourier transforms at each step in order to apply the filters.

### 3. FM network

Here we derive the expression of the Dirichlet kernel, which we use to perform the interpolation of the FGF

data from a cartesian lattice to logpolar coordinates.

Suppose we have a FGF sample on a square lattice of size L, denoted by {$\varphi(\mathbf{r}_i); \mathbf{r}_i \in \Omega_L$}. We take the convention $\mathbf{r}_i = (x_j, y_j)$,

$$x_j = -\frac{1}{2} - \frac{1}{2L} + \frac{j}{L}; \quad y_j = -\frac{1}{2} - \frac{1}{2L} + \frac{j}{L}$$

The discrete Fourier transform is given by

$$\varphi(\mathbf{k}_j) = \frac{1}{L} \sum_{\mathbf{x}_i \in \Omega_L} \varphi(\mathbf{r}_i) e^{i 2\pi \mathbf{k}_j \cdot \mathbf{x}_i},$$

with $\mathbf{k}_j \in \Omega_L$ the reciprocal square lattice.

We consider the analytic extension of the discrete Fourier transform, given by

$$\tilde{\varphi}(\mathbf{k}) = \frac{1}{L} \sum_{\mathbf{r}_i \in \Omega_L} \varphi(\mathbf{r}_i) e^{i 2\pi \mathbf{k} \cdot \mathbf{r}_i}.$$

To perform the interpolation between the square and log-polar lattice, we evaluate this field at the points of the logpolar lattice of size N, which we will denote $\mathbf{k}' \in \Omega'_N$.

$$\tilde{\varphi}(\mathbf{k}') = \frac{1}{L} \sum_{\mathbf{r}_i \in \Omega_L} \varphi(\mathbf{r}_i) e^{i 2\pi \mathbf{k}' \cdot \mathbf{r}_i}.$$

We then write $\varphi(\mathbf{r}_i)$ in terms of its Fourier coefficients on the square lattice

$$\tilde{\varphi}(\mathbf{k}') = \frac{1}{L^2} \sum_{\mathbf{r}_i \in \Omega_L} \sum_{\mathbf{k} \in \Omega_L} \tilde{\varphi}(\mathbf{k}) e^{i 2\pi (\mathbf{k}' - \mathbf{k}) \cdot \mathbf{r}_i}$$

and define the Dirichlet kernel $D$ as

$$\tilde{\varphi}(\mathbf{k}') = \sum_{\mathbf{k} \in \Omega_L} D(\mathbf{k}' - \mathbf{k}) \tilde{\varphi}(\mathbf{k}), \quad (A1)$$

so we have

$$D(\mathbf{k}' - \mathbf{k}) = \frac{1}{L^2} \sum_{\mathbf{r}_i \in \Omega_L} e^{i 2\pi (\mathbf{k}' - \mathbf{k}) \cdot \mathbf{r}_i}$$

We can compute the geometric sum, using the convention for the values of $\mathbf{r}_i$

$$D(\mathbf{q}) = D_{1d}(q_x) D_{1d}(q_y),$$

with the 1d Dirichlet kernel $D_{1d}$ defined as

$$D_{1d}(q) = \begin{cases} \frac{1}{L} \frac{\sin(\pi L q)}{\sin(\pi q)} & q \neq 0 \\ 1 & q = 0 \end{cases}$$

The Dirichlet kernel interpolation also allows us to build the rescaling of the FGF samples. Concretely, we apply equation A1 to interpolate between the original cartesian grid, and another cartesian grid rescaled by a factor s and rotated by an angle $\alpha$. The high frequencies, corresponding to $sk$ larger than the maximum $k$ on the lattice $k_{max}$, are then regenerated by drawing from the FGF distribution.
<figure><img src="image_17.png" /><figcaption>FIG. 16. Diagram of the U-Net architecture used in all experiments. The dark faces of the volumes corresponds to the image or coarse-grained image sizes, while the width represents the number of channels at each step. The right arrows correspond to convolutional layers, while the downwards and upwards ones are pooling and inverse pooling operations. The dashed arrows indicate the skip connections, where the representation of the data at a given scale is reused later on during the upscaling process.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>FIG. 17. Riesz transforms up to second order applied to a sample of the FGF. The orientation corresponding to each filter is clearly visible.</figcaption></figure>

## Appendix B: Experiment details

In this section we detail the experimental setup for both FGF tasks and for the ASM experiments, including details on all the models used and some further results.

### 1. FGF: phase mixing task

a. **U-Net.** We use a 4-layer U-Net with an initial convolutional layer that outputs 8 channels, and a channel distribution of 8 → 16 → 32 → 64 → 128 for the downwards and upwards layers. All convolutions use 3×3 complex-valued filters. A final convolution outputs one single channel from the 8 feature channels at the previous layer. The model is trained for 500 epochs with a learning rate of $10^{-3}$.

b. **Riesz network.** We use a 4-layer Riesz network with channel distribution of 1 → 32 → 40 → 48 → 1. The model is trained for 50 epochs with a learning rate of $10^{-3}$. Note that this model is trained for a shorter time due to the higher computational cost.

In Fig. 18 we show the train and validation loss for each of the three models. Fig. 19 shows a sample output of the three models and the corresponding absolute error. While the output of the Fourier embedding matches the target, there is a visible loss of the high-resolution details for both the Riesz network and the U-net, which matches the quantitative result stated in the main text, that they fail to overcome the spectral bias problem.
<figure><img src="image_19.png" /><figcaption>FIG. 18. Train and validation loss for the U-Net, Riesz network, and Fourier embedding network respectively in the FGF phase mixing task, with no extrapolation.</figcaption></figure>

## 2. FGF: spectral flow task

For the spectral flow task, we trained a U-Net and a Riesz network with the same parameters as the one used for the phase mixing task. Note that although in this case the dataset consists of only 250 samples, this is still enough to learn, as each sample contains 64 × 64 pixels. In Fig. 20 we show the train and validation loss for the three models (U-Net, Riesz network and FM network). We observe severe overfitting for the first two models. Even with our use of early stopping, the validation error is still very important.

Fig. 21 shows, for the three models, a sample output along with its absolute error. The Fourier-Mellin (FM) has the best performance out of the three, with the difference being concentrated around the edges and at the smaller scales, while the U-Net and Riesz network only manage to capture some of large scale structure.

Finally, in Fig. 22, for the FM network, we consider the normalized Mean Squared Error $E(|k|)$ (equation (48))

<figure><img src="image_20.png" /><figcaption>FIG. 19. Sample output of the Riesz network, U-net and Fourier embedding network on the FGF.</figcaption></figure>

in Fourier space as we vary the rescaling factor. As expected, we observe a decrease in performance as the rescaling becomes stronger (smaller $s$). Note that the noise region becomes larger as the rescaling factor becomes smaller, decreasing the effective amount of data available to the model in the learning process, and extending the region $k \in [sk_{max}, k_{max}]$ where it is impossible to perform better than random (since at these frequencies the ground truth target is generated and independent from the input).

# Appendix C: Scale equivariance

## 1. General properties

### a. Reminders

In the context of Neural Networks, the notion of (deterministic) scale equivariance (applied to images) has been introduced in [10, 11]. Here we aim to address the case of statistically self-similar fields and self-similar processes, and thus introduce the appropriate definitions, using less formal notations.
<figure><img src="image_21.png" /><figcaption>FIG. 20. Train and validation loss for the U-Net, Riesz network, and FM network in the FGF spectral flow task, with no extrapolation. Note that the loss is not normalized, and the first two correspond to the MSE loss in direct space, while the third one is the MSE loss in Fourier-Mellin space.</figcaption></figure>

For simplicity, we focus solely on dilations centered at the origin, excluding translations and rotations. Let us denote $T_s \in G$ an element of the group [46] of scale transformations $G$ with $s \in \mathbb{R}^+$ the scaling factor associated to $T_s$. For points or vectors $\mathbf{x} \in \mathbb{R}^d$, the element $T_s$ acts on $\mathbf{x}$ as follows:

$$T_s \cdot \mathbf{x} = s\mathbf{x} \qquad (C1)$$

For some field $\varphi : X \to \mathbb{R}$ with $X \subset \mathbb{R}^d$, we have a straightforward representation of the group action given by

$$T_s \cdot (\varphi(\mathbf{x})) \stackrel{\text{def}}{=} \varphi(s\mathbf{x}) = \varphi(T_s \cdot \mathbf{x}) \qquad (C2)$$

$$\text{i.e.} \quad T_s \cdot (\varphi(\cdot)) = \varphi(T_s \cdot \cdot) \qquad (C3)$$

Taking the Fourier transform of $\varphi(s\mathbf{x})$ we obtain the known property of the time scaling (in 1D), or more gen-

<figure><img src="image_22.png" /><figcaption>FIG. 21. Sample output of the Riesz network, U-net and FM network on the FGF spectral flow task. The absolute error with respect to the ground truth is shown with equal color scale between the 3 models to allow for comparison.</figcaption></figure>

erally in $d$ dimensions:

$$(T_s \cdot \tilde{\varphi})(\mathbf{k}) \stackrel{\text{def}}{=} \tilde{\varphi}_s(\mathbf{k}) = \frac{1}{s^d} \tilde{\varphi}\left(\frac{\mathbf{k}}{s}\right) \qquad (C4)$$

Now, in addition to these generic facts, we may assume our fields and/or processes to respect some statistical equivariance conditions. We recall that the field $\varphi$ is (statistically) self-similar with Hurst exponent $H \in \mathbb{R}$ iff:

$$\forall \mathbf{x}, \quad T_s \cdot (\varphi(\mathbf{x})) \stackrel{\text{def}}{=} \varphi(s\mathbf{x}) \stackrel{d}{=} s^H \varphi(\mathbf{x}) \qquad (C5)$$

where the equality is in distribution. This expression is equivalent to that of the covariance of the increments between points very close together:

$$\mathbb{E}_{\varphi} [ |\varphi(\mathbf{x} + \delta\mathbf{x}) - \varphi(\mathbf{x})|^2 ] \sim |\delta\mathbf{x}|^{2H}, \quad \text{for } |\delta\mathbf{x}| \sim 0$$

In Fourier, this scaling relation corresponds to the high-frequency limit of the energy spectral density (often simply called power spectrum in physics):

$$\mathbb{E}_{\varphi} [ |\tilde{\varphi}(\mathbf{k})|^2 ] \propto |\mathbf{k}|^{-(2H+d)}$$
Error vs. rescaling factor

<figure><img src="image_23.png" /><figcaption>FIG. 22. Normalized Mean Squared Error $E(|\mathbf{k}|)$ (equation (48)) for the FM network on the spectral flow FGF task, varying the rescaling factor between $s = 0.4$ and $s = 0.95$, and a rotation angle of $\pi/4$. In grey is the region corresponding to the noise regime for a scaling factor of $s = 0.9$ (the one used in the experiments shown in the main text). This region extends to the left as $s$ decreases. The white line at $E(|\mathbf{k}|) = 1 = 10^0$ corresponds to the constant predictor $\hat{\psi} = 0$.</figcaption></figure>

Note that for stationary stochastic processes, in general $\varphi$ is not square integrable and there is a nuance between energy density and power density (evaluated over a finite window). Here in practice we will always work on a finite grid with finite resolution, so the difference amounts to a constant prefactor, that we do not care about.

Concretely, for the FGF with Hurst exponent $H = (\beta - d)/2$, in Fourier and in $d = 2$ dimensions, the components $\varphi(\mathbf{k})$ are Gaussians with standard deviation $k^{-\frac{2H+d}{2}} = k^{-\beta/2}$, hence, taking the Fourier transform on both sides of (C5) and applying (C4) on the left, we have:

$$ \tilde{\varphi}_s = \frac{1}{s^d} \tilde{\varphi} \left( \frac{\mathbf{k}}{s} \right) \stackrel{d}{=} s^H \tilde{\varphi}(\mathbf{k}) \quad (C6) $$

which is equivalent to (8). It's worth noting that the scaling of the Fourier transform of the rescaled field is thus: $\tilde{\varphi}_s(\mathbf{k}) \stackrel{d}{=} s^H \tilde{\varphi}(\mathbf{k})$.

The Boltzmann distribution $P(\varphi)$ ought to be conserved under re-scaling, as we assume $\varphi$ to be self-similar. In particular, the stochastic integral $\int d^2\mathbf{k} |\mathbf{k}|^\beta |\tilde{\varphi}(\mathbf{k})|^2$ should be invariant under the change of variable $\mathbf{k}' = \mathbf{k}/s$, which is non trivial for a stochastic integral (since the scale at which we observe the field changes its scaling, an effect that competes with the Jacobian of the change of variable). One can indeed check that this property holds for the integral $I(\varphi) = \int d^2\mathbf{k} |\mathbf{k}|^\alpha |\tilde{\varphi}(\mathbf{k})|^2$ iff $\alpha = \beta$, using equation 8.

## 2. Scale equivariance for linear operators

For linear operators, imposing scale equivariance in addition to rotation invariance leads to a simple functional form which could be used as a basic building block of deep architectures, when combined with ReLU activation functions for instance. Let us recall here a standard argument used to obtain the form of scaling function. Considering the kernel $F$ defining the functional map (49), assuming rotational invariance imposes the form

$$ F(\mathbf{k}, \mathbf{k}') = F(k, k', \alpha - \alpha') $$

Rewriting (49) for $k \to k/s$ and making the change of variable $k' \to k'/s$ in the integral we get

$$ \psi(\mathbf{k}/s) = s^{-2} \int d^2\mathbf{k}' F_{\theta}(k/s, k'/s, \alpha - \alpha') \varphi(\mathbf{k}'/s) \quad (C7) $$

By definitions (50,51), we have $\varphi_s \sim \text{FGF}(\beta)$ and $\psi_s \sim \text{FGF}(\gamma)$ and these should therefore be related by the exact same relation (49) if this one is equivariant, so we have

$$ \begin{aligned} \psi_s(\mathbf{k}) &= \int d^2\mathbf{k} F(k, k', \alpha - \alpha') \varphi_s(\mathbf{k}') \\ &= s^{(\beta-\gamma)/2-2} \int d^2\mathbf{k}' F(k/s, k'/s, \alpha - \alpha') \varphi_s(\mathbf{k}') \end{aligned} $$

using (C7) in the second line. Since this is supposed to be valid for any field $\varphi \sim \text{FGF}(\beta)$ the model should verify

$$ F(k, k', \alpha - \alpha') = s^{(\beta-\gamma)/2-2} F(k/s, k'/s, \alpha - \alpha') \quad (C8) $$

for any $k, k'$ and $s$. Hence choosing $s = k$ leads to the general form given in (52) after symmetrizing by convention in $k$ and $k'$.

## 3. Fourier-Mellin basis

The decomposition of the operator $F = D_F^{-\gamma/2} D_{FM} D_F^{\beta/2}$ involves a the strictly scale invariant operator $D_{FM}$ which becomes diagonal in the “Fourier-Mellin” representation. This representation corresponds to use the Fourier-Mellin basis (53). When $F$ reduces to $D_{FM}$ we have the decomposition

$$ F(\mathbf{k}, \mathbf{k}') = \sum_{\mu \in \mathbb{Z}} \int d\lambda D(\lambda, \mu) \phi_{\lambda, \mu}(\mathbf{k}) \phi_{\lambda, \mu}(\mathbf{k}') $$

with

$$ D(\lambda, \mu) = \frac{1}{\sqrt{4\pi \log(L)}} \int_{-\log(L)}^{\log(L)} du \int_{0}^{2\pi} f(e^u, \alpha) e^{i(\lambda u + \mu \alpha)} $$

valid in the limit where $\log(L) \to \infty$. For finite $L$ some numerical problems can be encountered when using this basis in practice, as orthonormality of the Fourier-Mellin basis is strictly ensured with the definition (53), which is the reason why we have to resort to a log polar lattice interpolated from the square with a Dirichlet kernel.
# Appendix D: Asymptotic behavior of the linear model with RMT

## 1. Train and test error as function of time

In order to write the train and test error along the learning process it is convenient first to introduce the $M \times L^4$ feature matrix $J$ defined by its matrix elements indexed by $p$ and implicitly by $(\mathbf{k}, \mathbf{k}')$ (which can take $L^4$ possible values):

$$J_p(\mathbf{k}, \mathbf{k}') \stackrel{\text{def}}{=} \frac{1}{k'^{\beta/2}} f_p(\mathbf{k}, \mathbf{k}').$$

Given that, the population matrix (57) and the coefficients (56) rewrite

$$C = JJ^{\top} \quad \text{and} \quad w = J^{\top}\tilde{K}$$

where $\tilde{K}$ is the $L^4$ dimensional vector defined as

$$\tilde{K}(\mathbf{k}, \mathbf{k}') \stackrel{\text{def}}{=} \frac{1}{k'^{\beta/2}} K^{\parallel}(\mathbf{k}, \mathbf{k}').$$

In terms of all these quantities the train and test error take the following form as a function of learning time:

$$E_{\text{train}}(t) = \text{Tr} \left[ \hat{C}(\mathbb{I} - \hat{G}(t)\hat{C})^2 J^{\top\dagger} \tilde{K} \tilde{K}^{\top} J^{\dagger} \right] \\ + \sum_{k} \sigma_{\text{eff}}(k)^2 \left[ 1 - \frac{2}{N} \text{Tr}(\hat{G}(t)\hat{C}(k)) \right. \\ \left. + \frac{1}{N} \text{Tr}(\hat{G}(t)\hat{C}\hat{G}(t)\hat{C}(k)) \right] \quad (\text{D1})$$

$$E_{\text{test}}(t) = \text{Tr} \left[ (\mathbb{I} - \hat{C}\hat{G}(t)) C (\mathbb{I} - \hat{C}\hat{G}(t)) J^{\top\dagger} \tilde{K} \tilde{K}^{\top} J^{\dagger} \right] \\ + \sum_{k} \sigma_{\text{eff}}(k)^2 \left[ 1 + \frac{1}{N} \text{Tr}(\hat{G}(t)\hat{C}(k)\hat{G}(t)C) \right] \quad (\text{D2})$$

with

$$\hat{C}_{pq}(k) = \sum_{\mathbf{k}', \mathbf{k}'' \in \Omega^2} f_p(\mathbf{k}, \mathbf{k}') f_q(\mathbf{k}, \mathbf{k}'') \mathbb{E}_{\varphi \sim \widehat{\text{FGF}}(\beta)} [\varphi(\mathbf{k}') \varphi(\mathbf{k}'')]$$

(by rotation invariance this actually depends only on the modulus of $\mathbf{k}$) and where $\sigma_{\text{eff}}^2$ is the effective variance of the noise which emerges from the fact that the regression is misspecified. As mentioned in the plain text this noise comes from the contribution of $K^{\perp}$ in (58).

## 2. Dynamics of train and test errors

Train and test errors given by (D1,D2) can be analyzed asymptotically in the so-called proportional regime, namely when $N$, $N_s$ and $P$ go to infinity at fixed $\rho$ in 65. To this end some preparatory work needs to be done.

First, the longitudinal contribution represented by the bias term in (D1,D2), may be rewritten in a suitable form allowing in the asymptotic limit. Since the features are not learned, $\tilde{K}$ is not correlated with $\hat{C}$ neither aligned with $C$. As a result one of the contribution to the bias term reads e.g. in the train error:

$$\begin{aligned} \text{Tr} [\hat{G}(t) \hat{C} \hat{G}(t) J^{\top\dagger} \tilde{K} \tilde{K}^{\top} J^{\dagger}] & \approx \frac{\kappa}{M} \text{Tr} [\hat{G}(t) \hat{C} \hat{G}(t) J^{\top\dagger} J^{\dagger}] \\ &= \frac{\kappa}{M} \text{Tr} [\hat{G}(t) \hat{C} \hat{G}(t) C^{-1}] \end{aligned}$$

where $\approx$ stands for $O(\frac{1}{\sqrt{M}})$ errors. As can be noticed, a $k$ dependent term shows up in the variance term of both the train and test errors (D1,D2). This term renders the asymptotic analysis much more involved as will be detailed in the following. But this actually simplifies when considering a different setting of the learning where a preconditioner is used to flatten the contributions of all scales

$$\mathcal{L}(\theta) = \frac{1}{2} \mathbb{E}_{(\varphi, \psi) \sim \widehat{\text{FGF}}(\beta, \gamma)} \left[ \sum_{\mathbf{k} \in \Omega} k^{\gamma} |\psi_{\theta}[\varphi](\mathbf{k}) - \psi(\mathbf{k})|^2 \right],$$

and assume that the rescaled noise $k^{\gamma}\sigma_{\text{eff}}(k)^2 = \bar{\sigma}_{\text{eff}}$ is constant. Doing so requires to change accordingly the definitions of $\hat{C}$ and $\hat{Z}$ in (62,63) by an extra factor $k^{\gamma}$ inside the sum over $\mathbf{k}$. In turn the population matrix is now defined as

$$C_{pq} = \langle f_p, f_q \rangle_{\beta, \gamma}$$

corresponding to a change in the definition (55) of the scalar product between kernels to

$$\langle F_1, F_2 \rangle_{\beta, \gamma} \stackrel{\text{def}}{=} \int d^2 \mathbf{k} d^2 \mathbf{k}' \frac{k^{\gamma}}{k'^{\beta}} F_1(\mathbf{k}, \mathbf{k}') F_2(\mathbf{k}, \mathbf{k}')$$

and accordingly in the definitions of $\tilde{K}$ and $J$ above

$$J_p(\mathbf{k}, \mathbf{k}') \stackrel{\text{def}}{=} \frac{k^{\gamma/2}}{k'^{\beta/2}} f_p(\mathbf{k}, \mathbf{k}'),$$

$$\tilde{K}(\mathbf{k}, \mathbf{k}') \stackrel{\text{def}}{=} \frac{k^{\gamma/2}}{k'^{\beta/2}} K^{\parallel}(\mathbf{k}, \mathbf{k}'),$$

and finally the new definition for $\kappa$:

$$\kappa \stackrel{\text{def}}{=} \langle K^{\parallel}, K^{\parallel} \rangle_{\beta, \gamma}.$$

With all that, the train and test errors now takes the form

$$E_{\text{train}}(t) \approx \kappa \text{Tr} \left[ \hat{C}(\mathbb{I} - \hat{G}(t)\hat{C})^2 C^{-1} \right] \\ + \sigma_{\text{eff}}^2 \left[ 1 - \frac{1}{\rho} + \frac{1}{N} \text{Tr}(\mathbb{I} - \hat{G}(t)\hat{C})^2 \right] \quad (\text{D3})$$

$$E_{\text{test}}(t) \approx \kappa \text{Tr} \left[ (\mathbb{I} - \hat{C}\hat{G}(t)) C (\mathbb{I} - \hat{C}\hat{G}(t)) C^{-1} \right] \\ + \sigma_{\text{eff}}^2 \left[ 1 + \frac{1}{N} \text{Tr}(\hat{G}(t)\hat{C}\hat{G}(t)C) \right] \quad (\text{D4})$$
which can be straightforwardly analyzed asymptotically, at least when $C = \mathbb{I}$. Let us call

$$j(x,t) = \frac{1 - e^{-xt}}{x}$$

the function governing the dynamics of the vector parameter $\theta$. Using the Cauchy integral representation we may now express the train and test errors in terms of the resolvent

$$\hat{G}(z) = \frac{1}{z\mathbb{I} - \hat{C}}.$$

We have

$$E_{\text{train}}(t) \approx \frac{\kappa}{2i\pi} \oint_C dz z[1 - zj(z,t)]^2 \text{Tr}[\hat{G}(z)C^{-1}] \\ + \sigma_{\text{eff}}^2 (1 - \rho^{-1}) \\ + \frac{1}{2i\pi N} \oint_C dz[1 - zj(z,t)]^2 \text{Tr}[\hat{G}(z)] \quad (\text{D5})$$

$$E_{\text{test}}(t) \approx \frac{\kappa}{2i\pi} \oint_C dz_1[1 - z_1j(z_1,t)] \oint_C dz_2[1 - z_2j(z_2,t)] \\ \times \text{Tr}[\hat{G}(z_1)C\hat{G}(z_2)C^{-1}] \\ + \sigma_{\text{eff}}^2 \left(1 + \frac{1}{2i\pi N} \oint_C dz z j(z,t)^2 \text{Tr}[\hat{G}(z)C]\right) \quad (\text{D6})$$

and the contour $C$ is an anti-clockwise contour around the real positive axis. In this form we can obtain asymptotic expressions thanks to basic random matrix theory (RMT) formulae [47, 48] at least for the train error and for the test error when $C = \mathbb{I}$. For $C \neq \mathbb{I}$ the bias term of the test error requires more advances techniques to be taken care of, and we postpone this discussion for the moment and concentrate on the case $C = \mathbb{I}$.

### 3. Asymptotic limits of learning trajectories

Let us first derive the RMT equations in this context. This is a little bit non-standard because the covariant matrix takes the form of a Wishart matrix with non-iid samples. Such case has been addressed in [49] for instance. We have

$$\hat{C} = \frac{1}{N} \sum_{s=1, \mathbf{k} \in \Omega}^{N} \mathbf{x}^{(s)}(\mathbf{k}) \mathbf{x}^{(s)\top}(\mathbf{k}),$$

where

$$x_p^{(s)}(\mathbf{k}) = \sum_{\mathbf{k}' \in \Omega} f_p(\mathbf{k}, \mathbf{k}') \varphi^{(s)}(\mathbf{k}') , \quad p = 1, \dots, M, \quad s = 1, \dots, N. \quad (\text{D7})$$

It is a sum of $NL^2$ rank one matrices formed with the vectors $\mathbf{x}^{(s)}(\mathbf{k})$ indexed by $s$ and $\mathbf{k}$. These are not iid because of the $\mathbf{k}$ dependency. The essence of the method resides in finding if there exists a deterministic equivalent [50] $G(z)$ of the resolvent $\hat{G}(z)$. This would mean that for any sequence of $\mathbb{R}^M \times \mathbb{R}^M$ matrices $Q_M$ of Frobenius norm equal to one, $\text{Tr}[(G(z) - \hat{G}(z))Q_M]$ goes to zero typically like $\sqrt{\log(N)/N}$ (see e.g. [51]). If it exists, $G(z)$ should coincide with the average of $\hat{G}(z)$ when taking the proportional limit. As a 2-points function, this obeys a Dyson equation

$$G(z) = \frac{1}{z} + \frac{1}{z} \Sigma(z) G(z)$$

where $\Sigma(z)$ corresponds to the self-energy which in the context of RMT reduces to a sum over planar diagrams. Let us re-derive this informally with the Feynmann diagrams expansion [43, 52] when both the features and the data are correlated. This corresponds to a data matrix of the form

$$X = A^{\frac{1}{2}} Z B^{\frac{1}{2}}, \quad (\text{D8})$$

where $X$ is $M \times N$ matrix, $A$ is $M^2$ while $B$ is $N^2$, representing respectively the feature and data linear dependency. $Z$ is a $M \times N$ random matrix with iid elements. We are interested by the limit spectrum of

$$\hat{C} \stackrel{\text{def}}{=} \frac{1}{N} X X^{\top} \quad \text{and} \quad \hat{C}^{\star} \stackrel{\text{def}}{=} \frac{1}{N} X^{\top} X$$

Note that the relation between $A$, $B$ and the population matrix is

$$C = A \lim_{N \to \infty} \frac{1}{N} \text{Tr}[B].$$

Let us consider the resolvent corresponding to the $(N+M)^2$ matrix

$$M(z) \stackrel{\text{def}}{=} \begin{pmatrix} \sqrt{z} \mathbb{I}_M & -\frac{X}{\sqrt{N}} \\ -\frac{X^\top}{\sqrt{N}} & \sqrt{z} \mathbb{I}_N \end{pmatrix} = \sqrt{z} \mathbb{I}_{N+M} - \frac{\hat{S}}{\sqrt{N}}$$

with $\hat{S}$ denoting now the data matrix written in this two-component representation. Its inverse is given by

$$M(z)^{-1} = \sqrt{z} \begin{pmatrix} \hat{G}(z) & -\frac{1}{\sqrt{z}} \hat{G}(z) X \\ -\frac{1}{\sqrt{z}} \hat{G}^{\star}(z) X^{\top} & \hat{G}^{\star}(z) \end{pmatrix} \quad (\text{D9})$$

with

$$\hat{G}(z) = (z \mathbb{I}_M - \hat{C})^{-1}$$

$$\hat{G}^{\star}(z) = (z \mathbb{I}_N - \hat{C}^{\star})^{-1}$$

Expanding $M^{-1}$ in terms of the data matrix and averaging over $Z$, leads, thanks to the Wick theorem (see
e.g. [52]) to obtain the self energy $\Sigma(z)$ of the Dyson equation

$$G(z) = G_0(z) + G_0(z)\Sigma(z)G(z)$$

concerning the resolvent $G(z)$ of $M(z)$ in the asymptotic limit and where $G_0(z) = \frac{\mathbb{I}_{N+M}}{\sqrt{z}}$. Thanks to the planar diagrams approximation we get for the self energy

$$\Sigma(z) = \lim_{N,M \to \infty} \frac{1}{N} \mathbb{E}[\hat{S}G(z)\hat{S}] \\ = \lim_{N,M \to \infty} \begin{pmatrix} \frac{1}{N} \mathrm{Tr}[G^*(z)B]A & 0 \\ 0 & \frac{1}{N} \mathrm{Tr}[G(z)A]B \end{pmatrix}$$

where $G(z)$ and $G^*(z)$ are the deterministic equivalents respectively for the resolvent for $\hat{C}$ and $\hat{C}^*$. Since we know from (D9) that

$$G(z) = \sqrt{z} \begin{pmatrix} G(z) & 0 \\ 0 & G^*(z) \end{pmatrix}$$

we end up with the RMT equations

$$\Gamma(z) = \frac{1}{\rho} \int \nu(dx) \frac{x}{z - x\Gamma^*(z)} \\ \Gamma^*(z) = \int \nu^*(dx) \frac{x}{z - x\Gamma(z)}$$

where

$$\Gamma(z) = \lim_{N,M \to \infty} \frac{1}{N} \mathrm{Tr}[\hat{G}^*(z)A] \\ \Gamma^*(z) = \lim_{N,M \to \infty} \frac{1}{N} \mathrm{Tr}[\hat{G}(z)B]$$

and $\nu$ and $\nu^*$ are the spectral density respectively associated to A and B. If we specify this now to our problem, X is to be seen as a $M \times (NL^2)$ matrix with elements $x_p^{(s)}(\mathbf{k})$, $p = 1, \dots, M$ and $(s, \mathbf{k}) \in \{1, \dots, N\} \times \Omega$. From (D7) the dependency structure is in general more entangled than (D8) and leads to RMT equations where $\Gamma(z)$ and $\Gamma^*(z)$ have matrix form of respective sizes $M \times M$ and $L^2 \times L^2$, which cannot be given a sense in the asymptotic limit. In order to recover finite dimensional

mean-field equations we need to impose that dependency factorizes in the following way:

$$\mathbb{E}[x_p^{(s)}(\mathbf{k})x_q^{(s')}(\mathbf{k}')] = A_{pq}B(\mathbf{k}, \mathbf{k}')\mathbb{1}_{\{s=s'\}}$$

which corresponds to assume that the features obey the following relations

$$\sum_{\mathbf{k} \in \Omega} \frac{f_p(\mathbf{k}, \mathbf{k}_1) f_q(\mathbf{k}', \mathbf{k}_1)}{k_1^{\beta}} = A_{p,q} B(\mathbf{k}, \mathbf{k}')$$

This condition is obeyed for instance when the feature are rank one operators based on a set of orthonormal functions w.r.t. to the inner product associated to the density $\mu(d\mathbf{k}) = k^\beta$. In that case we have $A = \mathbb{I}$.

Along the contour of the Cauchy integrals in (D5,D6) we integrate over $z = y + i\epsilon$ with $\epsilon$ infinitesimal. In the limit $\epsilon \to 0$, both $\Gamma$ and $\Gamma^*$ may acquire a finite imaginary part which we write as

$$\lim_{\epsilon \to 0^\pm} \Gamma(z) = \Gamma_r(y) \pm \Gamma_i(y) \\ \lim_{\epsilon \to 0^\pm} \Gamma^*(z) = \Gamma_r^*(y) \pm \Gamma_i^*(y).$$

Letting $\bar{\nu}(dx)$ the asymptotic limit of the empirical spectrum in the proportional regime, its Stieltjes transform is given by the trace of the resolvent:

$$g(z) \stackrel{\text{def}}{=} \int \frac{\bar{\nu}(dx)}{z - x}.$$

Then the spectrum is given by the Stieltjes transform

$$g(y + i\epsilon) = g_r(y) + i\pi \frac{\epsilon}{|\epsilon|} \bar{\nu}(y),$$

which thanks to the identity

$$\Gamma(z)\Gamma^*(z) = \frac{1}{\rho}[1 + zg(z)]$$

rewrites

$$\bar{\nu}(y) = (1 - \rho)\mathbb{1}_{\{\rho < 1\}}\delta(y) + \rho \frac{\Gamma_i(y)\Gamma_r^*(y) + \Gamma_r(y)\Gamma_i^*(y)}{\pi y} \quad (D10)$$

In terms of these quantities we obtain equations (67,68) for the train and test errors up to some quite technically involved additional corrections when $C \neq \mathbb{I}$ to evaluate exactly the bias term of the test error, i.e. the first term of (D2).