# Convolutions on partially regular recurrent lattices

Ö.D. Gürcan, L. Manfredini

► To cite this version:

Ö.D. Gürcan, L. Manfredini. Convolutions on partially regular recurrent lattices. Communications in Nonlinear Science and Numerical Simulation, 2025, 152, pp.109262. 10.1016/j.cnsns.2025.109262. hal-05365037

HAL Id: hal-05365037

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

Submitted on 14 Nov 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Research paper

# Convolutions on partially regular recurrent lattices

Ö.D. Gürcan<sup>a,∗</sup>, L. Manfredini<sup>a</sup>

<sup>a</sup> Laboratoire de Physique des Plasmas, CNRS, Ecole Polytechnique, Sorbonne Université, Université Paris-Saclay, Observatoire de Paris, F-91120, Palaiseau, France

## ARTICLE INFO

**Keywords:**

Fast Fourier transform

Log-lattices

Convolution

Turbulence

## ABSTRACT

Partially regular recurrent lattices, are k-space grids that we propose, which consist of a central regular region that is extended using a recurrence relation, resulting in an asymptotically logarithmic lattice structure. Such a lattice can be used to model the turbulent cascade over a very large range of scales covered by the recurrent part, while keeping the large scale eddies mainly in the regular part. Here we propose a novel pseudo-spectral algorithm for computing the convolutions over such a lattice, using an overlapping partition of its different parts, using the fact that the interactions in the recurrent part of the lattice are limited, in each direction, to a small number of elements linked through the recurrence relation. We compare the results with a full grid, dense, fast Fourier transform (fft) based convolution, where the nonexistent elements on the full grid are set to zero, and show that the difference remains within a few orders of the machine precision. The algorithm in two dimensions uses one fft for the regular grid, and bunch of smaller ffs for the rest of the points, either elongated to match the length of the regular part of the grid in one dimension or even smaller ffs (typically 36×36) to include the interactions between recurrent parts of the lattice. The algorithm can be trivially generalized to arbitrary number of dimensions.

## 1. Introduction

Convolutions are ubiquitous in nonlinear problems in physics, numerical analysis [1] and practical applications such as data analysis [2,3] and deep learning [4], since they describe, in a general sense, interactions on a connected network. In particular, convolutions in k-space arise from a Fourier space description of a quadratic nonlinearity in real space, such as the advective nonlinearity in the study of fluids. One of the key challenges in understanding fluid behavior is that of turbulence [5], which presents itself especially at high Reynolds numbers, as a duality of chaotic evolution and hierarchical self-organization across a large range of scales. In direct numerical simulations of turbulence, especially using spectral or pseudo-spectral methods [6,7], the key nonlinear step is the computation of the convolutions. More generally, since turbulence can be formulated as a network of interacting nodes [8,9] representing wave-numbers in Fourier space, convolutions can naturally be used to describe interactions among these nodes, whether or not the underlying grid is regular.

Various numerical tools from direct numerical simulations, to reduced models such as shell models [10] are used, in order to model turbulent cascades [11]. Recently, a particular kind of reduced model, called logarithmic lattice have attracted growing interest [12–15], which aims to bridge the gap between these limiting approaches: Log lattices are based on shell model-like decompositions in each axis, so that for example a three dimensional log-lattice model based on the GOY model (Gledzer-Ohkitani-Yamada) for local interactions can be described as GOY^3, usually with a particular shell spacing, such 2 or the golden ratio.

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: ozgur.gurcan@lpp.polytechnique.fr (Ö.D. Gürcan).</footnote>
Recently, the current authors have shown [16] that a recurrent series, such as the Fibonacci series, which has an asymptotic logarithmic scaling equal to the golden ratio can be used to develop a shell model that is practically equivalent to the GOY model with a golden ratio shell spacing. In this work, it was suggested that these models can be generalized to multiple dimensions resulting in asymptotically logarithmic, recurrent lattices and that the recurrence relation can be applied to a regular grid, resulting in what we call a partially regular recurrent lattice (PRRL), which is basically a lattice that is very similar to a log-lattice, but following integer values and with a regular part in the middle.

We show that starting from the last two elements of a regular grid, and using the recurrence relation $k_n^l = k_{n-1}^l + k_{n-2}^l$, where $i = (x, y, z)$, corresponding to the Fibonacci series, one can construct such a lattice, whose logarithmic scaling factor (i.e. $k_n^l/k_{n-1}^l$) approaches the golden ratio. Note that both the size of the regular grid, and the number of iterations of the recurrence relation can be selected separately in each direction, resulting in very flexible grids. For example a regular grid in one direction and a grid that starts as regular and continues as a recurrent series becoming asymptotically logarithmic in the other direction(s), can be described using a PRRL. However at the simplest level, choosing $k_{0,1}^l = 1, 2$, or $k_{0,1}^l = 1, 3$ we get the Fibonacci and the Lucas sequence recurrent lattice models. These models have the same relation to the log-lattices that the recurrent shell models have to the GOY model with the golden ratio spacing.

Both shell models, and log-lattices have various shortcomings. Probably the most important among those, is that because such models have roughly constant number of Fourier modes as a function of scale, as opposed to a regular grid -that would have $k^d$ modes, where $d$ is the dimensionality- the statistics of these models, for example the equipartition solutions etc. are completely different [17]. Also the fact that a number of modes that fall within a shell is represented by a single complex variable, with a single phase, makes it such that they can not describe spatially coherent structures. Log lattices as well as their close cousins, the nested polyhedra models [18,19], also seems to fail to capture the turbulent intermittency, paralleling what happens with fractal decimation [20]. Nonetheless, these models are capable of describing certain aspects of the turbulent cascade, and both shell models and more recently log-lattices have been used to address many complex physical problems. A turbulence model based on the PRRL, which allows a central region with a regular grid, would be an important improvement over these models, especially in representing coherent structures, somewhat in the spirit of large eddy simulations. However, the recurrent part would have the same shortcomings as shell models and log-lattices.

