Rached Abdelkhalek
email: rached.abdelkhalek@total.com

Henri Calandra

Olivier Coulaud

Guillaume Latu

Jean Roman

FDTD Based Seismic Modeling and Reverse Time Migration on a GPU Cluster

Keywords: Seismic modeling, Reverse Time Migration, Finite Difference, CUDA

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Background 1.Seismic Modeling

Numerical seismic modeling aims at simulating seismic wave propagation in a geological medium in order to generate synthetic seismograms that are the seismograms that a set of sensors would record, given an assumed structure of the subsurface. Among the numerous approaches to seismic modeling, direct methods based on approximating the geological model by a numerical mesh are of particular interest. Infact, this approch can give very accurate results. However, a disadvantage of this approach is its high computational demand [START_REF] Claerbout | Imaging the Earth's Interior[END_REF].

Reverse Time Migration

Reverse Time Migration (RTM) [START_REF] Biondo | 3D Seismic Imaging[END_REF] is a technique for creating seismic images in areas of complex wave propagation, providing imaging of so called turning and prismatic waves. Its main limitation is its high computational cost. RTM is based on the simulation of waves propagation. Both source and receivers wave fields are propagated respectively forward and backward in time. These wave fields are then crosscorrelated for corresponding time steps in order to form the subsoil image.

Governing Equations

The three dimensional acoustic wave equation (1) links the pressure field u(x,y,z,t) to the density ρ(x, y, z) and the velocity c(x,y,z) [START_REF] Etgen | Computational methods for large-scale 3d acoustic finitedifference modeling: A tutorial[END_REF].

1 c 2 ρ ∂ 2 u ∂t 2 = ∇.( 1 ρ ∇u) (1)
Using finite difference (FD) methods to solve the wave equation is one way among others to tackle direct methods. The way this equation is derived among a regularly meshed domain is described in [START_REF]NVIDIA CUDA Compute Unified Device Architecture -Programming Guide[END_REF]: we write a cascaded first order spatial difference expression to compute the second time difference of the wave field:

∂ 2 u ∂t 2 = K[ ∂ - ∂x ( 1 ρ ∂ + u ∂x )+ ∂ - ∂z ( 1 ρ ∂ + u ∂y )+ ∂ - ∂z ( 1 ρ ∂ + u ∂z )] (2)
where K = c 2 ρ is the bulk modulus. The ∂ -and ∂ + symbols denote the spatial difference operators that are centered halfway between grid points either forward or backward in the direction of the spatial difference.

When the density is considered to be constant in all the domain, equation ( 1) is simplified to equation (3).

1 c 2 ∂ 2 u ∂t 2 = ∆u ( 3
)
This approximation is especially done during migration process. The discretization of this equation is done with a second order scheme in time and a eighth-order central differences in space.

GPU cluster testbed

Our GPU cluster testbed is composed of 10 Xeon bi-socket quadcore nodes (2.0 GHz clock frequency, 16 GB of RAM), coupled with 5 Nvidia Tesla S1070 servers. Each node is connected via one PCIe gen2 bus to the Tesla server. The Tesla server is composed of 4 T10 GPUs. Each pair is connected via a switch to a PCIe 2.0 connection. Thus each node has access to 2 GPUs via the same bus.

T10 GPUs can be seen as a set of 30 multiprocessors. Each multiprocessor comprises 8 streaming processors running in a single instruction multiple data (SIMD) like way. These processors can execute 3 single precision floating point operations per clock cycle. Our T10 GPUs have a clock rate of 1.44 GHz, providing thus a theoretical peak performance of 30 × 8 × 3 × 1.44 = 1 TFlops per GPU and 4 TFlops per S1070 blade. To exploit the GPU computing power and make it available for non graphics programmers, NVIDIA introduced the Compute Unified Device Architecture (CUDA) [START_REF] Carcione | Seismic modeling[END_REF]. CUDA defines a novel GPU architecture and a programming model based on a Single Program Multiple Data (SPMD) paradigm: the programmer writes a portion of code (the kernel ) to be executed by several threads in parallel on different data.

Performance and Numerical Results

Figure 2 shows the speedup and times obtained using our GPU implementation of the variable density modeling, for different subdomain decompositions. Each subdomain is treated on a host node. Reported times represent the whole application cost, including host-GPU communications. Speedup decreases with increasing number of subdomains because the computing time per subdomain decreases and the communication time becomes more predominant. Yet, only computing times are reduced when using GPUs. We report in figure 2 times and speedup for the Re- involves more host-GPU communications than the modeling. This reduces the obtained speedup.

The difference between synthetic seismograms produced using the CPU and the GPU implementations for a realistic data set used in production process is proportional to the wave field amplitude but remains very low (0.1% in terms of percentage error). This validates our implementation.

Figure 1 :Figure 2 :

 12 Figure 1: Variable Density Modeling averaged times in seconds for one iteration on a 521 × 254 × 1067 test case.


Hubert Leterme

Kévin Polisano

Valérie Perrier

Karteek Alahari

From CNNs to Shift-Invariant Twin Models Based on Complex Wavelets

We propose a novel antialiasing method to increase shift invariance and prediction accuracy in convolutional neural networks. Specifically, we replace the first-layer combination "real-valued convolutions → max pooling" (RMax) by "complex-valued convolutions → modulus" (CMod), which is stable to translations. To justify our approach, we claim that CMod and RMax produce comparable outputs when the convolution kernel is band-pass and oriented (Gaborlike filter). In this context, CMod can be considered as a stable alternative to RMax. Thus, prior to antialiasing, we force the convolution kernels to adopt such a Gaborlike structure. The corresponding architecture is called mathematical twin, because it employs a well-defined mathematical operator to mimic the behavior of the original, freely-trained model. Our antialiasing approach achieves superior accuracy on ImageNet and CIFAR-10 classification tasks, compared to prior methods based on low-pass filtering. Arguably, our approach's emphasis on retaining high-frequency details contributes to a better balance between shift invariance and information preservation, resulting in improved performance. Furthermore, it has a lower computational cost and memory footprint than concurrent work, making it a promising solution for practical implementation.

Introduction

Over the past decade, some progress has been made on understanding the strengths and limitations of convolutional neural networks (CNNs) for computer vision [START_REF] Lecun | Deep learning[END_REF]35,[START_REF] Zeiler | Visualizing and understanding convolutional networks[END_REF]. The ability of CNNs to embed input images into a feature space with linearly separable decision regions is a key factor to achieve high classification accuracy. An important property to reach this linear separability is the ability to discard or minimize non-discriminative image components. In particular, feature vectors are expected to be stable with respect to translations [35]. However, subsampling opera-tions, typically found in convolution and pooling layers, are an important source of instability-a phenomenon known as aliasing [START_REF] Azulay | Why do deep convolutional networks generalize so poorly to small image transformations[END_REF]. A few approaches have attempted to address this issue.

Blurpooled CNNs. Zhang [START_REF] Zhang | Making Convolutional Networks Shift-Invariant Again[END_REF] proposed to apply a low-pass blurring filter before each subsampling operation in CNNs. Specifically, 1. max pooling layers (Max → Sub) 1 are replaced by max-blur pooling (Max → Blur → Sub); 2. convolution layers followed by ReLU (Conv → Sub → ReLU) are blurred before subsampling (Conv → ReLU → Blur → Sub). 2 The combination Blur → Sub is referred to as blur pooling. This approach follows a well-known practice in signal processing, which involves low-pass filtering a high-frequency signal before subsampling, in order to avoid artifacts in reconstruction. Their approach improved the shift invariance as well as the accuracy of CNNs trained on ImageNet and CIFAR-10 datasets. However, this was achieved with a significant loss of information.

A question then arises: is it possible to design a nondestructive antialiasing method, and if so, does it further improve accuracy? In a more recent work, Zou et al. [START_REF] Zou | Delving Deeper into Anti-Aliasing in Con-vNets[END_REF] tackled this question through an adaptive antialiasing approach, called adaptive blur pooling. In this nonlinear setting, the blurring filter varies across channels and image locations, therefore preserving high-frequency information in strategic zones. Albeit achieving higher prediction accuracy, this approach remains fundamentally based on low-pass filtering. Consequently, features that are not blurred may still be unstable to translations. Furthermore, adaptive blur pooling requires additional memory, computational resources, and trainable parameters.

Proposed Approach. In this paper, we propose an alternative antialiasing approach based on complex-valued con-volutions, extracting high-frequency features that are stable to translations. We observed improved accuracy for ImageNet and CIFAR-10 classification, compared to the two antialiasing methods based on blur pooling [START_REF] Zhang | Making Convolutional Networks Shift-Invariant Again[END_REF][START_REF] Zou | Delving Deeper into Anti-Aliasing in Con-vNets[END_REF]. Furthermore, our approach offers significant advantages in terms of computational efficiency and memory usage, and does not induce any additional training, unlike adaptive blur pooling.

