On the Fourier analysis in the SO(3) space : EquiLoPO Network
Dmitrii Zhemchuzhnikov Sergei Grudinin


Toeplitz operators and zeros of square-integrable random holomorphic sections
    Alexander Drewitz, Bingxiao Liu and George Marinescu
    April 29, 2024
=============================================================================




Analyzing volumetric data with rotational invariance or equivariance is an active topic in current research. Existing deep-learning approaches utilize either group convolutional networks limited to discrete rotations or steerable convolutional networks with constrained filter structures. 
This work proposes a novel equivariant neural network architecture that achieves analytical Equivariance to Local Pattern Orientation on the continuous SO(3) group while allowing unconstrained trainable filters - EquiLoPO Network.
Our key innovations are a group convolutional operation leveraging irreducible representations as the Fourier basis and a local activation function in the SO(3) space that provides a well-defined mapping from input to output functions, preserving equivariance. By integrating these operations into a ResNet-style architecture, we propose a model that overcomes the limitations of prior methods. A comprehensive evaluation on diverse 3D medical imaging datasets from MedMNIST3D demonstrates the effectiveness of our approach, which consistently outperforms state of the art. This work suggests the benefits of true rotational equivariance on SO(3) and flexible unconstrained filters enabled by the local activation function,  providing a flexible framework for equivariant deep learning on volumetric data with potential applications across domains. Our code is publicly available at <https://gricad-gitlab.univ-grenoble-alpes.fr/GruLab/ILPO/-/tree/main/EquiLoPO>.




§ INTRODUCTION




Deep-learning methods have shown remarkable success in analyzing spatial data across various domains. However, in many real-world scenarios, the data can be presented in arbitrary orientations. Therefore, the output of the neural network should be invariant or equivariant to rotations of the input. While data augmentation can partially address this requirement, it leads to increased computational demands, especially for volumetric data 
with three rotation angles.

To tackle this challenge, researchers have developed equivariant neural network architectures that utilize rotationally equivariant operations. We can broadly categorize these  methods into two classes: group convolutional networks and steerable convolutional networks. Group convolutional networks achieve rotational equivariance by convolving data in both translational and rotational spaces, but they are typically limited to a discrete set of rotations. On the other hand, steerable convolutional networks employ filters that are analytically equivariant to continuous rotations, but they impose constraints on the filter structures.

In this work, we propose a novel equivariant neural network architecture that combines the strengths of both approaches. Our key contributions are:

    
  * We present a group convolutional network that achieves analytical equivariance with respect to the continuous rotational space SO(3) by leveraging irreducible representations as the Fourier basis.
    
  * Contrary to steerable convolutional networks, our approach does not impose constraints on the filter structures beyond the finite resolution limits.
    
  * We present a local activation function in the rotational space that provides a well-defined mapping from input function values to output function values, ensuring the preservation of equivariance properties while avoiding the reduction of the architecture to a steerable convolutional network.



By addressing the limitations of existing methods, our approach offers a powerful and flexible framework for analyzing volumetric data in a rotationally equivariant manner without the need for data augmentation or constraints on filter structures.



§ RELATED WORK



 §.§ Deep learning for volumetric data and equivariance


In recent years, deep learning has firmly established its place in the analysis of spatial data. Neural networks learn a hierarchy of features by recognizing spatial patterns at various levels. However, in real-world scenarios, multidimensional data are often given in arbitrary orientation and shift, as the coordinate system is not defined. In such cases, it is desirable for the output of the neural network to be independent of the rotation or shift of the input data. While dealing with shifts is straightforward through the use of convolution, which is inherently shift-equivariant, achieving equivariance with respect to rotations requires additional efforts. A primary solution for achieving this effect is augmentation of the training dataset with rotated samples <cit.> . Augmentation leads to an increase not only in the dataset size and consequently training time, but also in the number of network parameters needed to memorize patterns in multiple orientations. In volumetric cases, this increased demand for computational resources can be overwhelming, as there are three angles of rotation in 3D, compared to just one in 2D. For some data types, a canonical coordinate system may be defined <cit.>. However, in most real-world scenarios, such a coordinate system cannot be identified, and even within a canonical coordinate system, the same local patterns may be encountered in different orientations. These circumstances have led the community to focus on methods that are analytically invariant or equivariant to rotations, utilizing rotationally equivariant operations. 



Rotational equivariant methods for regular data can be divided into two groups: group convolution networks and steerable convolutional networks. Our method formally belongs to the first group but without a specific approach to activation in the Fourier space can be reduced to a method from the second group as it is shown in Appendix <ref>. 
Thus,
we will briefly describe below
group convolution networks, spherical harmonic networks and activation in the Fourier space.




 §.§ Group Convolution networks

The pioneering method from the first class was the Group Equivariant Convolutional Networks (G-CNNs) introduced by <cit.>, who proposed a general view on convolutions in different group spaces. Many more methods  were built up subsequently upon this approach <cit.>. Several implementations of Group Equivariant Networks were specifically adapted for regular volumetric data, e.g.,
CubeNet<cit.> and 3D G-CNN <cit.>.  
Methods of this class achieve rotational equivariance by convolving data not only in translational but also in the rotational space. 
The authors of these methods consider a discrete set of 90-degree rotations, which exhaustively describe the possible positions of a cubic pattern on a regular grid. 
However, equivariance with respect to this discrete group of rotations does not guarantee  equivariance on the continuous group SO(3).
Separable SE(3)-equivariant network <cit.> approximates equivariance in SO(3) but lacks analytical equivariance since the authors sample only a finite set of points in the SO(3) space. 
Analytical rotational equivariance in 3D can be achieved by using irreducible representations in the O(2) or the SO(3) group. 



 §.§ Spherical harmonics networks

 <cit.> introduced steerable networks, a class of methods that use analytically equivariant filters with respect to a particular group of transformations. 
The first two approaches that applied analytically-equivariant filters to the 3D data are the Tensor Field Networks (TFN) <cit.> and the N-Body Networks (NBNs) <cit.>. 
<cit.> presented a similar approach, the Clebsch-Gordan Nets, applied to data on a sphere. 
<cit.> presented steerable networks for the regular three-dimensional  data. The authors deduced a complete equivariant kernel basis where input and output irreducible features have arbitrary degrees. 
The filter must belong to the subspace of the equivariant kernels. This requirement may limit the expressiveness of the network <cit.>.
To summarize, group convolutional networks do not put constraints on filters but provide equivariance only on a discrete space of rotations. Spherical harmonics networks use irreducible representations to achieve analytical equivariance on the continuous rotational space but imply constraints on the filters.
Our method employs irreducible SO(3) representations (Wigner matrices)  as the Fourier basis in a group convolutional network. This approach allows us to obtain analytical equivariance and avoid constraints on the filter shape. In a general setup, such a convolution can be reduced to a steerable network, shown in Appendix <ref>, because an SO(3) irrep can be seen as a set of O(2) irreps. 
However, we introduce a novel activation function in SO(3) that prevents the reduction of the whole architecture to a steerable network.
Indeed, our activation 
is not permutationally invariant with respect to Wigner coefficients of different orders and the same degree and thus cannot be reduced to   a steerable network, as we demonstrate in Appendix <ref>.





 §.§ Activation Operators on Irreducible Representations