Note however that, the goal of the current paper, is not to develop such a turbulence model, but to describe how to deal with the convolutions on such a grid. In the following, we propose a very efficient algorithm [21] to compute convolutions on such a lattice, using ffts on overlapping partitions, and noting that the convolution contributions to recurrent wave-numbers can be compressed to a minimal series having the same interaction topology, basically transforming the computation of the contribution to a 36 element fft (including padding). The method is easier to explain in two dimensions but it generalizes trivially to higher dimensions. Note that while we consider mainly convolutions of Hermitian symmetric fields, as one would obtain from real to complex Fourier transforms and therefore treat only half of the $k_y$ direction in a 2D Fourier space, the general case of complex convolution can be trivially implemented by treating $k_y$ the same way as $k_x$.

The compression algorithm can be explained in 1D by considering a purely recurrent grid constructed from the recurrence relation $k_n = k_{n-2} + k_{n-1}$, with $k_0, k_1, k_2 = 0, 1, 3$, which continues following the Lucas numbers $k_n = L_n$. Since the convolution for any $k_n$ on such a recurrent series can be written explicitly as:

$$[gh]_{k_n} = g_0 h_{k_n} + g_{k_n} h_0 + g_{k_{n-2}} h_{k_{n-1}} + g_{k_{n-1}} h_{k_{n-2}} + g_{k_{n-1}}^* h_{k_{n+1}} + g_{k_{n+1}} h_{k_{n-1}}^* + g_{k_{n+1}}^* h_{k_{n+2}} + g_{k_{n+2}} h_{k_{n+1}}^* \quad (1)$$

similar to the nonlinear term of the Sabra model [22], is translation invariant on the Lucas series, we can move the values $g_{k_j}$ and $h_{k_j}$ to the first 6 elements Lucas series (really to $[0, 1, 3, 4, 7, 11]$), compute the convolutions and consider the middle element of the result (i.e. $k_3 = 4$) as the resulting contribution to $k_n$. The method is explained in detail in Section 3.1.

Note that since the largest wave-number of this series excluding the Nyquist wave-number should be 11, we get $N_y = 2 \times (11 + 1) = 24$, and multiplying this by $3/2$ gives us $N_{py} = 36$. This means that we need a padded resolution of 36 in order to contain all the elements involved in this convolution without causing aliasing. Obviously for a purely recurrent series in 1D, there is no obvious advantage to using ffts. However when the number of dimensions increase, the Hermitian symmetric form of Eq. (1) becomes rather complicated, and adding a regular grid region in the middle, especially considering the boundary between the regular and the recurrent grids, becomes practically impossible to compute by hand. The algorithm that we propose solves these problems.

The rest of the paper is organized as follows. Section 2 is dedicated to the general, mathematical formulation of the convolution on a PRRL, where in Section 2.1 we discuss how to construct the PRRL given its two parameters representing the size of the regular grid and the number of iterations, reminding the reader in Section 2.2 the definition of the fft based convolution that we use both for computing the contribution to and from the regular grid part and for the dense fft that we will later use for comparison. Section 2.3 discusses how the convolution can be partitioned and the contributions from each partition combined, both in general and in particular for the PRRL. Section 3 is dedicated to the numerical implementation, with Section 3.1 discussing the computation of compressed contributions for the recurrent partitions, which is at the heart of the algorithm that we propose in Section 3.2. Section 3.3 presents a comparison both in terms of speed and accuracy to a dense fft based algorithm and discusses expected theoretical computational complexity of the algorithm. Section 4 is conclusion.

## 2. PRRL and convolution

### 2.1. Constructing a PRRL

Consider a two dimensional regular grid, in the wave-number space corresponding to a spatial grid of $N_{rx}, N_{ry}$ elements over a box of size $2\pi, 2\pi$ for simplicity. Dropping the Nyquist mode, and using Hermitian symmetry, this gives a Fourier space grid of the form
$k_i^x = [0, 1, 2, ..., N_{rx}/2 - 1, 1 - N_{rx}/2, ..., -2, -1]$ and $k_j^y = [0, 1, 2, ..., N_{ry}/2 - 1]$. Considering $k_i^x$ as the first half of the $k_x$ grid, which has $N_{rx}/2 - 1$ as its last element, and $k_j^y$ as the full $k_y$ grid, we can extend the wave-number domain using the recurrence

$$k_i^j = k_{i-1}^j + k_{i-2}^j \qquad (2)$$

This gives makes

$$k_n^y = \left[ 0, 1, 2, \dots, \frac{N_{ry}}{2} - 1, N_{ry} - 3, \frac{3N_{ry}}{2} - 4, \frac{5N_{ry}}{2} - 7, \dots \right] \qquad (3)$$