Our proposed method replaces the first layers of a CNN:

Conv → Sub → Bias → ReLU → MaxPool, (1)
equivalently written as

Conv → Sub → MaxPool → Bias → ReLU, (2)
by the following combination:

CConv → Sub → Modulus → Bias → ReLU, (3)
where CConv denotes a convolution operator with a complex-valued kernel, whose real and imaginary parts approximately form a 2D Hilbert transform pair [START_REF] Havlicek | The analytic image[END_REF]. From (2) and (3), we introduce the two following operators:

RMax : Conv → Sub → MaxPool; (4)
CMod : CConv → Sub → Modulus. ( 5
)
Our method is motivated by the following theoretical claim. In a recent paper submission [START_REF] Leterme | On the Shift Invariance of Max Pooling Feature Maps in Convolutional Neural Networks[END_REF], we proved that 1. CMod is nearly invariant to translations, if the convolution kernel is band-pass and clearly oriented; 2. RMax and CMod produce comparable outputs, except for some filter frequencies regularly scattered across the Fourier domain. We then combined these two properties to establish a stability metric for RMax as a function of the convolution kernel's frequency vector. This work was essentially theoretical, with limited experiments conducted on a deterministic model solely based on the dual-tree complex wavelet packet transform (DT-CWPT). However, it lacked applications to tasks such as image classification. Building upon this theoretical study, in this paper, we consider the CMod operator as a proxy for RMax, extracting comparable, yet more stable features.

In compliance with the theory, the RMax-CMod substitution is only applied to the output channels associated with oriented band-pass filters, referred to as Gabor-like kernels. This kind of structure is known to arise spontaneously in the first layer of CNNs trained on image datasets such as Ima-geNet [38]. In this paper, we enforce this property by applying additional constraints to the original model, prior to antialiasing. Specifically, a predefined number of convolution kernels are guided to adopt Gabor-like structures, instead of letting the network learn them from scratch. For this purpose, we rely on the dual-tree complex wavelet packet transform (DT-CWPT) [2]. Throughout the paper, we refer to this constrained model as a mathematical twin, because it employs a well-defined mathematical operator to mimic the behavior of the original model. In this context, replacing RMax by CMod is straightforward, since the complexvalued filters are provided by DT-CWPT.

Other Related Work. Following the ideas developed for antialiasing, Chaman and Dokmanic [5] reached perfect shift invariance by using an adaptive, input-dependent subsampling grid, whereas previous models rely on fixed grids. This idea was harnessed by Xu et al. [37] to get shift equivariance in generative models. This approach is not intended to compete with other antialiasing methods, but rather to complement them at the subsampling stages.

Another aspect of shift invariance in CNNs is related to boundary effects. The fact that CNNs can encode the absolute position of an object in the image by exploiting boundary effects was discovered independently by Islam et al. [START_REF] Md | How Much Position Information Do Convolutional Neural Networks Encode?[END_REF], and Kayhan and Gemert [START_REF] Osman | On Translation Invariance in CNNs: Convolutional Layers Can Exploit Absolute Spatial Location[END_REF]. This phenomenon is left outside the scope of our paper.

Wavelet scattering networks (ScatterNets), by Bruna and Mallat [START_REF] Bruna | Invariant Scattering Convolution Networks[END_REF], perform cascading wavelet convolutions and nonlinear operations. They produce shift-invariant image representations that are stable to deformation and preserve high-frequency information. Further extensions include roto-translation invariant ScatterNets [27], hybrid Scatter-Nets combined with fully-trained layers [26] or dictionary learning [40], dual-tree complex wavelet ScatterNets [START_REF] Singh | Dual-Tree wavelet scattering network with parametric log transformation for object classification[END_REF], graph ScatterNets [START_REF] Zou | Graph convolutional neural networks via scattering[END_REF], learnable ScatterNets via feature map mixing [6] or parametric wavelet filters [START_REF] Gauthier | Parametric Scattering Networks[END_REF]. As deep learning architectures with well-defined mathematical properties, ScatterNets may be used as explanatory models for standard, freely-trained networks. However, there is no exact correspondence between the two types of architectures [34]. In contrast, our models are enhanced versions of existing networks, rather than ad hoc constructions. This allows to draw comparisons in terms of prediction accuracy and shift invariance.

Finally, we draw the reader's attention to the family of complex-valued CNNs (CVCNNs), which is the focus of a comprehensive survey [22]. These networks are well suited for tasks requiring the phase information to be propagated through the entire network, as done in the context of magnetic resonance imaging [7], polarimetric imaging [START_REF] Zhang | Complex-Valued Convolutional Neural Network and Its Application in Polarimetric SAR Image Classification[END_REF] or audio signals [START_REF] Trabelsi | Deep Complex Networks[END_REF]. However, for image recognition tasks, CVCNNs do not seem to perform better than standard CNNs, with equal number of trainable parameters [START_REF] Trabelsi | Deep Complex Networks[END_REF]. Conversely, our approach discards the phase information by computing the modulus. In summary, CVCNNs and our models, although both employing complex-valued convolutions, are suited for different contexts.

Proposed Approach

We first describe the general principles of our antialiasing approach based on complex convolutions. We then provide some details about the mathematical twin based on DT-CWPT, and explain how our method has been benchmarked against blur-pooling-based antialiased models.

We represent feature maps as 2D sequences, defined by straight capital letters: X ∈ S. Indexing is denoted by square brackets: for any 2D index n ∈ Z 2 , X[n] ∈ R or C. The cross-correlation between X and V ∈ S is defined by

(X ⋆ V)[n] := k∈Z 2 X[n + k] V[k].

Standard Architectures

A convolution layer with K input channels, L output channels and subsampling factor m ∈ N \ {0} is parameterized by a weight tensor V := (V lk ) l∈{1..L}, k∈{1..K} ∈ S L×K . For any multichannel input X := (X k ) k∈{1..K} ∈ S K , the corresponding output Y := (Y l ) l∈{1..L} ∈ S L is defined such that, for any output channel l ∈ {1 . . L},

Y l := K k=1 (X k ⋆ V lk ) ↓ m. (6)
For instance, in AlexNet and ResNet, K = 3 (RGB input images), L = 64, and m = 4 and 2, respectively. Next, a bias

b := (b 1 , • • • , b L ) ⊤ ∈ R L is applied to Y,
which is then transformed through nonlinear ReLU and max pooling operators. The activated outputs satisfy

A max l := MaxPool (ReLU(Y l + b l )) , (7)
where we have defined, for any Y ∈ S and any n ∈ Z 2 ,

ReLU(Y)[n] := max(0, Y[n]); (8) MaxPool(Y)[n] := max ∥k∥ ∞ ≤1 Y[2n + k]. (9)

Antialiasing Principle

We consider the first convolution layer of a CNN, as described in (6). As widely discussed in the literature [38], after training with ImageNet, a certain number of convolution kernels V lk spontaneously take the appearance of oriented waveforms with well-defined frequency and orientation (Gabor-like kernels), as illustrated in Fig. 1a. A visual representation of trained convolution kernels is provided in Appendix E. In this paper, we refer to the corresponding output channels l ∈ G ⊂ {1 . . L} as Gabor channels. The main idea is to substitute, for any l ∈ G, RMax by CMod, as explained hereafter. Following (2), expression (7) can be rewritten where Y max l is the output of an RMax operator as introduced in [START_REF] Bruna | Invariant Scattering Convolution Networks[END_REF]. More formally,

A max l = ReLU Y max l + b l , (10)
Y max l := MaxPool K k=1 (X k ⋆ V lk ) ↓ m . (11)
Then, following (3), the RMax-CMod substitution yields

A mod l = ReLU Y mod l + b l , (12)
where Y mod l is the output of a CMod operator (5), satisfying

Y mod l := K k=1 (X k ⋆ W lk ) ↓ (2m) . (13)
In the above expression, W lk is a complex-valued analytic kernel satisfying

W lk := V lk + iH(V lk ), (14)
where H denotes the two-dimensional Hilbert transform as introduced by Havlicek et al. [START_REF] Havlicek | The analytic image[END_REF]. It is defined, for any real-valued filter V ∈ S and any frequency vector

ξ := (ξ 1 , ξ 2 ) ⊤ ∈ [-π, π] 2 , by HV(ξ) := -i sgn(ξ 1 ) V(ξ). (15)
The Hilbert transform is designed such that the Fourier transform of W lk is entirely supported in the half-plane of nonnegative x-values. Therefore, since V lk is Gaborlike, the energy of W lk is concentrated within a small window in the Fourier domain, as depicted in Fig. 1d. Due to this property, the modulus operator provides a smooth envelope for complex-valued cross-correlations with W lk [19]. This leads to the output Y mod l (13) being nearly invariant to translations. Additionally, the subsampling factor in ( 13) is twice that in [START_REF] Havlicek | The analytic image[END_REF], to account for the factor-2 subsampling achieved through max pooling (9).