Activation functions play a crucial role in neural network architectures as they introduce nonlinear operations. 
When working with irreducible representations (irreps) in rotationally-equivariant neural networks, it is mandatory to choose activation functions that preserve the  equivariance properties.
Several approaches, listed below, have been proposed to apply activation functions to irreps.
Norm-based activation functions operate on the norm (magnitude) of each irrep, preserving the equivariance properties. The L2 norm of each irrep is computed, and a scalar activation function is applied to the norm. The activated norm is then used to scale the original irrep <cit.>. The same principle was used for the Fourier decomposition of a function in 3D by <cit.>.

Gated activation functions introduce learnable parameters to control the activation of each irrep <cit.>. A separate set of learnable weights is used to compute a gating signal, which is then applied to the irrep using element-wise multiplication.

Capsule networks use a special type of activation function called the squashing function <cit.>. The squashing function scales the magnitude of the output vectors (irreps) to be between 0 and 1 while preserving their direction.

Tensor Product (TP) activation is a learnable activation function that operates on the tensor product of irreps <cit.>. It applies a learnable set of weights to the tensor product of the input irreps and then projects the result back onto the original irrep basis.


The listed activation functions introduce non-linearity in equivariant neural networks. Unlike traditional interpretations, we view irreducible representations (irreps) as Fourier coefficients of the SO(3) space, offering a distinct perspective on their role in the network.

From this viewpoint,
we can classify activation functions as either global or local. 
Let f_in and f_out represent the input and output functions of an activation operation σ. We call an activation global if

    f_out = σ(f_in),

meaning that the value at any point in f_out depends on the values at all points in f_in. Conversely, an activation is local if 

    f_out(x) = σ(f_in(x)),


for any point x. This implies that the value at any point of the output function depends only on the value at the same point in the input function, establishing a direct and unambiguous mapping between input and output values in real space.


To the best of our knowledge, all the previously published activation methods in this domain are global.

This is suboptimal and can lead to noncompact representations and a lack of feature hierarchy learning, which our approach aims to address. By integrating Wigner matrices within a Fourier function framework in SO(3), our method maps values unambiguously in real rotational space, aiming to approximate the ReLU operator in real space. 

We further support the significance of locality with computational experiments, showing that models with local mapping significantly outperform those with global mapping in terms of accuracy.











§ THEORETICAL FRAMEWORK ON SPHERICAL HARMONICS, WIGNER, AND CLEBSCH-GORDAN COEFFICIENTS




This section provides an overview of spherical harmonics, Wigner coefficients, and Clebsch-Gordan coefficients, essential mathematical tools in the field of quantum mechanics for the analysis of angular momentum.



 §.§ Spherical Harmonics


Real spherical harmonics Y_l^m(Ω) : S^2→ℝ are a set of orthogonal functions defined on the surface of a sphere, denoted as S^2. They are useful in expanding functions defined over the sphere and appear extensively in the solution of partial differential equations in spherical coordinates.

The orthogonality of spherical harmonics is expressed as


    ∫_S^2  dΩ Y_l^k(Ω) Y_l'^k'(Ω) = δ_l l'δ_k k',


where δ stand for the Kronecker delta, signifying that spherical harmonics are orthogonal with respect to both the degree l and the order m.  Any square-integrable function f(Ω): S^2→ℝ can be decomposed into real spherical harmonics as

    f(Ω) = ∑_l = 0^∞∑_k = -l^l f^k_l Y^k_l(Ω),

where 

    f^k_l = ∫ 
        f(Ω)
        Y^k_l(Ω)
         d Ω

are the spherical harmonic expansion coefficients.
In practice, we use a fixed maximum expansion order L defined by the resolution of input data.



 §.§ Wigner Coefficients


Wigner matrices, denoted as D^l_mm' (ℛ), form  the irreducible representations of SO(3), the group of rotations in three-dimensional Euclidean space. These matrices are rotation operations for spherical harmonics, 


    Y_l^k_1(ℛΩ) = ∑_k_2 = -l^l D^l_k_1 k_2(ℛ) Y_l^k_2(Ω),


highlighting the transformation properties of spherical harmonics under rotation. The orthogonality of Wigner coefficients is held due to the following property,


    ∫_SO(3) d ℛ D^l_k_1 k_2 (ℛ) D^l'_k'_1 k'_2 (ℛ)  = 8 π^2/2 l +1δ_l l'δ_k_1 k'_1δ_k_2 k'_2.

Another property of a Wigner matrix is its unitarity,

    D^l_k_1 k_2 (ℛ^-1) = D^l_k_2 k_1 (ℛ).



Any square-integrable 
function f(ℛ) : SO(3) →ℝ can be decomposed into Wigner matrices as

    f(ℛ) = ∑_l = 0^∞∑_k_1 = -l^l∑_k_2 = -l^l f^l_k_1 k_2 D^l_k_1 k_2(ℛ),

where 

    f^l_k_1 k_2 = 2 l +1 /8 π^2∫  f(ℛ) D^l_k_1 k_2(ℛ) 
        d ℛ

are Wigner matrix decomposition coefficients.
Again, for practical considerations in Eq.  <ref> we will use a fixed maximum expansion order L .



 §.§ Clebsch-Gordan Coefficients