considering the $N_{ry}/2 - 2$ and $N_{ry}/2 - 1$ as the first and second elements of a series, we observe that the nth element that follows from this series would have the form $\frac{F_n}{2} N_{ry} - L_{n-1}$, where $F_n$ and $L_n$ are the Fibonacci and the Lucas sequences. The same extension algorithm can be applied to the $k_x$ grid, to extend both negative and positive wave-numbers so that in the end one obtains the grid structure shown in Fig. 1. We denote the number of such iterations as $N_{lx}$ and $N_{ly}$, where the index l is used for logarithmic, since the recurrent part of the lattice is asymptotically logarithmic. Note that a 2D array defined on this lattice would have $N_{kx} = N_{rx} - 1 + 2N_{lx}$ and $N_{ky} = N_{ry}/2 + N_{ly}$ elements, since Nyquist mode is explicitly excluded and $N_l$ really denotes the number of elements added in Fourier space. We take the indices of such an array in two dimensions to be defined from $-N_{rx}/2 - N_{lx} + 1$ to $N_{rx}/2 + N_{lx} - 1$ (inclusive) and from 0 to $N_{ry}/2 + N_{ly} - 1$ using the standard fft indexing scheme. Note that for a real field f, the Hermitian symmetry requires $f_{0,0}$ to be real and $f_{0,-j} = f_{0,j}^*$. This should be imposed explicitly if it is not the case for the initial data. Of course one can apply the algorithm to complex arrays by simply dropping the requirement of Hermitian symmetry and treating the y direction the same way as the x direction. Also everything here can be extended to arbitrary number of dimensions d, and to real fields by treating the last dimension as we treat y, and to complex fields by treating all the dimensions as we treat the x direction. This makes the 2D convolutions of real fields the minimum example that has all the ingredients to be generalized to compute any convolution in any number of dimensions.

Note that using the standard fft indexing, a 2D array of Fourier coefficients on a regular grid are stored conventionally as $f_{i,j} = f_{k_i^x, k_j^y}$. Using the same convention for the PRRL, is quite natural. One interesting side effect is that one can compute the ffts of these arrays directly, which gives a real field representation where the regular part of the grid describes from large scales up to $k_{N_r/2-1}$, and the smaller scales corresponding to the recurrent part is brought to larger scales, as $k_{N_r/2+\ell} = \frac{2\pi}{L_x} \left( \frac{F_{\ell+3}}{2} N_r - L_{\ell+2} \right)$ in each direction is replaced by $\frac{2\pi}{L_x} (N_r/2 + \ell)$. We use this regularly as a quick and dirty way to visualize these fields in real space.

<figure><img src="image_2.png" /><figcaption>Fig. 1. The lattice structure of a PRRL with $N_{r(x,y)} = 16$, $N_{l(x,y)} = 3$. Here $P^{(-1,-1)}$, which is at the bottom center and outlined in blue, is the fully regular part of the grid. The green boxes that are elongated vertically designate the partitions that are regular in y but recurrent in x denoted by $P^{(-n,-1)}$, while the orange boxes that are elongated horizontally correspond to partitions that are regular in x and recurrent in y, denoted by $P^{(-1,n)}$. Note also that $P^{(0,0)}$ is the value at $k_x = k_y = 0$ (i.e. the mean). Some partitions are explicitly noted in the figure (to the top right corner of each corresponding partition) with corresponding colors. Note how the regular grid contains 15 other partitions (basically all the $P^{(n,n)}$ with $n < 3$), while the fully recurrent partitions, such as $P^{(3,3)}$ shown above, consist of single points.</figcaption></figure>
## 2.2. Convolution

Convolution is the Fourier transform of a multiplication, which can be written in general in two dimensions as

$$f_{j,\ell} = \sum_{\ell', j'} g_{j',\ell'} h_{j-j',\ell-\ell'}, \qquad (4)$$

with limits of the sum determining its kind. For example if the Hermitian symmetry is used explicitly to consider $\ell$ to be only over positive integers, the *matched output* convolution sum becomes

$$f_{j\ell} = \sum_{j'=\max\left[1-\frac{N_x}{2}, 1-\frac{N_y}{2}-j\right]}^{\min\left[\frac{N_x}{2}-1-j, \frac{N_y}{2}-1\right]} \left( \sum_{\ell'=0}^{\frac{N_y}{2}-1-\ell} (g_{j',\ell'}^* h_{j+j',\ell+\ell'} + \text{c.c.}) + \sum_{\ell'=0}^{\ell} g_{j'+j,\ell'} h_{-j',\ell-\ell'} \right) \quad (5)$$

for a regular grid of size $N_x$ and $N_y$. The usual way one would compute such a convolution sum is to write it as

$$f_{j\ell} = C_{j\ell}^{\text{fft}}(g, h) = \text{RFFT}(\text{IRFFT}_p(g)\text{IRFFT}_p(h)) \qquad (6)$$

where *IRFFT*<sub>*p*</sub> denotes 3/2 padded complex to real inverse Fourier transforms and *RFFT* denotes real to complex forward Fourier transform followed by extracting the initial (unpadded) part of the Fourier space. This is what we call an fft based convolution in this paper, and apart from an ambiguity in normalization depending on the fft library one uses, the results are equivalent to the explicit convolution sums as written in (5). Note that the fft based convolution is much faster than a direct sum computed by the formula (5) even though the results are basically the same within numerical accuracy (see Section 3.3).

## 2.3. Convolution over a PRRL

As shown in Fig. 1, the points in a PRRL may belong to either the regular or the recurrent part of the grid in each direction. One alternative interpretation of the convolution sum is that it corresponds to a sum over all triadic interactions, as implied for example by the common notation:

$$f_k = \sum_{\Delta} g_{\mathbf{p}}^{*} h_{\mathbf{q}}^{*} \qquad (7)$$

where $\Delta$ denotes a sum over all $\mathbf{p}$ and $\mathbf{q}$ such that $\mathbf{p} + \mathbf{q} = -\mathbf{k}$, where the Hermitian symmetry should be used to replace the $g_{-\mathbf{p}} \to g_{\mathbf{p}}^*$ if $p_y$ or $q_y$ happens to be in the lower half plane. This perspective allows us to consider the convolution sum as being over a list of “triads”. For example when we compute the convolution over the regular part of the grid, we actually consider contributions of all the triads for which all the wave-numbers are in the regular part of the grid. If $\Delta'$ is the rest of the triads (i.e triads for which at least one leg is outside the regular part), we can write the convolution sum as

$$f_{j\ell} = C_{j\ell}^{\text{fft}}(g_{\text{reg}}, h_{\text{reg}}) + \sum_{\Delta'} g_{\mathbf{p}}^{*} h_{\mathbf{q}}^{*}, \qquad (8)$$

where subscript reg denote belonging to the regular part of the grid. However since one can still have two regular legs for the triads in the set $\Delta'$, generating a list of all such triads is still numerically challenging.

Instead we can consider an overlapping partitioning of the k-space, which in 1D can be written as:

$$P^{(n)} = \begin{cases} k_j \in K_{\text{reg}} & \text{if } n = -1 \\ 0 & \text{if } n = 0 \\ \frac{F_n}{2} N_r - L_{n-1} & \text{if } n > 0 \end{cases} \qquad (9)$$

where $K_{\text{reg}} = \{k_j \in \mathbb{Z} : 0 \le k_j < N_r/2\}$ is the 1D regular grid, so that when the partition acts on $g$, it gives the part of $g$ that corresponds to the partition, with the regular grid corresponding to $n = -1$, so that $P^{(-1)}g \to g_{\text{reg}}$. The partition generalizes to multiple dimensions such that $g^{(n,\ell)} \equiv P^{(n,\ell)}g$, by considering $P^{(n)}$ as defined above for the last direction (i.e. $y$ in 2D and $z$ in 3D) but using both positive and negative values $k_\ell = (F_n/2 N_r - L_{n-1}, -F_n/2 N_r + L_{n-1})$ for the $n > 0$ part of the other directions since they are also involved in the same classes of interactions. This gives the partitions shown in Fig. 1.

Since the partitions are overlapping, contributions of each convolution to each combination of partitions is to be considered avoiding double counting:

$$f_{j\ell}^{(n)} = \bigcup_{n'=1}^{N_\ell+2} \bigcup_{n''=-1}^{N_\ell+2} C_{j\ell}^{\text{fft}}(g^{(n')}, h^{(n'')})^{(n)}. \qquad (10)$$