Wavelet-Based Twin Models (WCNNs)

As explained in Sec. 2.2, introducing an imaginary part to the Gabor-like convolution kernels improves shift invariance. Our antialiasing method therefore restricts to the Gabor channels l ∈ G ⊂ {1 . . L}. However, G is unknown a priori: for a given output channel l ∈ {1 . . L}, whether V lk will become band-pass and oriented after training is unpredictable. Thus, we need a way to automatically separate the set G of Gabor channels from the set of remaining channels, denoted by F := {1 . . L} \ G. To this end, we built "mathematical twins" of standard CNNs, based on the dualtree wavelet packet transform (DT-CWPT). These models, which we call WCNNs in short, reproduce the behavior of freely-trained architectures with a higher degree of control and fewer trainable parameters. In this section, we present their general structure; a more detailed description is provided in Appendix A. For the purpose of readability, we assume that K = 3 (RGB input images).

We denote by L gabor := card(G) and L free := card(F) the number of Gabor and remaining channels, respectively. They are determined empirically from the trained CNNssee Tab. 1. In a twin WCNN architecture, the two groups of output channels are organized such that F = {1 . . L free } and G = {(L free + 1) . . L}. The first L free channels, which are outside the scope of our antialiasing approach, remain freely-trained as in the standard architecture. Regarding the L gabor Gabor channels, the convolution kernels V lk are constrained to satisfy the following requirements. First, all three RGB input channels are processed with the same filter, up to a multiplicative constant. More formally, there exists a luminance weight vector µ := (µ 1 , µ 2 , µ 3 ) ⊤ , with µ k ∈ [0, 1] and

3 k=1 µ k = 1, such that, ∀k ∈ {1 . . 3} , V lk = µ k V l , (16)
where V l := 3 k=1 V lk denotes the mean kernel. Furthermore, V l must be band-pass and oriented (Gabor-like filter). The following paragraphs explain how these two constraints are implemented in our WCNN architecture.

Monochrome Filters. Expression ( 16) is actually a property of standard CNNs: the oriented band-pass RGB kernels generally appear monochrome (see kernel visualization of freely-trained CNNs Appendix E). A numerical assessment of this property can by found in [START_REF] Leterme | On the Shift Invariance of Max Pooling Feature Maps in Convolutional Neural Networks[END_REF]. In WCNNs, this constraint is implemented with a trainable 1 × 1 convolution layer [START_REF] Lin | Network In Network[END_REF], parameterized by µ, computing the following luminance image:

X lum := 3 k=1 µ k X k . ( 17
)
Gabor-Like Kernels. To guarantee the Gabor-like property on V l , we implemented DT-CWPT, which is achieved through a series of subsampled convolutions. The number of decomposition stages J ∈ N \ {0} was chosen such that m = 2 J-1 , where, as a reminder, m denotes the subsampling factor as introduced in (6). DT-CWPT generates a set of filters W dt k ′ k ′ ∈{1..4×4 J } , which tiles the Fourier domain [-π, π] 2 into 4 × 4 J overlapping square windows. 3 Their real and imaginary parts approximately form a 2D Hilbert transform pair such as defined in [START_REF] Md | How Much Position Information Do Convolutional Neural Networks Encode?[END_REF]. The WCNN architecture is designed such that, for any Gabor channel l ∈ G, V l is the real part of one such filter:

∃k ′ ∈ 1 . . 4 × 4 J : V l = Re W dt k ′ . (18)
The output Y l introduced in (6) then becomes

Y l = X lum ⋆ V l ↓ 2 J-1 . ( 19
)
To summarize, a WCNN substitutes the freely-trained convolution (6) with a combination of ( 17) and ( 19), for any Gabor output channels l ∈ G. This combination is wrapped into a wavelet block, also referred to as WBlock in short. Technical details about its exact design are provided Appendix A. Note that the Fourier resolution of V lk increases with the subsampling factor m. This property is consistent with what is observed in freely-trained CNNs: in AlexNet, where m = 4, the Gabor-like filters are more localized in frequency (and less spatially localized) than in ResNet, where m = 2.

Figure 2 displays the kernels V ∈ S L×K , with K = 3 and L = 64, for the WCNN architectures based on AlexNet and ResNet, which we refer to as WAlexNet and WResNet, respectively. In these models, m = 4 and 2, which implies J = 3 and 2, respectively. The kernels are shown as RGB color images, after training with ImageNet, for both freelytrained and Gabor channels.

Antialiased WCNNs. Using the antialiasing principles presented in Sec. 2.2, we replace RMax (11) by CMod (13) for all Gabor channels l ∈ G. In the corresponding model, referred to as CWCNN, the wavelet block is replaced by a complex wavelet block (CWBlock), in which (19) becomes

Z l = X lum ⋆ W l ↓ 2 J , (20)
where W l is obtained by considering both real and imaginary parts of the DT-CWPT filter (18). Then, a modulus is applied to Z l , which yields Y mod l such as defined in ( 13), with W lk := µ k W l for any RGB channel k ∈ {1 . . 3}. Finally, we apply a bias and ReLU to Y mod l , following (12). A schematic representation of WAlexNet and its antialiased version, referred to as CWAlexNet, is provided in Figs. 3b and3c (top).

WCNNs with Blur Pooling

We benchmark our approach against the antialiasing methods proposed by Zhang [START_REF] Zhang | Making Convolutional Networks Shift-Invariant Again[END_REF] and Zou et al. [START_REF] Zou | Delving Deeper into Anti-Aliasing in Con-vNets[END_REF]. For each output channel l ∈ {1 . . 64}, the corresponding convolution kernel (V lk ) k∈{1..3} is displayed as an RGB image in the spatial domain (left), and its associated magnitude spectrum in the Fourier domain (right). The convolution kernels associated to the Gabor channels are displayed on the 4 and 3 last rows for WAlexNet and WResNet, respectively. For the sake of visual rendering, they have been respectively cropped to (11 × 11) and (7 × 7) to match the size of the freely-trained kernels.

To this end, we first consider a WCNN antialiased with static or adaptive blur pooling, respectively referred to as BlurWCNN and ABlurWCNN. A schematic representation of BlurWAlexNet is provided in Fig. 3b (bottom).

Then, we substitute the blurpooled Gabor channels (right branch of the diagram) with our own CMod-based approach. The corresponding models are respectively referred to as CBlurWCNN and CABlurWCNN. Again, a schematic representation of CBlurWAlexNet can be found in Fig. 3c (bottom). Note that, for a fair comparison, all models use blur pooling in the freely-trained channels as well as deeper layers. Therefore, several antialiasing methods are employed in different parts of the network.

Adaptation to ResNet: Batch Normalization

In many architectures including ResNet, the bias is computed after an operation called batch normalization (BN) [START_REF] Ioffe | Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift[END_REF]. In this context, (1) becomes

Conv → Sub → BN → Bias → ReLU → MaxPool. (21)
As detailed in Appendix B, the RMax-CMod substitution yields, analogously to (3),

CConv → Sub → Modulus → BN0 → Bias → ReLU, (22)
where BN0 refers to a special type of batch normalization without mean centering. A schematic representation of ResNet-based models, as done in Fig. 3 for AlexNet, is provided in Fig. 8. The bias and ReLU are placed after the modulus, following (3). In the bottom models, we compare Zhang's antialiasing approach (Fig. 3b) with ours (Fig. 3c) in the Gabor channels.

WAlexNet WResNet m (subsampling factor) 4 2 J (decomposition depth) 3 2 L free , L gabor (output channels) 32, 32 40, 24

Table 1. Experimental settings for our WCNN twin models. Other details are provided in Appendix C.

Experiments

Experiment Details

ImageNet. We built our WCNN and CWCNN twin models based on AlexNet [20] and ResNet-34 [12]. Their overall design is described in Sec. 2, along with setting details in Tab. As mentioned above, we compare blur-pooling-based antialiasing approach (Fig. 3b, bottom) with ours (Fig. 3c, bottom). To apply static or adaptive blur pooling to the WC-NNs, we proceed as follows. Following Zhang's implementation, the wavelet block is not antialiased if m = 2 as in ResNet, for computational reasons. However, when m = 4 as in AlexNet, a blur pooling layer is placed after ReLU, and the wavelet block's subsampling factor is divided by 2. Moreover, max pooling is replaced by max-blur pooling. The size of the blurring filters is set to 3, as recommended by Zhang [START_REF] Zhang | Making Convolutional Networks Shift-Invariant Again[END_REF]. Besides, DT-CWPT decompositions are performed with Q-shift orthogonal filters of length 10 as introduced by Kingsbury [18].