The Clebsch-Gordan coefficients facilitate the coupling of two angular momenta in quantum mechanics, leading to composite states with well-defined total angular momentum. 
The interconnection between Clebsch-Gordan coefficients, spherical harmonics, and Wigner matrices is illustrated through the integration of products of spherical harmonics and the transformation properties under rotation:


    ∫_S^2 d Ω Y^K_L(Ω) Y^k_l(Ω) Y^k'_l'(Ω) =   √((2 l + 1)(2 l' + 1)/4 π (2L+1))⟨ l 0 l' 0 | L 0⟩⟨ l k l' k' | L K ⟩



    ∫_SO(3) d ℛ D^l_k_1 k_2 (ℛ) D^l'_k'_1 k'_2 (ℛ)  D^L_K_1 K_2 (ℛ) =  8 π^2/2 L +1⟨ L K_1 | l k_1 l' k'_1 ⟩⟨ L K_2 | l k_2 l' k'_2  ⟩






§ SE(3) CONVOLUTION




 §.§ Motivation



<cit.> proposed a method for invariant pattern recognition in 3D CNNs. In each convolutional layer, they analyzed how the output varies with the orientation of the pattern. The authors examined a continuous distribution of the output in the rotation space at each voxel, where they conducted sampling and extracted the maximum value from the sampled points.

The method is equivariant with respect to the orientation of the input data and invariant to the orientation of data patterns. 
Formally, the method uses a convolutional operator that can be written down as

    h(r⃗) = max_ℛ∫_ℝ^3 f(r⃗ + r⃗_0) w(ℛr⃗_0) d r⃗_0,

where f(r⃗) is the input map and w(r⃗) is a trainable filter. In practice, ℝ^3 is discretized and the integration is replaced by a sum:

    h(r⃗) = max_ℛ∑_r⃗_0 f(r⃗ + r⃗_0) w(ℛr⃗_0).


However, the invariant pattern recognition approach is not the most expressive, as it cannot distinguish between different inputs. 
Figure <ref>A schematically demonstrates a toy example where a layer with invariant pattern recognition outputs identical feature maps for smiling and sad faces looking in different directions. 
A straightforward solution would be to retain not only the maximum output value  over the SO(3) rotation space, but also the distribution over the rotations. 
In other words, the result of such a convolutional operator will be a map in 6D: ℝ^3 × SO(3),

    h(r⃗,ℛ) = ∫_ℝ^3 f(r⃗ + r⃗_0) w(ℛr⃗_0) d r⃗_0,

where f(r⃗) and w(r⃗) are the input and the filter maps, respectively, or in a discrete form:

    h(r⃗,ℛ) = ∑_r⃗_0 f(r⃗ + r⃗_0) w(ℛr⃗_0).

From here on, without loss of generality, for brevity, we will only  use the integral form.
This operation is equivariant to both orientations of f(r⃗) and w(r⃗).


If f(r⃗) → f(ℛ_0r⃗), then

    h(r⃗,ℛ) →∫_ℝ^3 f(ℛ_0(r⃗ + r⃗_0)) w(ℛr⃗_0) d r⃗_0 = ∫_ℝ^3 f(ℛ_0r⃗ + r⃗_0) w(ℛℛ_0^-1r⃗_0) d r⃗_0 = h(ℛ_0r⃗,ℛℛ_0^-1).

Consequently, if w(r⃗) → w(ℛ_0r⃗), then

    h(r⃗,ℛ) →∫_ℝ^3 f(r⃗ + r⃗_0) w(ℛℛ_0 r⃗_0) d r⃗_0 = h(r⃗,ℛℛ_0).


It is worth noticing that the equivariant property in the two cases above holds in different manners. 
In the first case, both arguments r⃗ and ℛ are rotated. In the second case, rotation ℛ_0 is applied only to the second argument.  
Besides, the expression in Eq. <ref> shows how to coordinate rotation of two components of arguments of a 6D map.
Introduction of such a 6-dimensional feature map leads to the question how to treat such data. In this paper, we present  a novel equivariant convolution in 6D where both input and filter maps are in ℝ^3 × SO(3) and the set of rotations is continuous,


    h(r⃗, ℛ) = ∫_SO(3)dℛ_0 ∫_ℝ^3 dr⃗_0 f(r⃗ + r⃗_0, ℛ_0) w(ℛr⃗_0, ℛ_0 ℛ^-1)
        .

The usage of a continuous space of rotations ensures analytical equivariance of the convolution. Below we will show how the irreducible representation helps to perform an integration in the SO(3) rotation space.



 §.§ Integration in the SO(3) rotation space 





Let us substitute expansion coefficients from Eq. <ref>
into the convolution operation  h(r⃗, ℛ) in Eq. <ref>
and compute its expansion coefficients:


    h^l_1_k_1 k_2(r⃗) = 2 l_1 + 1/8 π^2∫_SO(3) dℛ  h(r⃗, ℛ) D^l_k_1 k_2 (ℛ) 
     = 2 l_1 + 1/8 π^2∫_SO(3) dℛ   D^l_1_k_1 k_2 (ℛ) ∫_SO(3)dℛ_0 ∫_ℝ^3 dr⃗_0 f(r⃗ + r⃗_0, ℛ_0) w(ℛr⃗_0, ℛ_0 ℛ^-1).

Let us now

decompose functions f(r⃗) and w(r⃗) using Eq. <ref>, Eq. <ref> and  Eq. <ref>,

    f(r⃗ + r⃗_0, ℛ_0) = ∑_l_2 = 0^L_in∑_k_3 = -l_2^l_2∑_k_4 = -l_2^l_2 f^l_2_k_3 k_4(r⃗+ r⃗_0) D^l_2_k_3 k_4 (ℛ_0),

and 

    w(ℛr⃗_0, ℛ_0 ℛ^-1) = ∑_l_3 = 0^L_in∑_k_5 = -l_2^l_2∑_k_6 = -l_2^l_2∑_k_7 = -l_2^l_2∑_l_4 = 0^L_filter∑_k_8 = -l_4^l_4∑_k_9 = -l_4^l_4 w^l_3 l_4_k_5 k_7 k_8(r) D^l_3_k_5 k_6 (ℛ_0)  D^l_3_k_7 k_6 (ℛ) D^l_4_k_8 k_9 (ℛ) Y_l_4^k_9 (Ω_r⃗_0),

where L_in and L_filter are 
the maximum expansion orders

of the filter in the rotational and the three-dimensional Euclidean spaces, respectively.Eq. <ref> also uses the unitarity property from Eq. <ref>.

Eq. <ref> then simplifies to

    h^l_1_k_1 k_2(r⃗) = ∫_ℝ^3 dr⃗_0 ∑_l_2 = 0^L_in∑_k_3 = -l_2^l_2∑_k_4 = -l_2^l_2 f^l_2_k_3 k_4(r⃗+ r⃗_0) S^l_1 l_2_k_1 k_2 k_3 k_4 (r⃗_0),

where 

    S^l_1 l_2_k_1 k_2 k_3 k_4 (r⃗_0) = 8 π^2/ 2 l_2 + 1∑_l_4 = 0^L_filter( ∑_k_7 = -l_2^l_2∑_k_8 = - l_4^l_4⟨ l_1 k_1 | l_2 k_7 l_4 k_8 ⟩ w^l_2 l_4_k_3 k_7 k_8(r_0) ) ( ∑_k_9 = - l_4^l_4⟨ l_1 k_2 | l_2 k_4 l_4 k_9 ⟩ Y^k_9_l_4 (Ω_r⃗_0)).

Appendix <ref> presents a more detailed derivation of Eq. <ref>. 

Parameters w^l_2 l_4_k_3 k_7 k_8(r_0) are trainable.

The output map expansion coefficients h^l_1_k_1 k_2(r⃗) have the maximum degree of L_out. 

We shall also stress that the

values of L_filter,L_in and L_out must satisfy the triangular inequality,

    L_in - L_filter≤ L_out≤ L_in + L_filter.


The computational complexity of the convolution in Eq. <ref> is O(N^3 L_filter^3 L_in^3 L_out^3 D_in D_out), where N is the linear size of the input data, and D_in and D_out are numbers of input and output channels. 

Appendix <ref> 

proves
the roto-translational equivariance with respect to the input data and the rotational equivariance with respect to the filter.





§ NON-LINEAR OPERATIONS



 §.§ The need for local activation function


In neural networks, linear operations typically alternate with nonlinear ones; the latter are often referred to as activations. Nonlinear operations are necessary to approximate complex dependencies in the data. 
However, from the perspective of learning hierarchies of patterns,  the two types of operations have different purposes.
Indeed, linear layers are tasked with detecting patterns and fetching a quantity proportional to the probability of finding a certain pattern. 

Nonlinear layers, in turn, penalize low pattern probabilities, allowing only high ones to pass and often serve to learn compact representations. 
An incorrectly chosen activation can lead to noise accumulation when passing through multiple layers and also suboptimal latent representations. 

A classical activation solution 

is ReLU and its variants, which, despite some drawbacks, are the most intuitively understandable activation functions. 
The vanilla ReLU operator passes only positive values of its input. 
In convolutional and other networks dealing with spatial data, the critical property that differentiates activation operations 
from the rest of the network is their locality,
i.e., the activation is applied to different points in the Euclidean space independently. 





In our case, as we have already mentioned, the output of the linear convolution is six-dimensional: ℝ^3 ×SO(3). 
Since the ℝ^3 dimension discretizes into voxels, we apply the activation to each voxel separately. 
The most interesting question, however, is how to design the activation function for the SO(3) data dimension.

We define the distribution of values in SO(3) in our architecture with 
Wigner matrix decomposition coefficients.

Each coefficient, similarly to the Fourier series, represents the global properties of the entire SO(3) distribution rather than an individual point in SO(3). 
However, we aim to design a local activation for SO(3)
in the space of Wigner coefficients. 






A previously used approach <cit.> consists in sampling the SO(3) space and requires a transformation from the 
Wigner
representation into the 
rotation

space, an activation in the 
SO(3)

space, and an inverse transformation into the 
Wigner

space.



However,
this approach loses the analytical SO(3) equivariance.
It may not be a problem since sampling a sufficiently large number of points in the rotation SO(3) space can guarantee effective equivariance. 
However, this strategy
inevitably leads to an increased computational complexity. 
Conversely, the proposed Wigner coefficient representation,
allows using few values to define a function in SO(3) and retaining analytical equivariance at the same time.
There are no analytical expressions for the decomposition coefficients of ReLU(f(ℛ)) given the spherical harmonics coefficients of f(ℛ). 
However, we will use the analytical expression of the coefficients of a product of two functions, 

given in Eq. <ref> in Appendix <ref>.

It enables us to find an analytical expression for the spherical harmonics coefficients for a polynomial applied to a function in SO(3). 





 §.§ Local activation in the Wigner space



It has been demonstrated that in NN architectures,
on a certain interval, the ReLU function can be approximated by a polynomial of a low degree <cit.>. 
Even the quadratic polynomial of degree two showed competitive results in the ResNet architectures <cit.>.
However, as mentioned in Appendix <ref>, the product of two functions defined in the SO(3) space with Wigner coefficients, has a higher resolution 
than each of the initial functions.

Subsequently, to avoid the loss of information, we need to increase the maximum expansion order of the product, given as the Wigner matrix expansion. 
For a function in SO(3) 
defined with Wigner coefficients
of the maximum degree L, the result of applying a polynomial of degree n will have the maximum expansion order of nL-1. 
Since the number of 
Wigner
expansion coefficients grows as the cube of the expansion order, to approximate the activation function, we chose the activation polynomial of the second degree. 

This paper studies multiple approximation strategies 
to this polynomial in the 
SO(3) space.

Generally, the ReLU approximation expression can be written down as

    f_act(r⃗) = D  P_2(f(r⃗)/D),

 where P_2(x) = c_0 + c_1 x + c_2 x^2 and D is a scaling factor. 
Below, we discuss several approaches for this approximation in the Wigner space.

 

 §.§ Adaptive coefficients

 
In the "adaptive coefficients" approach, we approximate the ReLU operator on a [-3 σ(f) + μ(f), 3σ(f) + μ(f)] interval, where
 (μ, σ) are the mean and the standard deviation of f(ℛ), respectively,

    μ = ∫_SO(3)   f(ℛ) dℛ/∫_SO(3) dℛ≡ f^0_00,



    σ = √(∫_SO(3) (f(ℛ)-μ)^2  d ℛ/∫_SO(3) d ℛ)≡√(∑_l = 1^L∑_k_1, k_2 = -l^l(f^l_k_1 k_2)^2/2l + 1).


There are three possible cases to consider. 
In the special case of 3 σ + μ < 0, we make an assumption  that ∀ℛ f(ℛ) < 0. 
Then, we put the polynomial approximation to y = 0.01 x, instead of  y = 0 to avoid zero gradients. 
In the second special case of -3 σ + μ > 0, 
we assume ∀ℛ f(ℛ) > 0 and set y = x. 
In the general case, the function f(ℛ) ranges both positive and negative values. 
For simplicity, we first divide the values of the function by 3 σ and then apply a polynomial that approximates ReLU on [-1 + k, 1 + k], where k = μ/3 σ, k ∈ [-1,1]. 
To find the optimal polynomial coefficients in this case, we formulate the following minimization problem,

    min_c_0(k), c_1(k), c_2(k) F(c_0, c_1, c_2, k),

where 

    F(c_0, c_1, c_2, k)  = ∫_-1 + k^0 (c_0 + c_1 x + c_2 x^2 )^2 d x + ∫_0^1 + k (c_0 + c_1x + c_2 x^2 - x)^2 d x.






Then, we can consistently deduce (see details in Appendix <ref>)

    c_2 = 15/32(k^4 - 2k^2 + 1)
    
            c_1 = 1/16 (-15 k^5 + 26 k^3 - 3k + 8)
    
            c_0 = 3/32 (5 k^6 - 9 k^4 + 3k^2 + 1)
        .



Figure <ref>  in Appendix <ref> shows the plot of these coefficients as a function of normalized mean k and the error of this approximation, respectively. 

Finally, to obtain the polynomial activation, we multiply the result by 3 σ,

    f_act(ℛ) = 3 σ P_2( f(ℛ)/3 σ).





 §.§ Constant coefficients


In the second strategy, 
we consider the
denominator D  equal to the third part of the function norm f_2. In practice, this means that the function range lies in the [-1, 1] interval. This is a special case of the interval from the previous approach at k = 0. Such a value of k leads to the following polynomial coefficient values,

    c_2 = 15/32; c_1 = 1/2; c_0 = 3/32→  P_2 = 3/32 + 1/2x + 15/32 x^2.

We then multiply the result of the polynomial function by  f_2/3,

    f_act(ℛ) = (f_2/3) P_2( f(ℛ)/(f_2/3)).





 §.§ Trainable coefficients


We have also studied an approach with trainable polynomial coefficients.
Here, the denominator D is the same as in the previous case, D = f_2/3. 
However, the 
three

polynomial coefficients  (c_0,c_1,c_2) are trainable values. 










 §.§ Global activation





Above, we introduced a local activation function that  requires a higher resolution and, consequently, an increased number of coefficients. This raises the question of whether the performance improvement offered by this operation justifies the associated rise in the computational cost. As a baseline for the activation, we explored an operation where the output resolution is identical to that of the input.

The formalism presented above  allows us to define the global activation,


    f_GlobAct(ℛ) = GlobAct(f(ℛ)) = σ(WSoftMax(f)+b) f(ℛ)
    ,


where we use SoftMax from Eq. <ref>, apply σ = sigmoid, and then multiply the result by the function f(ℛ). W and b are trainable parameters.
Linearity of the Wigner matrix decomposition gives us the following expression for the Wigner coefficients of the global activation function,

    [f_GlobAct]^l_k_1 k_2 =  σ(WSoftMax(f)+b) [f]^l_k_1 k_2.


This expression can be seen as a gating mechanism that
only passes SO(3)-distributions  with sufficiently high trainable positive values, at least for a single point in SO(3). Although the input and output resolutions of this operation are identical, we still employ a higher resolution (expansion order) in the ReLU approximation within the SoftMax operation.  





 §.§ Normalization

The fact that the mean and the standard deviation of a function in the rotational space can be expressed in terms of Wigner coefficients allows us to introduce an expression for the normalization of a function in SO(3):

    f_n(ℛ) = γf(ℛ) - μ/σ + β,

where μ and σ are the mean and the standard deviation of function f, respectively, and γ and β are trainable coefficients. 
We can also extend this operation for data in ℝ^3×SO(3), where the Euclidean component is discretized and characterized by a set of functions f_ijk(ℛ), where i, j, and k are voxel indices. The batch normalization is then defined as follows:

    [f_ijk]_n(ℛ) = γf_ijk(ℛ) - μ_ijk/σ_ijk + β,

where μ_ijk and σ_ijk are the mean and the standard deviation of a function f_ijk, respectively, and μ = 1/N∑_i,j,kμ_ijk and σ^2 = 1/N∑_i,j,kσ^2_ijk, with N being the number of voxels.








 §.§  Max pooling operation in the continuous SO(3) space






For a wide range of tasks, whether they involve global or local (voxel-wise) prediction, it is necessary to design an operation that reduces the representation of the function in the rotation space to a single value. 
For brevity, we may call it a pooling operation in SO(3). 
The most straightforward approach, average pooling of a function in the SO(3) space, entails discarding 
Wigner
coefficients 
of all degrees higher than zero. 
However, we aim to design an operation analogous to max pooling. 


Let us note that the polynomial approximation 
from Eq. <ref>
allows for the simulation  of the SoftMax effect in SO(3). 

More precisely, we can define SoftMax as

    SoftMax(f) = ∫_SO(3) w(ℛ) ReLU(f(ℛ)) dℛ,


where weight w(ℛ) = ReLU(f(ℛ))/∫_SO(3)ReLU(f(ℛ'))  dℛ'. 
The weight function estimates the ratio of ReLU(f(ℛ)) to the positive part of the function f(ℛ). 


Thus, we can also define the SoftMax function in SO(3) using 
only the 
Wigner matrix expansion coefficients
of f(ℛ)
as


    SoftMax_poly(f) =  
    
    ∫_SO(3)act(f(ℛ))^2 dℛ/∫_SO(3)act(f(ℛ)) dℛ = f_a^2_2/f_a_1 = ∑_l', k'_1, k'_28 π^2/2 l'+1 ([f_a]^l_k_1 k_2)^2/[f_a]^0_0 0.


We should note that there can be different approaches to SoftMax simulation, even with the usage of polynomial coefficients. 

In scenarios where the pooling operation follows the ReLU operation, there is no need to apply the activation function again. Instead, we can directly calculate the ratio of the 2-norm squared of the input function to the pooling operation to its 1-norm.














§ DATASETS, ARCHITECTURES AND TECHNICAL DETAILS



Our method is centered on applications involving regular volumetric data. Extending this method to irregular data would necessitate significant modifications that are beyond the scope of this paper. Consequently, we evaluated our method using a collection of voxelized 3D image datasets.




 §.§ Datasets


We assessed our designs on MedMNIST v2, a vast MNIST-like  collection of standardized biomedical images <cit.>. 
It is designed to support a variety of tasks, including binary and multi-class classification, and ordinal regression.
The dataset encompasses 

six sets comprising a total of 9,998 3D images.
All images are resized 

to 

28 × 28 × 28 voxels, each paired with its respective classification label. 



We used the train-validation-text split provided by the authors of the dataset (the proportion is 7:1:2).



 §.§ Baseline architectures


For our baselines, we employed the same model configurations that the creators of the dataset
utilized for testing on MedMNIST3D datasets<cit.>. 
These include various adaptations of ResNet <cit.>, featuring 2.5D/3D/ACS <cit.> convolutional layers, alongside open-source AutoML solutions such as auto-sklearn <cit.>, and AutoKeras <cit.>. 
We have also added results of the models that were tested on the collection: FPVT <cit.>, Moving Frame Net <cit.>, Regular SE(3) convolution <cit.> and ILPOResNet50 (<cit.>).


Table <ref> lists all the tested architectures.



 §.§ Trained models

We developed multiple architectures with various activation methods, each reflecting the layer sequence of ResNet-18 <cit.>. 

The final architecture, featuring reduced-resolution activation, is inspired by ResNet-50 <cit.>. 


Figure <ref> schematically illustrates our EquiLoPo architecture along with its main components. These include the initial block, the repetitive building block, pooling operators in SO(3) and 3D spaces, and a linear transformation at the end.
Table <ref> details the basic block for SE(3) data in ResNet-18.

The initial convolutional block, outlined in Table <ref>, initiates the architecture.
We also specifically adapted the Batch Normalization process for 6D (3D ×SO(3)) data.








It is important to note that the 

Bottleneck  block 

in ResNet-50 (the simplest building block of the architecture), which increases activation resolution, results in an eightfold increase in the maximum degree. This surge is attributed to triple consecutive activations without an intervening convolution, potentially escalating computational demands. Consequently, our exploration was limited to the ResNet-18 architecture, where the activation operation enhances output resolution. Figure <ref> illustrates the EquiLoPO ResNet-18 architecture.










Every filter in our trained networks is confined within a 3 × 3 × 3 volume in 3D. Filters are parameterized by weights w^l_2 l_4_k_3 k_7 k_8(r_i), where r_i ∈0,1,√(2), √(3) represents all possible radii from the center in the cubic space.


We used the ADAM optimizer <cit.>. In order to avoid overfitting, we apply the dropout operation after each activation operation.
The learning rate of the optimizer and the dropout rate are hyperparameters. 
All the models are trained for 100 epochs.

Table <ref> in Appendix <ref> lists hyperparameters 
optimized for validation data performance. Table <ref> presents memory and time consumption metrics for EquiLoPOResNet-18, ILPOResNet-18 and ResNet-18 models in the inference mode. 

The current implementation consumes significantly more memory than the vanilla ResNet architecture and requires longer execution time because the activation is coded in Python rather than C++. Transitioning the implementation to C++ could significantly reduce both memory and CPU footprints.






    







§ RESULTS





 §.§ The need for higher resolution in the activation function


We previously highlighted the necessity of ensuring that the activation coefficients of the output possess a maximum degree higher than that of the input to prevent information loss when applying the polynomial, 
specifically,  L_out = 2L_in. 
However, investigating the scenario where high output degrees are discarded, aligning the maximum output degree with the maximum input degree ( L_out = L_in) merits consideration.

Figure <ref> illustrates the impact of the increased resolution on the approximation accuracy of the SO(3) ReLU function. For this plot, we analyzed one of the rotational distributions from

the last convolution of 
the trained model  with local adaptive coefficients on the Vessel dataset reported in Table <ref>
(L_in = 2) to observe how the activation operator, with constant (subplots A-B) and adaptive (subplots C-D) polynomial coefficients, operates under two conditions:  
L_out = L_in and L_out = 2L_in. We sampled 10^6 points in SO(3) and computed the values of both input and output functions at these points, presenting the results on a density plot.


Subplots A and C offer illustrations for the scenario with a reduced resolution (L_out = L_in). Here, we note an inherent ambiguity between the values of the input and output signals. 
Indeed, each input function value maps to multiple output values.
Such an inadequate mapping is also
evidenced by
the absence of a sharp negative value cutoff.

Subplots B and D showcase the polynomial approximation with an enhanced output resolution. 
This setup eliminates the ambiguity in mapping in the real space, as each input function value uniquely corresponds to a single output function value. 
Moreover, this operator exhibits nonlinearity and closely approximates the ReLU function, depicted by a black dashed line. 
Consequently, negative input values transform into near-zero output values, a phenomenon clearly visible in the right-hand histogram, which features a prominent peak near y=0 and a sharp cutoff of negative values.



The extent to which increasing the resolution enhances the performance of the architecture is evident in Table <ref>. 
We conducted additional tests using a single architecture with a reduced resolution. 

This architecture employs local activation functions with trainable polynomial coefficients. 
As we can wee from the table, the model with reduced resolution consistently underperforms across all datasets compared to the full-resolution model that also utilizes locally trainable coefficients.










 §.§ Constant and adaptive polynomial coefficients strategies in the activation



In the case when the input function's  mean is close to zero, the adaptive coefficients are almost identical to the constant ones and the two strategies have very similar effect as suggested by Eq. <ref>.


In practice, zero-mean functions are frequent

but there are also cases when the mean is shifted. 
Figure <ref> demonstrates how the activation acts on a function with the mean significantly shifted to the negative values. 

Upper subplots of the four pictures demonstrate distribution of the input function defined in SO(3). The center of the distribution is clearly shifted to the left.
In this case, the polynomial with adaptive coefficients (subplot D) approximates the ReLU function more precisely than the polynomial with constant coefficients (subplot B) 
because the approximation is closer to the ReLU function.
Indeed, according to Eq. <ref>,
the adaptive coefficients strategy gives twofold less approximation error than the constant coefficients one.

Table <ref> presents the performance metrics for both strategies across the dataset collection. The choice between local and global activation significantly influences the performance. 
Therefore, our comparative analysis focuses more on contrasting local activations with adaptive versus constant polynomial coefficients. The strategy employing adaptive coefficients outperforms the one with constant coefficients on all the datasets, except for Adrenal. This outcome is expected, as adaptive coefficients typically provide a superior approximation of the activation function.



 §.§ Local and global activation


We evaluated global and local activation layers across three distinct strategies for the coefficients of the approximating activation polynomial. 
Table <ref> indicates that for the trainable coefficients, the performance of global activation is comparable to that of the local activation. However, for both constant and adaptive coefficients, the local activation demonstrates superior results across all datasets.




 §.§ Analysis of the SoftMax pooling in SO(3)


To examine the behavior of the SoftMax 

operation, we analyzed the outputs from the last convolution layer of the trained model with locally adaptive coefficients on the Vessel dataset reported in Table <ref>. 
We extracted the SO(3) function coefficients for all the maps, voxels, and channels, and then computed the SoftMax for these functions using Eq. <ref>. 
We sampled 10^6 points in the SO(3) space and reconstructed the functions at these points. 
Figure <ref> demonstrates the relationship between SoftMax and the sampled maximum. 
It is evident that there is high 
correlation between the two functions in the positive region of the sampled maximum. 

Negative values in the sampled maximum lead to SoftMax values of zero. Near zero, uncertainty increases along with a wide range of SoftMax values. This phenomenon appears due to the polynomial approximation error. 
Specifically, the accuracy of the SoftMax approximation hinges on the condition that the mean of the 
output of the ReLU operation approximation
must be  greater or equal to zero. Additionally, the mean of the SO(3) distribution, or its L_1-norm, should be substantially smaller than the  L_2-norm. While this condition holds in most instances, there are exceptions where it does not.
These cases are clearly visible in Figure <ref> as those providing high uncertainty in the SoftMax estimation near zero sampled maximum.















 §.§ Comparison with other methods




Table <ref> lists a detailed comparison of our models' performance against the baseline models: various adaptations of ResNet <cit.>, featuring 2.5D/3D/ACS <cit.> convolutional layers, alongside open-source AutoML solutions such as auto-sklearn <cit.>, and AutoKeras <cit.>. 
We have also added results of the models that were tested on the collection: FPVT <cit.>, Moving Frame Net <cit.>, Regular SE(3) convolution <cit.> and ILPOResNet50 (<cit.>).

In this context, the classes are imbalanced, meaning that to measure the predictive precision of a model, we need other metrics besides accuracy (ACC).  Consequently, we also consider the AUC-ROC (AUC) metric, which provides a more informative assessment of  imbalanced datasets.

According to the metrics, we can observe that our models either outperform or demonstrate state-of-the-art levels rounded to the third significant digit
on nearly all metrics, except for the accuracy on the organ dataset. The organ dataset is largely resolved, with all tested methods achieving a very high level of AUC. Additionally, the dataset might benefit from more channels (more than the four we used) within the model's layers for improved performance. However, one of our goals was to limit the number of parameters in our model.


Among the state-of-the-art methods, only the one proposed by <cit.> has fewer parameters. The ILPOResNet model is based on the same principles of pattern definition and detection but carries out recognition in an invariant manner while remaining in 3D. In contrast, the presented equivariant model offers a more expressive architecture, as evidenced by its performance on the collection datasets, but requires more parameters.

When considering the various types of activation we tested, three stand out with the best performance: local activation with trainable coefficients, local activation with adaptive coefficients, and global activation with trainable coefficients. Tables <ref>, <ref>, and <ref> in Appendix <ref>  show the results of experiments on architectures with these types of activations but a smaller number of building blocks (1 or 2 instead of 8). Local activation with adaptive coefficients demonstrates the highest robustness since, even in small models, it performs relatively well.


Here, a reasonable question may arise whether the local activation is necessary if the global activation already demonstrates comparable results. We shall notice, however, that global activation, in the current implementation, employs SoftMax in the argument of the multiplier function, which in turn utilizes the ReLU approximation, a form of local activation. 
In Appendix <ref>, we tested an alternative strategy for global activation, replacing SoftMax with the 2-norm. As evidenced by the results presented in Table <ref>, this alternative strategy yields significantly inferior performance compared to the original approach with SoftMax.  











§ CONCLUSION


This work introduces a novel equivariant neural network architecture that achieves analytical rotational equivariance on the continuous SO(3) group while retaining the flexibility of unconstrained trainable filters. Our key innovations are a group convolutional operation that leverages irreducible representations as the Fourier basis and a local activation function in the SO(3) space that provides a well-defined mapping from input to output function values in the real space.

By employing these operations within ResNet-style architectures, we created an equivariant model that overcomes the limitations of existing group convolutional networks, restricted to discrete rotation groups and steerable convolutional networks, which impose constraints on the filter structures.

We thoroughly evaluated our model on the 
diverse
MedMNIST3D collection of 3D medical imaging datasets. Across multiple datasets, our equivariant architecture consistently outperformed state-of-the-art methods, demonstrating its effectiveness in analyzing volumetric data while respecting rotational symmetries.

The superior performance of our model highlights the importance of true rotational equivariance on the continuous SO(3) group, as well as the benefits of a local activation function that preserves equivariance properties while maintaining the flexibility of unconstrained filters.

Our work provides a powerful and flexible framework for equivariant deep learning on volumetric data, paving the way for future research in this domain. Potential future directions include exploring opportunities for the local activation operation within other networks using irreducible representations in 
various spaces. 



§ ACKNOWLEDGEMENTS

This work is partially funded by MIAI@Grenoble Alpes (ANR-19-P3IA-0003).



    
plainnat






§ FILTER EXPRESSION


In this section, we reveal how the expression in Eq. <ref> is obtained:

    S^l_1 l_2_k_1 k_2 k_3 k_4 (r⃗_0) =  2 l_1 + 1/8 π^2∑_l_3 = 0^L_in∑_k_5 = -l_3^l_3∑_k_6 = -l_3^l_3∑_k_7 = -l_3^l_3∑_l_4 = 0^L_filter∑_k_8 = - l_4^l_4∑_k_9 = - l_4^l_4( ∫_SO(3)  d ℛ_0 D^l_2_k_3 k_4 (ℛ_0)D^l_3_k_5 k_6 (ℛ_0) ) 
    ( ∫_SO(3) d ℛ  D^l_1_k_1 k_2(ℛ) D^l_3_k_7 k_6 (ℛ) D^l_4_k_8 k_9 (ℛ)  ) w^l_3 l_4_k_5 k_7 k_8(r) Y^k_9_l_4(Ω_r⃗)  =  2 l_1 + 1/8 π^28 π^2/ 2 l_2 + 18 π^2/ 2 l_1 + 1
    ∑_l_3 = 0^L_in∑_k_5 = -l_3^l_3∑_k_6 = -l_3^l_3∑_k_7 = -l_3^l_3∑_l_4 = 0^L_filter∑_k_8 = - l_4^l_4∑_k_9 = - l_4^l_4δ_l_2 l_3δ_k_3 k_5δ_k_4 k_6⟨ l_1 k_1 | l_3 k_7 l_4 k_8 ⟩⟨ l_1 k_2 | l_3 k_6 l_4 k_9⟩ w^l_3 l_4_k_5 k_7 k_8(r) Y^k_9_l_4(Ω_r⃗)
    =8 π^2/ 2 l_2 + 1∑_l_4 = 0^L_filter( ∑_k_7 = -l_2^l_2∑_k_8 = - l_4^l_4⟨ l_1 k_1 | l_2 k_7 l_4 k_8 ⟩ w^l_2 l_4_k_3 k_7 k_8(r_0) ) ( ∑_k_9 = - l_4^l_4⟨ l_1 k_2 | l_2 k_4 l_4 k_9 ⟩ Y^k_9_l_4 (Ω_r⃗_0))




§ EQUIVARIANCE OF THE SE(3) CONVOLUTION





 §.§ Equivariance with respect to the input function

Let f(r⃗, ℛ) be rotated by ℛ_1. Then consider the output of the 
convolution in Eq. <ref>:


    ∫_SO(3)dℛ_0 ∫_ℝ^3 dr⃗_0 f(ℛ_1r⃗ + ℛ_1r⃗_0, ℛ_0ℛ_1^-1) w(ℛr⃗_0, ℛ_0 ℛ^-1) = ∫_SO(3)dℛ_0 ∫_ℝ^3 dr⃗_0 f(ℛ_1r⃗ + r⃗_0, ℛ_0ℛ_1^-1) w(ℛℛ_1^-1r⃗_0, ℛ_0 ℛ^-1) 
    
        = ∫_SO(3)dℛ_0 ∫_ℝ^3 dr⃗_0 f(ℛ_1r⃗ + r⃗_0, ℛ_0) w(ℛℛ_1^-1r⃗_0, ℛ_0 ℛ_1 ℛ^-1) = h(ℛ_1r⃗, ℛℛ_1^-1)




 §.§ Equivariance with respect to the filter function

Let w(r⃗, ℛ) be rotated by ℛ_1. Then consider the output of the convolution in Eq. <ref>:


    ∫_SO(3)dℛ_0 ∫_ℝ^3 dr⃗_0 f(r⃗ + r⃗_0, ℛ_0) w( ℛℛ_1 r⃗_0, ℛ_0 ℛ_1^-1ℛ^-1) =  h(r⃗, ℛℛ_1)





§ PRODUCT OF TWO FUNCTIONS IN SO(3)


 The Wigner matrix decomposition coefficients of a product of two functions defined in SO(3) can be expressed through coefficients of two multiplier functions. 


Let functions f_1(ℛ) and f_2(ℛ) have coefficients [f_1]^l_1_k_1 k_2 and [f_1]^l_1_k_1 k_2 and f_prod(ℛ) be the product of f_1(ℛ) and f_2(ℛ):

    f_prod(ℛ) = f_1(ℛ)f_2(ℛ).

Then

    [f_prod]^l_3_k_5 k_6 = ∑_l_1=0^L_1∑_l_2=0^L_2∑_k_1 = -l_1^l_1∑_k_2 = -l_1^l_1∑_k_3 = -l_2^l_2∑_k_4 = -l_2^l_2⟨ l_3 k_5| l_1 k_1 l_2 k_3 ⟩⟨ l_3 k_6| l_1 k_2 l_2 k_4 ⟩ [f_1]^l_1_k_1 k_2 [f_2]^l_2_k_3 k_4

where [f_prod]^l_3_k_5 k_6 are coefficients of f_prod(ℛ), L_1 and L_2 are maximal degrees of f_1 and f_2 respectively. In order to avoid loss of information the maximal degree of f_prod is L_3 = L_1 + L_2. Thus, the product operation requires increase of resolution of data. 




§ COMPARISON WITH STEERABLE NETWORKS AND IMPORTANCE OF ACTIVATION 


As we indicated in the main text, the EquiLoPo convolution is defined as


    h^l_1 d_out_k_1 k_2(r⃗) = ∫_ℝ^3 dr⃗_0 ∑_l_2 = 0^L_in∑_k_3 = -l_2^l_2∑_k_4 = -l_2^l_2 f^l_2 d_in_k_3 k_4(r⃗+ r⃗_0) [S_EquiLoPo]^l_1 l_2 d_in d_out_k_1 k_2 k_3 k_4 (r⃗_0),

where 

    [S_EquiLoPo]^l_1 l_2  d_in d_out_k_1 k_2 k_3 k_4 (r⃗_0) = 8 π^2/ 2 l_2 + 1∑_l_4 = 0^L_filter   [p_EquiLoPo]^l_1 l_2 l_4  d_in d_out_k_1 k_3(r_0) ( ∑_k_9 = - l_4^l_4⟨ l_1 k_2 | l_2 k_4 l_4 k_9 ⟩ Y^k_9_l_4 (Ω_r⃗_0)),

and

    [p_EquiLoPo]^l_1 l_2 l_4  d_in d_out_k_1 k_3(r_0) = ∑_k_7 = -l_2^l_2∑_k_8 = - l_4^l_4⟨ l_1 k_1 | l_2 k_7 l_4 k_8 ⟩ [w_EquiLoPo]^l_2 l_4  d_in d_out_k_3 k_7 k_8(r_0).

In our implementation, we set coefficients [w_EquiLoPo]^l_2 l_4  d_in d_out_k_3 k_7 k_8(r_0) to be trainable but it is equivalent to the case when  [p_EquiLoPo]^l_1 l_2 l_4  d_in d_out_k_1 k_3(r_0) is trainable because

    ∑_l_1 = |l_2 - l_4|^l_2+l_4∑_k_1 = - l_1^l_1⟨ l_1 k_1 | l_2 k_7 l_4 k_8 ⟩    [p_EquiLoPo]^l_1 l_2 l_4  d_in d_out_k_1 k_3(r_0) =   [w_EquiLoPo]^l_2 l_4  d_in d_out_k_3 k_7 k_8(r_0).

Steerable Networks in the 3-dimensional case <cit.> employ the following  convolution operator:

    h^l_1  d_out_k_2(r⃗) = ∫_ℝ^3 dr⃗_0 ∑_l_2 = 0^L_in∑_k_4 = -l_2^l_2 f^l_2  d_in_k_4(r⃗+ r⃗_0)  [S_Steerable]^l_1 l_2  d_in d_out_k_2 k_4 (r⃗_0),

where 

    [S_Steerable]^l_1 l_2  d_in d_out_k_2 k_4 (r⃗_0) = ∑_l_4 = 0^L_filter [p_Steerable]^l_1 l_2 l_4  d_in d_out(r_0)( ∑_k_9 = - l_4^l_4⟨ l_1 k_2 | l_2 k_4 l_4 k_9 ⟩ Y^k_9_l_4 (Ω_r⃗_0)),

and coefficients [p_Steerable]^l_1 l_2 l_4  d_in d_out(r_0) are trainable. 


If we consider convolution in isolation from other operations in a network, then our convolution operator is equivalent to the Steerable network convolution where the number of input and output degree features are D_in^l_1 = D_in (2l_1 + 1) and D_out^l_2 = D_out (2l_2 + 1) respectively:

    [p_Steerable]^l_1 l_2 l_4  (d_in(2 l_1 + 1) + k_1) (d_out(2 l_2 + 1) + k_1)(r_0) = [p_EquiLoPo]^l_1 l_2 l_4  d_in d_out_k_1 k_3(r_0).





The difference between our method and Steerable Networks lies in the interpretation of the coefficients, namely we treat them as expansion coefficients in SO(3). This interpretation results in a special approach to activation. We propose such operation on the expansion coefficients that would correspond to the ReLU operator in the rotational space.











§ DERIVATION OF EXPRESSIONS FOR ADAPTIVE COEFFICIENTS



The solution of Eq. <ref> satisfies the following conditions,

    ∂ F / ∂ c_0 =  4 c_0  + 4 c_1 k + 4 c_2 (k^2 + 1/3) - (k+1)^2 = 0 
    ∂ F / ∂ c_1 = 4 c_0 k +  4 c_1 (k^2 + 1/3) + 4 c_2 (k^3 + k) - 2/3 (k+1)^3 = 0 
    ∂ F / ∂ c_2 = 4 c_0 (k^2 + 1/3) +  4 c_1 (k^3 + k) + 4 c_2 (k^4 + 2 k^2 + 1/5) - 1/2 (k+1)^4 = 0 
    .

Multiplying the first equation by k and (k^2 + 1/3) and subtracting it from  the second  and third equations, respectively, we get

    c_1 + 2 c_2 k = -1/4 (k^3 - 3k -2) 
    
            c_1 k + 2 c_2 (k^2 + 1/15) = -1/16 (3k^4 - 10k^2 - 8k +1)
        .


Figure <ref> shows the plot of these coefficients as a function of normalized mean k and the error of this approximation, respectively. 












§ HYPERPARAMETERS



Table <ref> shows optimal hyperparameters for the trained networks. All the networks, except for the last one, follow the ResNet-18 architecture. The last model is based on the ResNet-50 version.







§ PERFORMANCE WITH FEWER BUILDING BLOCKS



In the main experiments, we evaluated the proposed neural network architecture using the standard ResNet-18 configuration with 8 building blocks. To further investigate the method's performance and scalability, we conducted additional experiments by reducing the number of building blocks to 1 and 2, respectively, for the three activation functions that demonstrated the best performance in the main experiments. These were local activation with trainable coefficients(Table <ref>), local activation with adaptive coefficients(Table <ref>), and global activation with trainable coefficients(Table <ref>).


As the number of building blocks decreases from 8 to 2 and then to 1, there is a general trend of degrading the performance across all datasets and activation functions. The extent of performance degradation varies across datasets. For example, the OrganMNIST3D dataset exhibits a more significant drop in accuracy (ACC) when reducing the number of blocks compared to other datasets like NoduleMNIST3D or VesselMNIST3D. 

The three activation functions (local activation with trainable coefficients, local activation with adaptive coefficients, and global activation with trainable coefficients) show different levels of robustness to the reduction in building blocks. Local activation with adaptive coefficients (Table 5) maintains relatively high performance even with fewer blocks, compared to the other two activation functions.











§ ALTERNATIVE STRATEGY FOR THE GLOBAL ACTIVATION



In this section, we test an alternative strategy for global activation. Instead of using Eq. <ref>, we use

    f_GlobAct(ℛ) = GlobAct(f(ℛ)) = σ(W|f|_2+b) f(ℛ)
    ,

where we replace the SoftMax function with the 2-norm of the function.

Table <ref> compares the performance of the ResNet-18-like architecture with global activation using trainable coefficients and the two different strategies: SoftMax and 2-norm. The results show that using the SoftMax function consistently outperforms the 2-norm alternative across all datasets.