# 3. Numerical implementation

## 3.1. Compressed contributions on recurrent partitions

On the recurrent part of the grid, the interactions are limited to $n-2, n-1, n, n+1, n+2$, and 0 in each direction. Since the last two elements of the regular grid in each direction is also considered as being part of the recurrent series that constitute the lattice, they
<figure><img src="image_3.png" /><figcaption>Fig. 2. Lattice compression in two dimensions for computing $P^{(7,7)}$, shown here for a PRRL with $N_{r(x,y)} = 4$, $N_{l(x,y)} = 8$. Note that $P^{(7,7)}g = g^{(7,7)} = g_{\pm 6,6}$ which corresponds to the wave-numbers $k_x, k_y = \pm 13, 13$. The contribution to these two nodes from the convolutions shown above are identical, since the interactions which involve these wave-numbers are exactly the same. This means that even for $k_n$ values which can be very large, we can copy the corresponding elements to the compressed lattice shown below and compute an fft based convolution, and then extract the two elements corresponding to $P^{(n,n)}$. This computes all the triad interactions with the wave-number $k_{\pm nx}, k_{\pm ny}$ of the lattice, including the interactions with the $k_x = 0$ and $k_y = 0$ axes as well as $k_{00}$.</figcaption></figure>

include the interactions between the regular grid and the recurrent series. The key idea in this paper that simplifies the computation significantly is that we can contract the ffts to a minimal grid using the self-similarity implied by the recurrence relation.

Consider the 1D case, and the consecutive wave-numbers given by the series:

$$k_n = \frac{F_n}{2} N_r - L_{n-1} \qquad (11)$$

with possibly very large values of $n$ and $N_r$. Assume for the time being that $n \ge 3$, so that the intersection between the regular grid and the partition $n$ is null. We can then express the convolution as:

$$f_j^{(n)} = C_j \left( \left[ g^{(0)} + \sum_{n'=n-2}^{n+2} g^{(n')} \right], \left[ h^{(0)} + \sum_{n''=n-2}^{n+2} h^{(n'')} \right] \right)^{(n)} \qquad (12)$$

In order to compute this convolution using ffts we would basically set all the wave-numbers apart from 0, $k_{n-2}$, $k_{n-1}$, $k_n$, $k_{n+1}$, $k_{n+2}$ to zero for both $g$ and $h$, and compute the resulting convolutions using Eq. (6), with padding etc. However for large $k_n$ values this requires computing large ffts with a lot of zeros.

Instead of using the $k_n$'s we can also use the smallest series that satisfy the same recurrence relation as the $k_n$'s as long as this doesn't introduce any additional interactions (see Fig. 2). In other words, by mapping the wave-numbers $k_n$ to a list of small integers, with the same interaction topology as the list of $k_n$'s, we get the same convolution. In order to find this minimum series we consider
<caption>**Table 1**</caption>

input, mapping, valid and output indices for computing the convolution of a recurrent 1D Hermitian symmetric field and the corresponding wave-numbers for the input and output indices. The values of the input arrays at the *input* indices are copied to minimum padded temporary arrays at the *mapping* indices (see algorithm 1 for how this works). The size of these temporary arrays (in real space) has to be 36 in order to contain 11 non-zero independent nodes + padding. After the fft based convolution is computed, it is only the *valid* index (or indices in the more general case), is extracted as the valid part of the resulting contribution of this partition to the total convolution at the *output* index.