Our models were trained on the ImageNet ILSVRC2012 dataset [START_REF] Russakovsky | ImageNet Large Scale Visual Recognition Challenge[END_REF], following the standard procedure provided by PyTorch [28]. 4 Moreover, we set aside 100K images from the training set-100 per class-in order to compute the top-1 error rate after each training epoch ("validation set"). CIFAR-10. We also trained ResNet-18, ResNet-34 and their variants on the CIFAR-10 dataset. Training was performed on 300 epochs, with an initial learning rate equal to 0.1, decreased by a factor of 10 every 100 epochs. As for ImageNet, we set aside 5 000 images out of 50K to compute accuracy during the training phase. Given the images of small size in this dataset (32 × 32 pixels), feature extraction can be performed efficiently with a reduced number of layers. For this reason, the first layers (1) arguably have a higher influence on the overall predictive power. We therefore expect to clearly highlight the benefits of our approach on this specific task.

Evaluation Metrics

Classification Accuracy. Classification accuracy was computed on the standard ImageNet evaluation set (50K images). We followed the ten-crops procedure [20]: predictions are made over 10 patches extracted from each input image, and the softmax outputs are averaged to get the overall prediction. We also considered center crops of size 224 for one-crop evaluation. In both cases, we used top-1-5 error rates. For CIFAR-10 evaluation (10K images), we measured the top-1 error rate with one-and ten-crops.

Measuring Shift Invariance. For each image in the Im-ageNet evaluation set, we extracted several patches of size 224, each of which being shifted by 0.5 pixel along a given axis. We then compared their outputs in order to measure the model's robustness to shifts. This was done by computing the Kullback-Leibler (KL) divergence between output vectors-which, under certain hypotheses, can be interpreted as probability distributions [3, pp. 205-206]. This metric is intended for visual representation.

In addition, we measured the mean flip rate (mFR) between predictions [START_REF] Hendrycks | Benchmarking Neural Network Robustness to Common Corruptions and Perturbations[END_REF], as done by Zhang [START_REF] Zhang | Making Convolutional Networks Shift-Invariant Again[END_REF] in its blurpooled models. For each direction (vertical, horizontal and diagonal), we measured the mean frequency upon which two shifted input images yield different top-1 predictions, for shift distances varying from 1 to 8 pixels. We then normalized the results with respect to AlexNet's mFR, and averaged over the three directions. This metric is also referred to as consistency.

We repeated the procedure for the models trained on CIFAR-10. This time, we extracted patches of size 32 × 32 from the evaluation set, and computed mFR for shifts varying from 1 to 4 pixels. Besides, normalization was performed with respect to ResNet-18's mFR.

Results and Discussion

Validation and Test Accuracy.

Top-1 accuracy of AlexNet-based models along training with ImageNet is plotted in Fig. 4. Similar plots are provided for ResNet in Fig. 10. In addition, error rates of AlexNet-and ResNetbased architectures, computed on the evaluation sets, are provided in Tab. 2 for ImageNet and Tab. 3 for CIFAR-10.

Our CMod-based approach significantly outperforms the baselines for AlexNet: CWCNN vs WCNN (blue diamonds), and CBlurWCNN vs BlurWCNN (red stars). Remarkably, the exclusive application of CMod-based antialiasing to the Gabor channels (CWCNN, solid blue line) is sufficient to match the performance of blur-pooling-based antialiasing (BlurWCNN, dashed red line), which, in contrast, is implemented throughout the entire network. Positive results are also obtained for ResNet-based models Table 3. Evaluation metrics on CIFAR-10 (%): top-1 error rate using one-and ten-crops methods ("1crp" and "10crp"); and mFR measuring consistency ("shft"). Models: † static and ‡ adaptive blur pooling; * CMod-based antialiasing (our approach).

trained on ImageNet (Tab. 2, bottom). However, adaptive blur pooling, when applied to the Gabor channels (ABlur-WCNN), yields similar or marginally higher accuracy than our approach (CABlurWCNN). Nevertheless, our method is computationally more efficient, requires less memory (see "Computational Resources" below for more details), and do not demand additional training, unlike adaptive blur pooling. To better support this claim, we conducted ablation studies which we present in Appendix D.3. Finally, when trained on CIFAR-10 (see Tab. 3), our CMod-based antialiased models built upon ResNet-18 and 34 achieve significant gains in accuracy over non-antialiased models, as well as models antialiased with both blur-pooling-based methods.

As a side note, the training curves for WCNNs (colored dashed lines) closely follow those of standard CNNs (gray dotted lines). This is an expected result, since the former models are designed to mimic the behavior of the latter.

Shift Invariance (KL Divergence). The mean KL divergence between outputs of shifted images are plotted in Fig. 5 for AlexNet trained on ImageNet. Moreover, the mean flip rate for shifted inputs (consistency) is reported in Tab. 2 for ImageNet (AlexNet and ResNet-34) and Tab. 3 for CIFAR-10 (ResNet-18 and 34).

In both AlexNet and WAlexNet, the initial convolution layer's output undergoes a one-pixel shift for every fourpixel shift in the input image. Consequently, any divergence between the output vectors is due to the instability of subsequent layers to one-pixel shifts. In contrast, instabilities which are accountable to the initial layer are observed for shifts that are not multiples of 4. Likewise, for input shifts of eight pixels, the max pooling's output is shifted by exactly one pixel, resulting in even higher stability. In CWAlexNet, the same eight-to-one-pixel ratio occurs to the modulus layer's output, which explains why the two curves meet for 8-pixel shifts. However, the RMax-CMod substitution has greatly reduced first-layer instabilities, resulting in a flattened curve and avoiding the "bumps" observed for non-antialiased models.

On the other hand, BlurAlexNet and BlurWAlexNet exhibit considerably flattened curves compared to nonantialiased models, as well as CWAlexNet. This demonstrates the effectiveness of blur-pooling-based method in enhancing shift invariance. Applying CMod-based antialiasing instead of blur pooling on the Gabor channels (CBlurWAlexNet) actually degrades shift invariance (except for shifts smaller than 1.5 pixels), as evidenced by the bell-shaped curve. Nevertheless, the corresponding classifier is significantly more accurate. This is not surprising, as our approach prioritizes the conservation of high-frequency details, which are important for classification. An extreme reduction of shift variance using a large blur pooling filter would indeed result in a significant loss of accuracy. Therefore, our work achieves a better balance between shift invariance and information preservation. To gain further insights into this phenomenon, we conducted experiments by varying the size of the blurring filters in AlexNet-based models. Figure 6 shows the relationship between consistency and prediction accuracy on ImageNet (custom validation set), for different filter sizes ranging from 1 (no blur pooling) to 7 (heavy loss of high-frequency information). 5We find that a near-optimal trade-off is achieved when the filter size is set to 2 or 3. Furthermore, at equivalent consistency levels, CBlurWCNN outperforms BlurWCNN in terms of accuracy. Note however that the optimal version of CBlurWCNN is not more consistent than the optimal version of BlurWCNN, despite achieving higher accuracy.

Computational Resources. Table 4 compares the computational resources and memory footprint required for each Figure 6. Classification accuracy (ten-crops) vs consistency, measuring the stability of predictions to small input shifts, for AlexNet-based models (the lower the better for both axes). For each of the three architectures, we increased the blurring filter size from 1 (i.e., no blur pooling) to 7. The blue diamonds (no blur pooling) and red stars (blur pooling with filters of size 3) correspond to the models from Figs. 4 and5. At equivalent consistency levels, our CMod-based approach (solid line) yields higher accuracy. antialiasing method. To achieve this, we considered models with freely-trained convolutions, because the goal was to evaluate the computational performances of the various antialiasing approaches, excluding implementation tricks based on DT-CWPT from the scope of analysis. More details are provided in Appendices F and G.

Method

Conclusion