<table><thead><tr><th></th><th>input</th><th>mapping</th><th>valid</th><th>output</th></tr></thead><tbody><tr><td>indices</td><td>[0, n - 2, n - 1, n, n + 1, n + 2]</td><td>[0, 1, 3, 4, 7, 11]</td><td>4</td><td>n</td></tr><tr><td>wave-numbers</td><td>[0, k<sub>n-2</sub>, k<sub>n-1</sub>, k<sub>n</sub>, k<sub>n+1</sub>, k<sub>n+2</sub>]</td><td>-</td><td>-</td><td>k<sub>n</sub></td></tr></tbody></table>

the first few Lucas numbers, with a 0 instead of 2 for the first element, giving the series [0, 1, 3, 4, 7, 11]. This means that we can map [0, k<sub>n-2</sub>, k<sub>n-1</sub>, k<sub>n</sub>, k<sub>n+1</sub>, k<sub>n+2</sub>] → [0, 1, 3, 4, 7, 11], and so use a padded fft of size $N_p = 36$ (e.g. $N_p/2 + 1 \ge (11+1) \times 3/2$) in order to compute the contribution to $f^{(n)}$, regardless of the actual value of $k_n$ (see Table 1). In multiple directions, one has partitions that are regular in one direction and recurrent in the other. In this case, one contracts the recurrent direction while keeping the regular direction dense, since the convolution works independently along separate directions. Note that using the first few Fibonacci numbers instead of the Lucas numbers does not work as well, since they would also have a spurious $1+1=2$ interaction in addition to the interactions arising from the recurrence relations. One could however use [0, 2, 3, 5, 8, 13], but that would require a slightly larger padded fft, so we conclude that the series we consider above is really the minimal one.

<caption>**Algorithm 1** Convolution algorithm using an overlapping partition list.</caption>

fk = fill(0)
for each partition:
    gkp[mapping] = gk(input)
    hkp[mapping] = hk(input)
    fkp = irfft(rfft(gkp)*rfft(hkp))
    fk(output] += fkp[valid]

## 3.2. The convolution algorithm

The algorithm consists of two stages: the initialization phase, where a partition list is generated, and the execution phase where the convolution is computed using the partition list. The partition list can be computed once at the initialization phase and reused for multiple convolutions on the same PRRL. In order to simplify different steps of the execution, each partition keeps track of four separate indices: i) the *input* indices, which are the indices that contribute to the convolution contribution that we compute ii) the *mapping* indices, which are the indices on the minimum (padded) array to which these input indices are to be mapped, iii) *valid* indices, which are the indices on the minimum padded array that correspond to the partition for which we compute the convolution contribution and iv) *output* indices: which are the indices of the partition to which this contribution is to be added. Each index list contains the negative indices as well apart from the last direction of the fft for Hermitian symmetric fields. Note that one can even allocate and keep the temporary padded arrays with the correct sizes in the partition list in order to avoid repeated allocations.

The input indices, basically correspond to the indices of the $n'$ partitions that appear in Eq. (12), while the output indices correspond to the partition $n$ for which the contribution is being computed. While for $n \ge 3$, Eq. (12) gives the full convolution, for the partitions that intersect with the regular grid, it gives only the contribution to the convolution that comes from the irregular part. In order to exclude double counting, the $n=0$ mode, should be excluded while constructing the input indices for these partitions, since $n=0$ is already part of the regular grid. Similarly when considering the contribution to $n=0$, we need to construct the input indices using all the recurrent partitions, excluding $n=1$ and $n=2$, since the contributions from those are already computed as part of the regular grid. Finally for the last few partitions, we can consider only the necessary part of the minimum series and use slightly shorter Fourier transforms. Once the partition list is constructed the convolution can be computed using the simple algorithm described in pseudo-code in Algorithm 1.

## 3.3. Performance

We consider the PRRL as a Fourier space decimation. Therefore the natural comparison for the convolution algorithm that we propose is to a standard dense fft based convolution where the non-existing part of the lattice is set to zero. The comparison can be seen in Fig. 4 for white noise with random phase initial data for $g$ and $h$, without any power law $k$-dependence. Note that this is one of the more challenging cases, since a falling power law $k$-dependence would decrease the contribution of the high-$k$ part of the spectrum, and result in a reduced error. Note that here we keep the Fourier coefficients for the recurrent wave-numbers as consequent
<figure><img src="image_4.png" /><figcaption>Fig. 3. Lattice compression and the different indices that are needed to compute the partial contribution to the convolution for a) $n = 0$ on the top left b) $n = 1$ on the bottom left c) $n = 2$ on the top right and d) $n > 2$ on the bottom right. Here what is shown explicitly is how we treat $P^{(n-1)}$ for a 2D convolution, where $k_n$ etc. really denote $k_{nx}$. However considering this for a given direction corresponding to $n$, while treating the other directions either as regular or recurrent (hence using this scheme in that direction as well) we get the general algorithm for $n$ dimensions. In each case the list of indices corresponding to the input and mapping are explicitly given. The valid and output indices are denoted in red in the lists of mapping and input indices respectively. Note that the $n = 0$ is excluded in the computation of $n = 1$ and $n = 2$, while those two are in turn excluded from the computation of $n = 0$, since all these three partitions are also part of the regular grid, and the interactions among them are already included in the contribution from the regular grid. Note also that for the regular axis (here y axis) all the indices are the same and goes from 0 to $N_r/2 - 1$ [would be $-(N_r/2 - 1)$ to $(N_r/2 - 1)$ if it wasn't the last direction for a real convolution].</figcaption></figure>

elements in a compact initial array (i.e. $g_{N_r+\ell} = g_{k_{N_r+\ell}}$). One can see the resulting convolution in the bottom left plot of Fig. 4, and the error in the bottom right plot. Note that the error is within two orders of magnitude above machine precision.

In order to perform speed comparisons, we consider the two algorithms above, and since the algorithm is not (yet) parallelized, the comparisons are performed up to sizes that are possible for the dense fft algorithm on a single workstation. The code is written in python for readability, and the performance study is done on a single work station [Intel i9-10900K (10 Core) CPU @ 3.70GHz]. As the difference in speed is substantial for relevant cases, a more detailed study does not seem to be necessary to show the superiority of our algorithm at this level. Maybe an alternative algorithm using sparse matrices may be developed for a proper comparison in the future when the algorithm is also parallelized. This is out of scope of the current paper, which focuses mainly on the mathematical basis behind the algorithm, and not its efficient implementation.

It is true that comparing to dense fft based convolution in terms of speed is somewhat dubious, especially at very high resolutions where the memory access becomes the main bottleneck. One alternative maybe to go back to direct convolution sums. However it becomes rather cumbersome to keep track of interactions in an explicit convolution sum, and unless one separates the regular part as we have done above, the regular grid part of the explicit convolution becomes extremely slow unless the regular grid is very small.

We can estimate the theoretical scaling by considering the number of partitions and the size of each partition in order to compute the number and sizes of ffs that needs to be computed. In a given direction there are $N_l + 4$ partitions ($N_l$ partitions for the recurrent part plus $n = -1, 0, 1, 2$ corresponding to the regular grid, and its boundaries). Among those, there are $N_l + 3$ ffs of size 36 (or less for boundary terms) and one padded fft of size $3N_r/2$ for the regular grid part. Thus, For a $d$ dimensional convolution, the number of degrees of freedom, $N_d$, for the padded fast Fourier transform can be written as:

$$N_d = \left( 36 \times (N_l + 3) + \frac{3}{2} \times N_r \right)^d \quad (13)$$
<figure><img src="image_5.png" /><figcaption>Fig. 4. Error of the convolution algorithm for $N_{r(x,y)} = 64$, $N_{l(x,y)} = 8$ as compared to a dense fft based algorithm with non-existing wave-numbers set to zero. Here top left and right plots show $g_{j\ell}$ and $h_{j\ell}$ respectively, while the bottom left plot shows the resulting convolution. Note that the recurrent parts of the grid are visible as blue stripes on this plot, because these regions have less number of interactions compared to the other regions. For random initial values as presented here, multiplying the value of $f_{j\ell}$ in the horizontal region by $\sqrt{N_{rx}}/4$, and the vertical one by $\sqrt{N_{ry}}/4$ (so that the intersection is multiplied by $\sqrt{N_{rx} \times N_{ry}}/16$) would get rid of those stripes. However as can be seen on the bottom right plot, the difference between the $f_{j\ell}$ as computed by our algorithm vs. a masked dense fft algorithm is about two orders of magnitudes of machine (double) precision.</figcaption></figure>

Assuming that we can combine multiple ffts efficiently, the scaling of the algorithm would be $N_d \log N_d$. However, it is a bit worse here, since in order to have a cleaner comparison, we chose to perform the ffts in series which results. For 2D, this seems to result in an approximate empirical scaling of the form (see Fig. 5):

$$t(N_l) \propto (N_l + 4) \times N_l \log N_l \qquad (14)$$

where $N_1 = (36 \times (N_l + 3) + \frac{3}{2} \times N_r)$ from (13). Note that in comparison, for a dense fft the scaling would be $N_d^{\text{dense}} \log N_d^{\text{dense}}$ with the number of degrees of freedom for the dense fft scaling as $N_d^{\text{dense}} \approx (F_{N_r} N_r)^d \approx (N_r \phi^{N_r})^d$. In contrast, the computation time for the explicit convolution sum would scale as $(N_r + N_l)^{2d}$. Note that comparing the two terms in Eq. (13), we find that if $N_l \ll N_r/24 - 3$, the cost of additional ffts for the recurrent part of the lattice basically becomes negligible, and the cost becomes approximately independent of $N_l$. This means that for large enough $N_r$ and $N_l$, our algorithm is the only feasible way to compute a convolution.

### 3.4. How PRRL reduces real flows

In a real life application, one would have a flow, or a vorticity field that is to be described using a PRRL. In order to see how the PRRL reduction works on a realistic vorticity field, we take the results of direct numerical simulations, of a Hasegawa-Wakatani system [23], which is basically an instability driven beta plane turbulence model [24] known to make a transition into a zonal flow dominated regime [25]. Here we show a zonal flow dominated state, because it has clear large scale structure, which would be difficult to capture using a standard log-lattice reduction.

We thus take the result of direct numerical simulation and project it onto a PRRL in order to see what is lost. In terms of real space details, a PRRL is only as good as its regular part. In other words, we could basically only keep the regular part, and set everything else to zero (i.e. a low pass filter), and the spatial structure would not change visibly. Fig. 6a-c shows, if we take a 2730 × 2730 simulation result, and project it onto a (64 + 7) × (64 + 7) PRRL, we basically have a 64 × 64 resolution version of the original vorticity field in real space. However in terms of the k-spectrum, the PRRL is capable of describing the form of wave-number spectrum (i.e. it's power law form etc.) all the way up to smallest scales (see Fig. 6d). It can not however describe how the energy in these small scales are spatially localized, and thus, one would not expect it to be able to describe intermittency, at least beyond what can already be described by its low resolution regular grid part. Note that our convolution algorithm, does not introduce any additional qualitative reduction on top of this. However currently, since we do not yet have a working turbulence model using the method discussed in this paper, we are unsure about the effect of a sudden reduction of the number of triads at the boundary between the regular and the recurrent parts on the dynamics of the cascade. One may imagine that it would lead to a bottleneck effect. We observed for example in Fig. 4 that using random $g_k$ and $h_k$ fields, the recurrent part of the convolution (hence the nonlinear term) were substantially lower, because there are much less triads in these regions. We noted that this can be remedied by multiplying the value of $f_k$ in the recurrent regions by $\sqrt{N_r}/4$ in the recurrent direction. This suggests that it may be possible, to renormalize each triad that is kept with a weight that is proportional to the number of triads that are dropped, and so each of the modes that are kept in the recurrent part, can be interpreted
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Running times (in secs) for the dense fft based convolution (in orange) vs. our algorithm (in blue) for a PRRL of size $N_r[x,y] = 256$ and $N_l[x,y] = N_l$ as a function of $N_l$. This suggests that the dense fft based algorithm scales roughly as $(\phi^{N_r})^2$ (shown as the dashed orange line), while our algorithm seems to scale as $(N_l + 4) \times N_l \log N_l$ (shown as the black dashed line) where $N_l = (36 \times (N_l + 3) + \frac{3}{2} \times N_r)$ is the degrees of freedom in one dimension with the empirical numerical factor 4 relevant mainly for small (i.e. less than 512) regular grids. On the other hand, if the regular grid is large, the algorithm becomes independent of $N_l$ (see Eq. 13 and the discussion that follows). Also, the $t \propto N_l$ scaling is shown as a dashed red line. Note that the largest wave-number for the $(256 + 8) \times (256 + 8)$ case is 6964, which necessitates a padded resolution of $20898 \times 20898$. Beyond this resolution, the dense fft method becomes impractical. Note that the largest wave-number for the $(256 + 8) \times (256 + 8)$ case is 6964, which necessitates a padded resolution of $20898 \times 20898$. Beyond this resolution, the dense fft method becomes impractical.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** The effect of PRRL reduction on a realistic vorticity field dominated by zonal flows in real and Fourier space. From left to right, a) the DNS with a padded resolution of 4096 (which gives an actual resolution of $2730^2$), b) the same data projected onto a PRRL with $(64 + 7)^2$, c) the same data with a low pass filter that keeps only the first $64^2$ Fourier modes and finally d) the spectral energy density $E(k)$, computed using regular bins in logarithmic k-space. This shows that the real space resolution of the PRRL is only as good as its regular part. In contrast, for the k-spectrum, all three agree for large scales, while the low pass filter drops to zero around its maximum wave-number (i.e. $k_x = k_y = k_{64}$), the PRRL is able to describe the full spectrum up to the smallest scales.</figcaption></figure>