The mathematical twins introduced in this paper serve a proof of concept for our CMod-based antialiasing approach. However, its range of application extends well beyond DT-CWPT filters. While we focused on the first convolution layer, it is important to note that such initial layers play a critical role in CNNs by extracting low-level geometric features such as edges, corners or textures. Therefore, a specific attention is required for their design. In contrast, deeper layers are more focused on capturing highlevel structures that conventional image processing tools are poorly suited for, as pointed out by Oyallon et al. in the context of hybrid scattering networks [26]. Furthermore, while our approach is tailored for CNN architectures, which were chosen to make a fair comparison to related methods, it has potential for broader applicability. There is a growing interest in using self-attention mechanisms in computer vision [START_REF] Dosovitskiy | An Image is Worth 16x16 Words: Transformers for Image Recognition at Scale[END_REF] to capture complex, long-range dependencies among image representations. Recent work on vision transformers has proposed using the first layers of a CNN as a "convolutional token embedding" [10, [START_REF] Wu | CvT: Introducing Convolutions to Vision Transformers[END_REF][START_REF] Yuan | Incorporating Convolution Designs Into Visual Transformers[END_REF], effectively reintroducing inductive biases to the architecture, such as locality and weight sharing. By applying our antialiasing method to this embedding, we can provide self-attention modules with nearly shift-invariant inputs, which could be highly beneficial in improving the performance of vision transformers, especially when the amount of available data is limited.

A. Design of WCNNs: Technical Details

In this section, we provide complements to the description of the mathematical twin (WCNN) introduced in Sec. 2.3. We explained that, for each Gabor channel l ∈ G, the average kernel V l is the real part of a DT-CWPT filter, as described in (18). We now explain how the filter selection is done; in other words, how k ′ is chosen among 1 . . 4 × 4 J . Since input images are real-valued, we restrict to the filters with bandwidth located in the half-plane of positive x-values. For the sake of concision, we denote by K dt := 2 × 4 J the number of such filters.

For any RGB image X ∈ S 3 , a luminance image X lum ∈ S is computed following [START_REF] Osman | On Translation Invariance in CNNs: Convolutional Layers Can Exploit Absolute Spatial Location[END_REF], using a 1 × 1 convolution layer. Then, DT-CWPT is performed on X lum . We denote by D := (D k ) k∈{1..K dt } the tensor containing the real part of the DT-CWPT feature maps:

D k = X lum ⋆ Re W (J) k ↓ 2 J-1 . (23)
For the sake of computational efficiency, DT-CWPT is performed with a succession of subsampled separable convolutions and linear combinations of real-valued wavelet packet feature maps [START_REF] Selesnick | The dual-tree complex wavelet transform[END_REF]. To match the subsampling factor m := 2 J-1 of the standard model, the last decomposition stage is performed without subsampling.

Filter Selection. The number of dual-tree feature maps K dt may be greater than the number of Gabor channels L gabor . In that case, we therefore want to select filters that contribute the most to the network's predictive power. First, the low-frequency feature maps D 0 and D (4 J +1) are discarded. Then, a subset of K ′ dt < K dt feature maps is manually selected and permuted in order to form clusters in the Fourier domain. Considering a (truncated) permutation matrix Σ ∈ R K ′ dt ×K dt , the output of this transformation, denoted by D ′ ∈ S K ′ dt , is defined by:

D ′ := Σ D. (24)
The feature maps D ′ are then sliced into Q groups of channels D (q) ∈ S Kq , each of them corresponding to a cluster of band-pass dual-tree filters with neighboring frequencies and orientations. On the other hand, the output of the wavelet block, Y gabor := (Y l ) l∈{L free +1..L} ∈ S L gabor , where Y l has been introduced in (6), is also sliced into Q groups of channels Y (q) ∈ S Lq . Then, for each group q ∈ {1 . . Q}, an affine mapping between D (q) and Y (q) is performed. It is characterized by a trainable matrix

A (q) := α (q) 1 , • • • , α (q) Lq ⊤ ∈ R Lq×Kq such that, for any l ∈ {1 . . L q }, Y (q) l := α (q)⊤ l • D (q) . (25)
As in the color mixing stage, this operation is implemented as a 1 × 1 convolution layer.

A schematic representation of the real-and complexvalued wavelet blocks can be found in Fig. 7.

Sparse Regularization. For any group q ∈ {1 . . Q} and output channel l ∈ {1 . . L q }, we want the model to select one and only one wavelet packet feature map within the q-th group. In other words, each row vector α q) contains no more than one nonzero element, such that (25) becomes for some (unknown) value of k ∈ {1 . . K q }. To enforce this property during training, we add a mixed-norm l 1 /l ∞regularizer [25] to the loss function to penalize non-sparse feature map mixing as follows:

(q) l := α (q) l, 1 , • • • , α (q) l, Kq ⊤ of A (
Y (q) l = α (q) lk D (q) k (26)
L := L 0 + Q q=1 λ q Lq l=1 α (q) l 1 α (q) l ∞ -1 , (27)
where L 0 denotes the standard cross-entropy loss and λ ∈ R Q denotes a vector of regularization hyperparameters. Note that the unit bias in (27) serves for interpretability of the regularized loss (L = L 0 in the desired configuration) but has no impact on training.

B. Batch Normalization in ResNet

This section complements Sec. 2.5 in the main paper. In many recent architectures including ResNet, the bias (see Fig. 3) is replaced by an affine batch normalization layer (BN). In this section, we show how to adapt our approach to this context.

A BN layer is parameterized by trainable weight and bias vectors, respectively denoted by a and b ∈ R L . In the remaining of the section, we consider input images X as a stack of discrete stochastic processes. Then, expression (7) is replaced by

A l := MaxPool ReLU a l • Y l -E m [Y l ] V m [Y l ]+ε +b l , ( 28
)
with Y l satisfying (6) (output of the first convolution layer).

In the above expression, we have introduced

E m (Y l ) ∈ R and V m (Y l ) ∈ R + ,
E m [Y l ] := m 2 N 2 n∈Z 2 E[Y l [n]]; (29)
V m [Y l ] := m 2 N 2 n∈Z 2 V[Y l [n]]. (30)
In practice, estimators are computed over a minibatch of images, hence the layer's denomination. Besides, ε > 0 is a small constant added to the denominator for numerical stability. For the sake of concision, we now assume that a = 1.

Extensions to other multiplicative factors is straightforward. Let l ∈ G denote a Gabor channel. Then, recall that Y l satisfies (19) (output of the WBlock), with

V l := Re W l , (31)
where W l denotes one of the Gabor-like filters spawned by DT-CWPT. The following proposition states that, if the kernel's bandwidth is small enough, then the output of the convolution layer sums to zero.

Proposition 1. We assume that the Fourier transform of W l is supported in a region of size κ×κ which does not contain the origin (Gabor-like filter). If, moreover, κ ≤ 2π m , then

n∈Z 2 Y l [n] = 0. (32)
Proof. This proposition takes advantage of Shannon's sampling theorem. A similar reasoning can be found in the proof of Theorem 1 in [START_REF] Leterme | On the Shift Invariance of Max Pooling Feature Maps in Convolutional Neural Networks[END_REF].

In practice, the power spectrum of DT-CWPT filters cannot be exactly zero on regions with nonzero measure, since they are finitely supported. However, we can reasonably assume that it is concentrated within a region of size π/2 J-1 = π/m, as explained in [START_REF] Leterme | On the Shift Invariance of Max Pooling Feature Maps in Convolutional Neural Networks[END_REF]. Therefore, since we have discarded low-pass filters, the conditions of Prop. 1 are approximately met for W l .

We now assume that (32) is satisfied. Moreover, we assume that

E[Y l [n]] is constant for any n ∈ supp(Y l ). Aside from boundary effects, this is true if E[X lum [n]
] is constant for any n ∈ supp(X lum ). 7 We then get, for any n ∈ Z 2 , E[Y l [n]] = 0. Therefore, interchanging max pooling and ReLU yields the normalized version of (10):

A max l = ReLU Y max l E m [Y 2 l ] + ε + b l . ( 33
) (a) ResNet (b) WResNet (baseline) (c) CWResNet (proposed approach)
Figure 8. First layers of ResNet and its variants, corresponding to a convolution layer followed by ReLU and max pooling. The models are framed according to the same colors and line styles as in Fig. 10. The bias module from Fig. 3 has been replaced by an affine batch normalization layer ("BN + Bias", or "BN0 + Bias" when placed after Modulus-see Appendix B). Top: ResNet without blur pooling. Middle: Zhang's "blurpooled" models [START_REF] Zhang | Making Convolutional Networks Shift-Invariant Again[END_REF]. Bottom: Zou et al.'s approach, using adaptive blur pooling [START_REF] Zou | Delving Deeper into Anti-Aliasing in Con-vNets[END_REF].

As in Sec. 2.2, we replace Y max l by Y mod l for any Gabor channel l ∈ G, which yields the normalized version of (12):

A mod l := ReLU Y mod l E m [Y 2 l ] + ε + b l . (34)
Implementing (34) as a deep learning architecture is cumbersome because Y l needs to be explicitly computed and kept in memory, in addition to Y mod l . Instead, we want to express the second-order moment E m [Y 2 l ] (in the denominator) as a function of Y mod l . To this end, we state the following proposition.

Proposition 2. If we restrict the conditions of Prop. 1 to κ ≤ π/m, we have

∥Y l ∥ 2 2 = 2 Y mod l 2 2 . ( 35
)
Proof. This result, once again, takes advantage of Shannon's sampling theorem. The proof of our Proposition 3 in [START_REF] Leterme | On the Shift Invariance of Max Pooling Feature Maps in Convolutional Neural Networks[END_REF] is based on similar arguments.

As for Prop. 1, the conditions of Prop. 2 are approximately met. We therefore assume that (35) is satisfied, and (34) becomes

A mod l := ReLU   Y mod l 1 2 E 2m [Y mod l 2 ] + ε + b l   . ( 36
)
In the case of ResNet, the bias layer (Bias) is therefore preceded by a batch normalization layer without mean centering satisfying [START_REF] Wu | CvT: Introducing Convolutions to Vision Transformers[END_REF], which we call BN0. The second-order moment of Y mod l is computed on feature maps which are twice smaller than Y l in both directions-hence the index "2m" in [START_REF] Wu | CvT: Introducing Convolutions to Vision Transformers[END_REF], which is the subsampling factor for the CMod operator.

Schematic representations of RMax-and CMod-based CWResNet are provided in Fig. 8.

C. Experimental Settings

In this section, we provide further information that complements the experimental details presented in Sec. 3.1 and Tab. 1.

As explained in Sec. 2.3, the decomposition depth J is chosen such that m = 2 J-1 (subsampling factor). Since m = 4 in AlexNet and 2 in ResNet, we get J = 3 and 2, respectively. Therefore, the number of dual-tree filters K dt := 2 × 4 J is equal to 128 and 32, respectively.

We then manually selected K ′ dt < K dt filters. In particular, we removed the two low-pass filters, which are outside the scope of our theoretical study. Besides, for computational reasons, in WAlexNet we removed 32 "extremely" high-frequency filters which are clearly absent from the standard model (see Fig. 9a). Finally, in WRes-Net we removed the 14 filters whose bandwidths outreach the boundaries of the Fourier domain [-π, π] 2 (see Fig. 9b).

These filters indeed have a poorly-defined orientation, since a small fraction of their energy is located at the far end of the Fourier domain [2, see Fig. 1, "Proposed DT-CWPT"]. Therefore, they somewhat exhibit a checkerboard pattern. 8As explained in Appendix A, once the DT-CWPT feature maps have been manually selected, the output D ′ ∈ S K ′ dt is sliced into Q groups of channels D (q) ∈ S Kq . For each group q, a depthwise linear mapping from D (q) to a bunch of output channels Y (q) ∈ S Lq is performed. Finally, the wavelet block's output feature maps Y gabor ∈ S L gabor are obtained by concatenating the outputs Y (q) depthwise, for (a) WAlexNet (J = 3) (b) WResNet (J = 2) Figure 9. Mapping scheme from DT-CWPT feature maps D ∈ S K dt to the wavelet block's output Y gabor ∈ S L gabor . Each wavelet feature map is symbolized by a small square in the Fourier domain, where its energy is mainly located. The gray areas show the feature maps which have been manually removed. Elsewhere, each group of feature maps D (q) ∈ S Kq is symbolized by a dark frame-in (b), Kq is always equal to 1. For each group q ∈ {1 . . Q}, a number indicates how many output channels Lq are assigned to it. The colored numbers in (a) refer to groups on which we have applied l ∞ /l 1 -regularization. Note that, when inputs are real-valued, only the half-plane of positive x-values is considered.

any q ∈ {1 . . Q}. Figure 9 shows how the above grouping is made, and how many output channels L q each group q is assigned to.

During training, the above process aims at selecting one single DT-CWPT feature map among each group. This is achieved through mixed-norm l ∞ /l 1 regularization, as introduced in (27). The regularization hyperparameters λ q have been chosen empirically. If they are too small, then regularization will not be effective. On the contrary, if they Model Filt. frequency Reg. param.

WAlexNet

[π/8, π/4[ - [π/4, π/2[ 4.1 • 10 -3 [π/2, π[ 3.2 • 10 -4
WResNet any -Table 5. Regularization hyperparameters λq for each group q of DT-CWPT feature maps. The groups with only one feature map do not need any regularization since this feature map is automatically selected. The second and third rows of WAlexNet correspond to the blue and magenta groups in Fig. 9a, respectively.

are too large, then the regularization term will become predominant, forcing the trainable parameter vector α (q) l to randomly collapse to 0 except for one element. The chosen values of λ q are displayed in Tab. 5.

Finally, the split L free -L gabor between the freely-trained and Gabor channels, provided in the last row of Tab. 1 (main document), have been empirically determined from the standard models. More specifically, considering standard AlexNet and ResNet-34 trained on ImageNet (see Fig. 12), we determined the characteristics of each convolution kernel: frequency, orientation, and coherence index (which indicates whether an orientation is clearly defined). This was done by computing the tensor structure [START_REF] Jahne | Practical Handbook on Image Processing for Scientific and Technical Applications[END_REF]. Then, by applying proper thresholds, we isolated the Gabor-like kernels from the others, yielding the approximate values of L free and L gabor . Furthermore, this procedure allowed us to draw a rough estimate of the distribution of the Gabor-like filters in the Fourier domain, which was helpful to design the mapping scheme shown in Fig. 9.

D. Additional Figures and Tables

D.1. Validation Error Along Training

As for AlexNet in Fig. 4, Top-1 validation curves of ResNet-based models along training with ImageNet and CIFAR-10 are plotted in Fig. 10. In many cases, accuracy increases when applying our CMod-based antialiasing method (dashed versus solid lines). However, in some cases, the gain on the validation set (a subset of images put aside from the training set) is marginal, for instance BlurWResNet-34 versus CBlurWResNet-34 trained on Im-ageNet (Fig. 10a, red stars), or non-existent, for instance WResNet-18 versus CWResNet-18 trained on CIFAR-10 (Fig. 10b, blue diamonds). Yet, notable gains are observed on the evaluation set (a separate dataset provided by Ima-geNet or CIFAR), as reported in Tabs. 2 and 3. This indicates a better generalization capability of our approach, compared to the static blur pooling-based method, or no antialiasing. In particular, for ImageNet classification (Tab. 2, bottom), simply replacing the blurpooled Gabor channels in the first layer with our CMod-based approach (BlurWCNN versus CBlurWCNN) produces improvements on a similar order of magnitude as adaptive blur pooling, which is applied throughout the whole network (ABlurWCNN).

D.2. Accuracy vs Consistency

Figure 11 shows the relationship between consistency and prediction accuracy of AlexNet and ResNet-based models on ImageNet, for different filter sizes ranging from 1 (no blur pooling) to 7 (heavy loss of high-frequency information). The data for AlexNet on the validation set are displayed in the main document, Fig. 6. As recommended by Zhang [START_REF] Zhang | Making Convolutional Networks Shift-Invariant Again[END_REF], the optimal trade-off is generally achieved when the blurring filter size is equal to 3. Moreover, in either case, at equivalent level of consistency, replacing blur pooling by our CMod-based antialiasing approach in the Gabor channels increases accuracy.

D.3. Ablation Study

In Sec. 3.3, we showed that our models outperform the baselines, except when adaptive blur pooling is used as an antialiasing method. In that case, although applying our CMod-based approach on the Gabor channels instead of adaptive blur pooling leads to slightly degraded accuracy, it comes with far less computational resources and memory footprint, as shown in Tab. 4. However, one may wonder: what happens if we do not employ any antialiasing on the Gabor channels (and keep on using blur pooling in the rest of the network)? Do we reach similar accuracy as well? To answer this question, we trained such ablated models, in both static and adaptive blur-pooling situations.

We found that removing antialiasing from the Gabor channels (either blur pooling or CMod) generally leads to decreased accuracy and consistency, as shown in Tab. 6. However, when using adaptive blur pooling, the primary improvements in performance come from other parts of the network-as evidenced when comparing the first and last rows. This limits the influence of antialiasing on the Gabor channels. Yet, although adaptive blur pooling yields the The freely-trained models, upon which the mathematical twins are built, appear in faint gray. Legend: † static and ‡ adaptive blur pooling; * CMod-based antialiasing (our approach). most important gains in performance, it comes at a substantial increase in computational resources due to the need for Figure 11. Classification accuracy (ten-crops) vs consistency, measuring the stability of predictions to small input shifts (the lower the better both axes). The metrics have been computed on ImageNet-1K, on both validation set (100K images set aside from the training set) and test set (50K images provided as a separate dataset). For each model (BlurCNN, BlurWCNN and CBlurWCNN), we increased the blurring filter size from 1 (i.e., no blur pooling) to 7. The blue diamonds (no blur pooling) and red stars (blur pooling with filters of size 3) correspond to the models from Fig. 4 for AlexNet and Fig. 10a for ResNet (models trained after 90 epochs). multiple such antialiasing layers throughout the network.

E. Kernel Visualization (Standard Models)

The convolution kernels V := (V lk ) l∈{1..64}, k∈{1..3} ∈ S 64×3 satisfying [START_REF] Havlicek | The analytic image[END_REF] Table 7. Computation time per operation, determined experimentally using PyTorch (CPU computations). They have been normalized with respect to the computation time of an addition.

them with WCNN kernels in Fig. 2, we can notice that, up to a few exceptions, the freely-trained channels (4 and 5 first rows for AlexNet and ResNet, respectively) have been specialized to lower-frequency filters (mono-or bicolor blobs).

F. Computational Cost

This section provides technical details about our estimation of the computational cost (FLOPs), such as reported in Tab. 4, for one input image and one Gabor channel. This metric was estimated in the case of standard 2D convolutions.

F.1. Average Computation Time per Operation

The estimated computation time for each type of operation has been reported in Tab. 7.

F.2. Computational Cost per Layer

In the following paragraphs, L ∈ N \ {0} denotes the number of output channels (depth) and N ′ ∈ N \ {0} denotes the size of output feature maps (height and width). However, note that N ′ is not necessary the same for all layers. For instance, in standard ResNet, the output of the first convolution layer is of size N ′ = 112, whereas the output of the subsequent max pooling layer is of size N ′ = 56. For each type of layer, we calculate the number of FLOPs required to produce a single output channel l ∈ {1 . . L}. Moreover, we assume, without loss of generality, that the model processes one input image at a time.

Convolution Layers. Inputs of size (K × N × N ) (input channels, height and width); outputs of size (L × N ′ × N ′ ). For each output unit, a convolution layer with kernels of size (N filt × N filt ) requires KN 2 filt multiplications and KN 2 filt -1 additions. Therefore, the computational cost per output channel is equal to

T conv = N ′ 2 (KN 2 filt -1) • t s + KN 2 filt • t p . (37)
Complex Convolution Layers. Inputs of size (K × N × N ); complex-valued outputs of size (L × N ′ × N ′ ). For each output unit, a complex-valued convolution layer requires 2 × KN 2 filt multiplications and 2 × (KN 2 filt -1) additions. Computational cost per output channel:

T C conv = 2N ′ 2 (KN 2 filt -1) • t s + KN 2 filt • t p . ( 38
)
Note that, in our implementations, the complex-valued convolution layers are less expensive than the real-valued ones, because the output size N ′ is twice smaller, due to the larger subsampling factor.

Bias and ReLU. Inputs and outputs of size (L×N ′ ×N ′ ).

One evaluation for each output unit:

T bias = N ′ 2 t s and T relu = N ′ 2 t relu . (39)
Max Pooling. Outputs of size (L×N ′ ×N ′ ), with N ′ depending on whether subsampling is performed at this stage (no subsampling when followed by a blur pooling layer). One evaluation for each output unit:

T max = N ′ 2 t max . (40)
Modulus Pooling. Complex-valued inputs and realvalued outputs of size (L × N ′ × N ′ ). One evaluation for each output unit:

T mod = N ′ 2 t mod . (41)
Batch Normalization. Inputs and outputs of size (L × N ′ × N ′ ). A batch normalization (BN) layer, described in (28), can be split into several stages.

1. Mean: N ′ 2 additions.

2. Standard deviation: N ′ 2 multiplications, N ′ 2 additions (second moment), N ′ 2 additions (subtract squared mean).

3. Final value: N ′ 2 additions (subtract mean), 2N ′ 2 multiplications (divide by standard deviation and multiplicative coefficient).

Overall, the computational cost per image and output channel of a BN layer is equal to

T bn = N ′ 2 (4 t s + 3 t p ) . (42)
T blur = N ′ 2 (N 2 b -1) • t s + N 2 b • t p . (43)
Adaptive Blur Pooling. Inputs of size (L × 2N ′ × 2N ′ ); outputs of size (L × N ′ × N ′ ). An adaptive blur pooling layer [START_REF] Zou | Delving Deeper into Anti-Aliasing in Con-vNets[END_REF] with filters of size (N b × N b ) splits the L output channels into Q := L/L g groups of L g channels that share the same blurring filters. The adaptive blur pooling layer can be decomposed into the following stages.

1. Generation of blurring filters using a convolution layer with trainable kernels of size (N b ×N b ): inputs of size

(L × 2N ′ × 2N ′ ), outputs of size (QN 2 b × N ′ × N ′ ).
For each output unit, this stage requires LN 2 b multiplications and LN 2 b -1 additions. The computational cost divided by the number L of channels is therefore equal to

T conv ablur = N ′ 2 N 2 b L g × (LN 2 b -1) • t s + LN 2 b • t p . ( 44
)
Note that, despite being expressed on a per-channel basis, the above computational cost depends on the number L of output channels. This is due to the asymptotic complexity of this stage in O(L 2 ).

2. Batch normalization, inputs and outputs of size

(QN 2 b × N ′ × N ′ ): T bn ablur = N ′ 2 N 2 b L g (4 t s + 3 t p ) . (45)
3. Softmax along the depthwise dimension:

T sftmx ablur = N ′ 2 N 2 b L g (t e + t s + t p ). (46)
4. Blur pooling of input feature maps, using the filter generated at stages (1)-( 3): inputs of size (L×2N ′ ×2N ′ ), outputs of size (L × N ′ × N ′ ). The computational cost per output channel is identical to the static blur pooling layer, even though the weights may vary across channels and spatial locations:

T blur = N ′ 2 (N 2 b -1) • t s + N 2 b • t p . (47)
Overall, the computational cost of an adaptive blur pooling layer per input image and output channel is equal to

T ablur = N ′ 2 N 2 b L g (L + 1)N 2 b + 3 • t s + (L + 1)N 2 b + 4 • t p + t e . (48)
We notice that an adaptive blur pooling layer has an asymptotic complexity in O(N 4 b ), versus O(N 2 b ) for static blur pooling.

F.3. Application to AlexNet-and ResNet-based Models

Since they are normalized by the computational cost of standard models, the FLOPs reported in Tab. 4 only depend on the size of the convolution kernels and blur pooling filters, respectively denoted by N filt and N b ∈ N \ {0}. In addition, the computational cost of the adaptive blur pooling layer depend on the number of output channels L as well as the number of output channels per group L g .

In practice, N filt is respectively equal to 11 and 7 for AlexNet-and ResNet-based models. Moreover, N b = 3, L = 64 and L g = 8. Actually, the computational cost is largely determined by the convolution layers, including step (1) of adaptive blur pooling.

G. Memory Footprint

This section provides technical details about our estimation of the memory footprint for one input image and one output channel, such as reported in Tab. 4. This metric is generally difficult to estimate, and is very implementationdependent. Hereafter, we consider the size of the output tensors, as well as intermediate tensors saved by torch. autograd for the backward pass. However, we didn't take into account the tensors containing the trainable parameters. To get the size of intermediate tensors, we used the Python package PyTorchViz. 9 These tensors are saved according to the following rules.

• Convolution (Conv), batch normalization (BN), Bias, max pooling (MaxPool or Max), blur pooling (Blur-Pool), and Modulus: the input tensors are saved, not the output. When Bias follows Conv or BN, no intermediate tensor is saved.

• ReLU, Softmax: the output tensors are saved, not the input.

• If an intermediate tensor is saved at both the output of a layer and the input of the next layer, its memory is not duplicated. An exception is Modulus, which stores the input feature maps as complex numbers.

• MaxPool or Max: a tensor of indices is kept in memory, indicating the position of the maximum values. The tensors are stored as 64-bit integers, so they weight twice as much as conventional float-32 tensors.

• BN: four 1D tensors of length L are kept in memory: computed mean and variance, and running mean and variance. For BN0 (see Appendix B), where variance is not computed, only two tensors are kept in memory.

In the following paragraphs, we denote by L the number of output channels, N the size of input images (height and width), m the subsampling factor of the baseline models (4 for AlexNet, 2 for ResNet), N b the blurring filter size (set to 3 in practice). For each model, a table contains the size of all saved intermediate or output tensors. For example, the values associated to "Layer1 → Layer2" correspond to the depth (number of channel), height and width of the intermediate tensor between Layer1 and Layer2. Then, the memory footprint for each output channel is equal to

=⇒ S std = 7 4
N 2 m 2 .

  The down arrow refers to subsampling: for any m ∈ N * , (X ↓ m)[n] := X[mn].

Figure 1 .

 1 Figure 1. (a), (b): Real and imaginary parts of a Gabor-like convolution kernel W lk := V lk + iH(V lk ), forming a 2D Hilbert transform pair. (c), (d): Power spectra (energy of the Fourier transform) of V lk and W lk , respectively.

  Figure 2. Convolution kernels in the first layer of WAlexNet (a) and WResNet-34 (b), after training with ImageNet-1K.For each output channel l ∈ {1 . . 64}, the corresponding convolution kernel (V lk ) k∈{1..3} is displayed as an RGB image in the spatial domain (left), and its associated magnitude spectrum in the Fourier domain (right). The convolution kernels associated to the Gabor channels are displayed on the 4 and 3 last rows for WAlexNet and WResNet, respectively. For the sake of visual rendering, they have been respectively cropped to (11 × 11) and (7 × 7) to match the size of the freely-trained kernels.

Figure 3 .

 3 Figure 3. First layers of AlexNet and its variants, corresponding to a convolution layer followed by ReLU and max pooling[START_REF] Azulay | Why do deep convolutional networks generalize so poorly to small image transformations[END_REF]. The models are framed according to the same colors and line styles as in Figs.4 and 5. The green modules are the ones containing trainable parameters; the orange and purple modules represent static linear and nonlinear operators, respectively. The numbers between each module represent the depth (number of channels), height and width of each output. Fig.3a: freely-trained models. Top: standard AlexNet. Bottom: Zhang's "blurpooled" AlexNet. Fig.3b: mathematical twins (WAlexNet) reproducing the behavior of standard (top) and blurpooled (bottom) AlexNet. The left side of each diagram corresponds to the L free := 32 freely-trained output channels, whereas the right side displays the L gabor := 32 remaining channels, where freely-trained convolutions have been replaced by a wavelet block (WBlock) as described in Sec. 2.3. Fig.3c: CMod-based antialiased WAlexNet, where WBlock has been replaced by CWBlock, and max pooling by a modulus. The bias and ReLU are placed after the modulus, following (3). In the bottom models, we compare Zhang's antialiasing approach (Fig.3b) with ours (Fig.3c) in the Gabor channels.

Figure 4 .

 4 Figure 4. Evolution of the top-1 validation error (one-crop) along training with ImageNet, for AlexNet-based models, as described in Fig.3. The freely-trained models, upon which the mathematical twins are built, appear in faint gray. Legend: † blur pooling; * CMod-based antialiasing (our approach).

Figure 5 .

 5 Figure 5. AlexNet-based models: mean KL divergence between the outputs of a reference image versus shifted images. Legend: † blur pooling; * CMod-based antialiasing (our approach).

Figure 7 .

 7 Figure 7. Detail of a wavelet block with J = 3 as in AlexNet, in its RMax (left) and CMod (right) versions. DT-RWPT corresponds to the real part of DT-CWPT.

  Figure10. Evolution of the top-1 validation error (one-crop) along training, for ResNet-based models, as described in Fig.8. The freely-trained models, upon which the mathematical twins are built, appear in faint gray. Legend: † static and ‡ adaptive blur pooling; * CMod-based antialiasing (our approach).

  (a) AlexNet, test set (50K images) (b) ResNet-34, validation set (100K images) (c) ResNet-34, test set (50K images)

  are shown in Fig. 12 for freely-trained AlexNet and ResNet-34 trained on ImageNet. The kernels are shown as RGB color images. When comparing ts =

Figure 12 .

 12 Figure 12. Convolution kernels in the first layer of freely-trained AlexNet (a) and ResNet-34 (b), after training with ImageNet-1K. For each output channel l ∈ {1 . . 64}, the corresponding convolution kernel (V lk ) k∈{1..3} is displayed as an RGB image in the spatial domain (left), and its associated magnitude spectrum in the Fourier domain (right).

G. 1 .

 1 AlexNet-based Models No Antialiasing. Conv → Bias → ReLU → MaxPool. ReLU → MaxPool L

9m 2 .

 2 https://github.com/szagoruyko/pytorchvizStatic Blur Pooling. Conv → Bias → ReLU → BlurPool → Max → BlurPool. Approach. CConv → Modulus → Bias → ReLU. ResNet-based Models No Antialiasing. Conv → BN → Bias → ReLU → MaxPool. Static Blur Pooling. Conv → BN → Bias → ReLU → Max → BlurPool. Adaptive Blur Pooling. Conv → BN → Bias → ReLU → Max → ABlurPool. Approach. CConv → Modulus → BN0 → Bias → ReLU.

Table 2 .

 2 Evaluation metrics on ImageNet (%): the lower the better. Models: † static and ‡ adaptive blur pooling;

	Model	One-crop top-1 top-5 top-1 top-5 Ten-crops	Shifts mFR
				AlexNet
	CNN	45.3	22.2	41.3	19.3	100.0
	WCNN	44.9	21.8	40.8	19.0	101.4
	CWCNN *	44.3	21.3	40.2	18.5	88.0
	BlurCNN †	44.4	21.6	40.7	18.7	63.8
	BlurWCNN †	44.3	21.4	40.5	18.5	63.1
	CBlurWCNN † *	43.3	20.5	39.6	17.9	69.4
				ResNet-34
	CNN	27.6	9.2	24.8	7.7	78.1
	WCNN	27.4	9.2	24.7	7.6	77.2
	CWCNN *	27.2	9.0	24.4	7.4	73.1
	BlurCNN †	26.7	8.6	24.0	7.2	61.2
	BlurWCNN †	26.7	8.6	24.1	7.3	65.2
	CBlurWCNN † *	26.5	8.4	23.7	7.0	62.5
	ABlurCNN ‡	26.1	8.3	23.5	7.0	60.8
	ABlurWCNN ‡	26.0	8.2	23.6	6.9	62.1
	CABlurWCNN ‡ *	26.1	8.2	23.7	7.0	63.1
	Model	ResNet-18 1crp 10crp shft 1crp 10crps shft ResNet-34
	CNN 14.9	10.8 100.0 15.2	10.9 100.3
	WCNN 14.2	10.3 92.4 14.5	10.5 99.2
	CWCNN * 13.8	9.6 88.8 12.9	9.2 93.0
	BlurCNN † 14.2	10.4 87.7 15.7	11.6	88.2
	BlurWCNN † 13.1	9.7 84.6 13.2	9.9 85.6
	CBlurWCNN † * 12.3	8.9 85.7 12.4	9.1 83.7
	ABlurCNN ‡ 14.6	11.0 90.9 16.3	12.8 91.9
	ABlurWCNN ‡ 14.5	11.0 86.5 14.0	10.4 93.3
	CABlurWCNN ‡ * 12.8	9.7 81.7 12.8	9.2 86.6

* CMod-based antialiasing (our approach).

Table 4 .

 4 Computational cost and memory footprint required for each antialiasing method, per Gabor channel. The values are nor-

		Computational cost	Memory footprint
		AlexNet	ResNet AlexNet ResNet
	No antialiasing (ref)	1 .0	1 .0	1 .0	1 .0
	BlurPool [42]	4.0	1.0	4.7	1.9
	ABlurPool [45]	-	2.1	-	2.0
	CMod (ours)	0.5	0.5	0.6	0.4
	malized relative to non-antialiased AlexNet or ResNet. Computa-
	tional cost: FLOPs for computing Y max l	(11) or Y mod l	(13). Mem-
	ory footprint: size of the intermediate and output tensors saved by
	PyTorch for the backward pass.

  which respectively denote the mean expected value and variance of Y

l [n], for indices n contained in the support of Y l , denoted by supp(Y l ). Let us denote by N ∈ N\{0} the support size of input images. Therefore, if the filter's support size N filt is much smaller that N , then supp(Y l ) is roughly of size N/m. We thus define the above quantities as follows:

Table 6 .

 6 Ablation study on WResNet. Evaluation metrics on Ima-geNet (%): the lower the better. Models: † static and ‡ adaptive blur pooling; * CMod-based antialiasing on the Gabor channels (our approach); ♢ ablated model: no antialiasing on the Gabor channels (neither blur pooling nor CMod); ♦ no antialiasing anywhere.

	Model	One-crop top-1 top-5 top-1 top-5 Ten-crops	Shifts mFR
	WCNN ♦	27.4	9.2	24.7	7.6	77.2
	BlurWCNN †	26.7	8.6	24.1	7.3	65.2
	CBlurWCNN † *	26.5	8.4	23.7	7.0	62.5
	→ BlurWCNN †♢	26.8	8.6	24.3	7.1
	ABlurWCNN ‡	26.0	8.2	23.6	6.9	62.1
	CABlurWCNN ‡ *	26.1	8.2	23.7	7.0	63.1
	→ ABlurWCNN †♢	26.4	8.4	23.9	7.0	70.3

Sub and Conv stand for "convolution" and "subsampling" respectively.

2 ReLU is computed before blurring; otherwise the network would simply perform on low-resolution images.

Figure 1 actually illustrates one of these filters, with J = 3.

PyTorch "examples" repository available at https://github. com/pytorch/examples/tree/main/imagenet

Similar plots can be found in Zhang's paper[START_REF] Zhang | Making Convolutional Networks Shift-Invariant Again[END_REF].

https://gricad.univ-grenoble-alpes.fr

This property is a rough approximation for images of natural scenes or man-made objects. In practice, the main subject is generally located at the center, the sky at the top, etc. These are sources of variability for color and luminance distributions across images, as discussed by Torralba and Oliva[START_REF] Torralba | Statistics of natural image categories[END_REF].

Note that the same procedure could have been applied to WAlexNet, but it was deemed unnecessary because the boundary filters were spontaneously discarded during training.

Acknowledgments

This work has been partially supported by the LabEx PERSYVAL-Lab (ANR-11-LABX-0025-01) funded by the French program Investissement d'avenir, as well as the ANR grants MIAI (ANR-19-P3IA-0003) and AVENUE (ANR-18-CE23-0011). Most of the computations presented in this paper were performed using the GRICAD infrastructure, 6 which is supported by Grenoble research communities.