as representing a lumping together of the modes rather than decimation. However, this is merely a numerical observation, and it should be taken as a suggestion for the direction of future studies, and not as an established renormalization procedure.

## 4. Conclusion

The concept of PRRL is introduced as a lattice that is constructed by extending a regular grid in the wave-number domain through a recurrence relation of the form $k_n = k_{n-1} + k_{n-2}$, parametrized by two integers $N_r$, the number of regular elements and $N_l$, the number of asymptotically logarithmic, recurrent elements in each direction.

Noting that a convolution can be partitioned [see Eqs. (10) and (12)], we propose using an overlapping partitioning scheme where $n = -1$ corresponds to the regular grid, $n = 0$ to the $k = 0$ component, and $n > 1$ corresponds to the recurrent wave-numbers with $k_n = k_{n-1} + k_{n-2}$ starting with $k_{N_r/2-2}$ and $k_{N_r/2-1}$ for $n = 1$ and $2$ respectively, including negative $n$'s to correspond to negative wave-numbers unless Hermitian symmetry is desired in that direction, which is usually the last direction as in a standard complex to real fft
formulation. This partitioning scheme allows us to compute the contributions to each partition from the others. While for $n = -1$, we compute the contributions from the regular grid by considering a standard fft based convolution of the regular grid, for $n = 0$, $n = 1$ and $n = 2$ and $n > 2$ should be considered separately. For $n > 2$, we know that it is only $n \pm 1$, $n \pm 2$ and $n = 0$ that will contribute, and therefore we can consider a reduced problem with only those partitions, and compute the contribution to partition $n$ from such a computation. Furthermore since each contribution involving a recurrent wave-number can be mapped to a minimal recurrent series of only 11 elements, the computation of the contribution of any recurrent series can be reduced to an fft based convolution of only 36 (padded) elements. Repeating this procedure for each direction, we can compute the convolution in any number of dimensions. It is suggested that as the PRRL is initialized, we can already generate 4 different lists of indices corresponding to input, mapping, valid and output indices (see Fig. 3), which can be used to compute the contribution to the convolution of each different partition using the Algorithm 1.

We have demonstrated that the proposed algorithm works as intended in reproducing the convolution as computed by a dense fft based algorithm with non-existing wave-numbers explicitly set to zero, and is orders of magnitude faster than both the fft based algorithm and an explicit computation of the convolution for medium values of $N_r$ and $N_l$. We showed also that while the theoretical estimate for the algorithmic complexity goes as $N_d \log N_d$, with $N_d$ given by Eq. (13), the empirical runtime scaling is slightly worse, following the form given in Eq. 14. We also note that for $N_l \ll N_r/24 - 3$, the runtime basically becomes independent of $N_l$. Note that since $N_l$ can't be very large either (i.e., $N_l = 30$ would result in the largest wave-number being larger than $10^6$), basically having a regular grid resolution of 4096 more or less puts us in this regime.

Finally an interesting observation was that when we initialize the two convolved fields $g_k$ and $h_k$ as random variables with random phases, it was observed that the recurrent parts of the convolution (i.e. $f_k$) had substantially lower values, primarily because the points in these regions have much less triads compared to the points in the regular part of the grid. We observed that multiplying the value of $f_k$ in these regions by $\sqrt{N_r}/4$, would renormalize them, in the sense that the resulting noisy $f_k$ would be at the same level with the regular grid region. Obviously this factor depends on the form of the convolved fields $g_k$ and $h_k$, but incorporating such a weight factor may allow us to consider the Fourier space decimation implied in a PRRL, rather as a sampling of the underlying (presumably regular) grid.

In short, we developed this really efficient method for computing convolutions on asymptotically logarithmic lattices, which also have the advantage of including interactions with the $k = 0$ axes, and the development of reduced models would benefit greatly from this algorithm. In addition to developing cascade models in 2D and 3D using this algorithm, the future directions include generalizing the method described to other recurrence relations, series and combinations of series described in Ref. 16, as well as implementation of similar algorithms for developing PRRLs that asymptotically become nested polyhedra instead of log-lattices.

### CRediT authorship contribution statement

Ö.D. Gürcan: Conceptualization, Investigation, Methodology, Writing – original draft, Writing – review & editing, Software, Visualization; L. Manfredini: Conceptualization, Investigation.

### Data availability

Data will be made available on request.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgments

The authors would like to thank Dr. Benjamin Favier of IRPHE-Aix Marseille University for bringing this issue to our attention and the pointing out the interest in the numerical community for such a study. The authors would also like to thank the Isaac Newton Institute for Mathematical Sciences, Cambridge, for support and hospitality during the program “Anti-diffusive dynamics: from subcellular to astrophysical scales” during which these discussions took place. This work has benefited from a grant managed by the Agence Nationale de la Recherche (ANR), as part of the program “Investissements d’Avenir” under the reference (ANR-18-EURE-0014) and has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (BE) (Grant Agreement No. 101052200—EUROfusion) and within the framework of the French Research Federation for Fusion Studies.

### References

[1] Nussbaumer HJ. Fast Fourier transform and convolution algorithms. Berlin, Heidelberg; Springer-Verlag; 1982. 2nd ed.

[2] Press W, Teukolsky S, Vetterling W, Flannery B. Numerical recipes: the art of scientific computing; Cambridge University Press; 2007. 3rd ed.

[3] Wiener N. Extrapolation, interpolation, and smoothing of stationary time series: with engineering applications. The MIT Press; 1949. https://direct.mit.edu/book-pdf/2313079/book_9780262257190.pdf.

[4] Lecun Y, Bengio Y, Hinton G. Deep learning. Nature 2015;521:436.
[5] Frisch U. Turbulence: the legacy of A. N. Kolmogorov. Cambridge, UK: Cambridge University Press; 1995.

[6] Orszag SA. Comparison of Pseudospectral and Spectral Approximation. Stud Appl Math 1972;51:253.

[7] Canuto C, Hussaini MY, Quarteroni A, Zang TA. Spectral methods in fluid dynamics. Berlin, Heidelberg; Springer-Verlag; 1988.

[8] Gurcan OD, Li Y, Morel P. Turbulence as a Network of Fourier Modes. Mathematics 2020;8:530.

[9] Gurcan OD. Wave-number space networks in plasma turbulence. Rev Mod Plasma Phys 2023;7:20.

[10] Biferale L. Shell models of energy cascade in turbulence. Ann Rev Fluid Mech 2003;35:441.

[11] Alexakis A, Biferale L. Cascades and transitions in turbulent flows. Phys Rep 2018;767.

[12] Campolina CS, Mailybaev AA. Fluid dynamics on logarithmic lattices. Nonlinearity 2021;34:4684.

[13] Barral A, Dubrulle B. J Fluid Mech 2023;962 A2.

[14] Costa G, Barral A, Dubrulle B. Reversible Navier-Stokes equation on logarithmic lattices. Phys Rev E 2023;107:65106.

[15] Campolina CS, Mailybaev AA. Logarithmic lattice models for flows with boundaries. Physica D Nonlinear Phenomena 2025;472:134473.

[16] Manfredini L, Gurcan OD. Shell models on recurrent sequences: Fibonacci, Padovan, and other series. Phys Rev E 2025;111:25103.

[17] Aurell E, Boffetta G, Crisanti A, Frick P, Paladin G, Vulpiani A. Statistical mechanics of shell models for two-dimensional turbulence. Phys Rev E 1994;50:4705.

[18] Gurcan OD. Nested polyhedra model of turbulence. Phys Rev E 2017;95:63102.

[19] Gurcan OD. Nested polyhedra model of isotropic magnetohydrodynamic turbulence. Phys Rev E 2018;97:63111.

[20] Frisch U, Pomyalov A, Procaccia I, Ray SS. Turbulence in Noninteger Dimensions by Fractal Fourier Decimation. Phys Rev Lett 2012;108:74501.

[21] Gurcan OD, Prrlconv. 2025. https://github.com/gurcani/PRRLConv.

[22] L'vov VS, Podivilov E, Pomyalov A, Procaccia I, Vandembroucq D. Improved shell model of turbulence. Phys Rev E 1998;58:1811.

[23] Hasegawa A, Wakatani M. Plasma edge turbulence. Phys Rev Lett 1983;50:682.

[24] Gurcan OD. Internally driven β-plane plasma turbulence using the Hasegawa-Wakatani system. 2024. Physics.plasm-ph.

[25] Guillon PL, Gurcan OD. Phase transition from turbulence to zonal flows in the Hasegawa-Wakatani system. Phys Plasmas 2025;32:12306.