\section{Introduction}
Hadronic jets are the most abundant products of proton collisions at the LHC.
They appear as a result of hard QCD processes, as well as decay products in many rare signatures.
Moreover, many LHC analyses involve searching for highly boosted hadronically-decaying objects, in which two or more collimated decay products deposit energy into an single overlapping region of the detector~\cite{Butterworth:2008iy,Plehn:2009rk,Abdesselam:2010pt,Shimmin:2016vlc}.
The extremely large QCD cross sections often lead to a scenario where the rate of background jets overwhelm the signal.
Therefore, classifying the particle origin of jets is of special importance.

In general, discrimination of candidates in a particular kinematic phase space is based on features of the \textit{jet substructure}.
That is, we consider the detailed patterns of radiation observed via the tens to hundreds of consitutent measurements that comprise each jet.
For decades, there has been much progress in developing analytical variables to encode distinctive structural features, such as jet shape, n-subjettiness, energy correlation functions, \textit{etc}.~\cite{Gallicchio:2012ez,Komiske:2017aww,Thaler:2010tr,Larkoski:2014gra}
Many of these are built with specific physical motivations; for example, certain variables are designed to distinguish between the number of sub-jet ``prongs'', since it is expected that boosted bosons decaying to a pair of quarks would tend to have two subjets, while a top quark may have three.

Some of these observables are designed to also satisfy certain theoretical guarantees, such as Infrared- and Collinear (IRC) safety~\cite{Banfi:2006hf,Dasgupta:2007wa,Moult:2016cvt}.
Such quantities are insensitive to the effects of soft radiation and colliniar splitting, allowing for accurate measurements free of systematics arising from nonperturbative modeling.
These developments in jet substructure have also yielded practical benefits for LHC experimentalists, with the development of improved pileup suppression and jet grooming algorithms~\cite{Soyez:2012hv,Larkoski:2014wba}.

In addition to its role in surpressing QCD backgrounds, jet substructure allows experiments to probe the difference between jets initiated by quarks vs. gluons.
Although formally defining jet flavor at the hadron level is a hazy concept~\cite{Komiske:2018vkc,Banfi:2006hf}, it is certainly the case that individual jets can be described as ``quark-like'' or ``gluon-like'' based on flavor-specific distributions.
Probing the differential cross sections of these jet categories can improve models for parton showers, parton distribution functions, as well as deriving constraints on QCD couplings.~\cite{Komiske:2018vkc,Gallicchio:2011xq,Dasgupta:2007wa,Gallicchio:2012ez,Larkoski:2014pca}.

More recently, there has been a surge of development in the application of modern deep learning techniques to the area of jet substructure and tagging~\cite{Larkoski:2017jix,Guest:2018yhq}.
This is a problem for which deep learning (DL) is well suited, as the input space of observations is naturally very high dimensional, and existing LHC data and simulations provide a very large volume of data from which to learn.

Many different ``off the shelf'' models developed for and by the field of industrial and academic ML research have already adapted to LHC specific data.
These models range from standard dense networks~\cite{Luo:2017ncs}, to 2D image-based convolutional networks (CNNs)~\cite{Cogan:2014oua,deOliveira:2015xxd}, to recursive networks (RNNs)~\cite{Egan:2017ojy,Cheng:2017rdo,Louppe:2017ipp}.

One tricket facet in the application of DL to problems in jet substructure is the variable size of the input space.
Typical jets at the LHC are comprised of tens to hundreds of pseudo-momentum constituents; yet, a network with a fixed architecture must be able to handle differently-sized jets on equal footing in order to model the data distribution.
In the jet image approach, this collection of consitituents is flattened into a fixed binning scheme, obviating the issue at the cost of generating very sparse images with which standard CNNs may struggle.
RNNs, which naturally operate on variable-length sequences; however, most such models require a ``temopral ordering'' which acts as an inductive bias that does not map meaningfully onto the problem at hand.
The remaining approach is to consider the inputs as an unordered set. These so-called ``point cloud'' methods involve pooling operations which are invariant w.r.t. particle ordering.
The simplest of these architectures is based on Deep Sets~\cite{2017arXiv170306114Z}, where each input point is encoded independently to a latent representation.
Another possibility that has been recently explored~\cite{Shlomi2020GraphNN,Moreno:2019neq,Moreno:2019bmu} is to include pairwise relational features between inputs, which can be modeled with Graph Neural Networks.

In this paper, we focus on the former point cloud approach. In particular, we build upon the work of ~\cite{Komiske:2018cqr} by generalizing the EnergyFlow and ParticleFlow network models to incorporate rotational equivariance.
This leads to performance that is competitive with graph network models, while retaining most of the simplicity of the original formulation.
We also demonstrate how preserving this natural spatial symmetry can lead to more interpretable results.







\section{Deep Sets and Flow Nets}

We will now analyze the connection between Deep Sets, Particle Flow and Energy Flow networks.

\subsection{Deep Sets as Projections}
The Deep Sets theorem~\cite{2017arXiv170306114Z} establishes that certain classes of permutation-invariant functions  $u : 2^X \rightarrow \mathcal{Y}$ acting on sets $S$ of inputs may be universally approximated by a function of the form:
\begin{equation}
\label{eq:deepsets}
u(S) \approx f(\{s_1,\dots,s_k\}) = F\left(\sum_{i=1}^{|S|} \Phi(s_i) \right)
\end{equation}
for suitable functions $\Phi : X \rightarrow \mathbb{R}^L$ and $F : \mathbb{R}^L \rightarrow \mathcal{Y}$ are some functions of sufficient complexity (usually represented by neural networks).
The hyperparameter $L$ is the \textit{latent space dimension}, which must be of sufficient size but is not generally known in advance.

Strictly speaking, the universal approximator property is only proved for functions $u$ on sets of fixed-size, i.e. on $X^k$. Moreover, the Deep Sets theorem relies on results for polynomial approximation, and applies only to functions $u$ which are continuous with respect to the Euclidean metric on $X^k\cong \mathbb{R}^{k\times d}$.
Some practical limitations of this as well as weaker conditions for universality in more general settings have been developed in~\cite{NEURIPS2018_f5f8590c}.
Nonetheless, the architecture of Eq.~\ref{eq:deepsets} has yielded competetive results, even when the underlying space $X$ encodes features such as color, size, energy, etc, and consequently has meaningful topology at all.
In our present study, we will ignore the universal approximator property.

To facilitate our study of symmetry, let us suppose the feature space $X$ has measure $\mu$, and the functions $f,g \in H$ are elements of the Hilbert space $\mathcal{H}=L^2(X,\mu)$ of square integrable functions
We may write the inner product:
\begin{equation}
    \left<f, g\right> = \int_X f^* g d\mu \,.
\end{equation}

We also define improper integral of the Dirac delta function on this space in the usual way\footnote{This might require something stricter than a measure space for $X$, possibly a locally compact group?}, 
\begin{equation}
    \int_X \delta(x - y)f(x)d\mu(x) =: f(y) \,.
\end{equation}

Now we propose that rather than considering $S$ to be a set of features $s\in X$, we instead view it as a generalized function $\mathcal{S}_S : X \rightarrow \mathbb{R}$.
We make this association via the functional defined as $\mathcal{S}_{(\cdot)} : 2^X \times X \rightarrow \mathbb{R}$ where
\begin{equation}
    \mathcal{S}_S(x) \mapsto \sum_{i=1}^{|S|}\delta(x - s_i) \,.
\end{equation}
This object (also known as an empirical distribution) can be viewed as a ``function'' representing the spatial density of the observed elements $s\in S$.

Given these definitions, we can evaluate the projection of $\mathcal{S}_S$ onto an arbitrary function $\Phi\in \mathcal{H}$:
\begin{equation}
\label{eq:delta-projection}
    \left<\Phi, \mathcal{S}_S \right> = \int_X \Phi^*(x))\mathcal{S}_S(x)d\mu = \sum_{i=1}^{|S|}\Phi^*(s_i) \,.
\end{equation}

We immediately recognize this result as the inner term of Ref.~\ref{eq:deepsets}, i.e., we may write:
\begin{equation}
    u(S) \approx F\left( \left< \Phi^*, \mathcal{S}_S \right>\right) \,.
\end{equation}
In other words, the Deep Sets networks can be thought of as a projection of the input onto a collection of learned feature maps.
For later convenience, we have generalized $\Phi$ to be an arbitrary complex-valued function.

\subsection{Energy Flow Networks}
In Ref.~\cite{Komiske:2018cqr}, the Particle Flow Network (PFN) is defined to be exactly as in Eq.~\ref{eq:deepsets}, with the input space $x=(p^\mu, q_0,\dots,q_n) \in X$, where $p^\mu$ is a particle's 4-momentum, and $q_i$ are arbitrary additional per-particle features such as charge or species.
The Energy Flow Network (EFN) is defined as a special case, where $\Phi$ is restricted such that
\begin{equation}
\label{eq:efn-linear}
\Phi(p^\mu,q_0,\dots,q_n) = z\  \Phi_{efn}(\hat{p}^\mu) \,,
\end{equation}
that is, $\Phi$ depends only on the momentum direction $\hat{p}^\mu$ and an energy scale\footnote{Both the particle energy $E$ and transverse momentum $p_\mathrm{T}$ are reasonable choices for the scalar energy; we shall assume $z=p_\mathrm{T}$, as it is invariant under translation in $\phi$ and $\eta$.} $z$.
$\Phi$ is also held to be linear in $z$; this is the basis for the EFN's IRC-safety.

Although arbitrary spatial directions have the topology of a sphere, $\hat{p}^\mu \in S_2$, momenta represented in detector coordinates $\vec{x} = (\eta, \phi)$ locally resemble $\mathbb{R}^2$ on scales $|\delta \vec{x}| \ll \pi$.
In order to better utilize this simple topology, we redefine our feature space to be $W = X \times V$, where $X=\mathbb{R}^2$ represents local spatial direction, and $V=\mathbb{R^+}$ represents energy.

This suggests defining a new functional, $\mathcal{E}_{(\cdot)} : 2^W \times X \rightarrow V$, which maps a jet's constituents $S = \{(\vec{x}_i, z_i)\}$ to a function defined by:
\begin{equation}
    \mathcal{E}_S(\vec{x}) = \sum_{i=1}^{|S|} z_i \delta(\vec{x}-\vec{x}_i) \,.
\end{equation}
Since $\vec{x}$ now has a meaningful topology, we can easily interpret this quantity as the directional distribution of energy for the observed jet $S$, i.e. this quantity is the jet's \textit{empirical energy flow}.
We can now replace the summation term in the EFN with the projection of this energy flow onto arbitrary functions $\Phi$ in the well-understood space of functions that are square-integrable on $\mathbb{R}^2$:
\begin{equation}
    u_{efn}(S) = F\left(\sum_{i=1}^{|S|} z_i \Phi^*(\vec{x}_i) \right) = F\left( \left< \Phi, \mathcal{E}_S \right> \right) \,.
\end{equation}

We emphasize the point that while the appearance of Eq.~\ref{eq:efn-linear} suggests that the EFN is a restriction of the functional form of $\Phi$, an alternative interpretation is that the EFN is a \textit{generalization} of a deep sets-like projection on the topological space $X$.
It is a generalization in the sense that we have promoted the empirical distribution of observed points $\sum \delta(\vec{x} - \vec{x_i})$ to an energy-valued scalar field $\mathcal{E}_S(\vec{x})$ over $X$.
It is a scalar field in that the values of the field at point $\vec{x}$ are invariant under group transformations on $X$.
However, the field may take values in a vector space such as $\mathbb{R}^n$, as we shall see in the next section.

For example, the triplet of RGB color channels in a 2D image is a scalar under rotations; rotating the coordinate system will reorient the image, but does not affect the local color of any pixels.
There are interesting scenarios in which the particle cloud have non-scalar features, for example position and momentum.
In this case, a rotation of the position space will also locally affect the orientation of momentum vectors.
The representational considerations for neural networks operating on such data has been explored recently~\cite{2018arXiv180208219T,pmlr-v119-bogatskiy20a}, but in the present work we focus on the scalar-field case.

\subsection{Particle Flow Networks}
Adopting this viewpoint of the particle point-cloud as a scalar field, we propose an alternative form for the PFN; this definition again involves defining a functional $\mathcal{J}_S : 2^W \times X \rightarrow V'$
\begin{equation}
\label{eq:zmap}
    \mathcal{J}_S(\vec{x}) \mapsto \sum_{i=1}^{|S|} Z(\vec{q}_i) \delta(\vec{x} - \vec{x}_i) \,,
\end{equation}
where $V$ now represents the space of any non-directional features (energy, charge, mass, etc.), and $Z : V \rightarrow V'$ is an arbitrary (nonlinear) function mapping from those features to a latent space $V' \subseteq \mathbb{C}^L$.
Note that the learnable filter(s) $\Phi$ remain an arbitrary function on $X$, as this is the relevant space over which the inner product is defined.

The map $Z$ allows for the transformation of arbitrary non-spatial features so that
$\mathcal{J}_S$ represents a scalar field taking values in the $L$-dimensional latent space $V'$.
This map could be a statically chosen function, or itself a trainable neural network.
The case $Z(\vec{q})=z$ corresponds to the EFN architecture.
With this representation, our network is defined as
\begin{equation}
    u_{pfn}(S) = F\left(\sum_{i=1}^{|S|} \Phi_{pfn}(\vec{x}_i, \vec{q}_i) \right) = F\left(\left< \Phi, \mathcal{J}_S \right> \right)
\end{equation}

Ultimately, adopting the scalar field representation of Eq.~\ref{eq:zmap} corresponds to simply requiring the function $\Phi_{pfn}(\hat{p}^\mu, \vec{q}) = \Phi(\hat{p}^\mu)Z(\vec{q})$ be separable in the direction features.
Indeed, preventing the learnable function from observing the joint distribution of spatial and non-spatial features seems a suboptimal restriction.
As we will see, however, there are a number of advantages to this form.

\section{Equivariant Projection Filters}
EFNs and PFNs work by learning 2D pattern maps, and calculating latent features by correlating inputs with these filters.
However, due to continuous factors of variation in the data, it could take infinitely many filters in order to optimally match even a one-dimensional data manifold.
As an example, consider a jet containing two distinct subjets, arising from a quark pair within the jet, illustrated in Fig.~\ref{fig:cartoon}.
Due to an (approximate) physical symmetry about the jet axis, the orientation of these subjets relative to the $+\eta$ axis is essentially randomly distributed from $[0, 2\pi)$.
That is, in most cases we expect jets which differ only by an $SO(2)$ rotation about the axis to be physically identical, and should have equal density in the data distribution.
Therefore, in order for a PFN to respond identically each one of these orientations, it must either 1) learn infinitely many filters corresponding to each rotation, 2) learn limit resolution in such a way that filters are not sensitive to orientation, or 3) respond in unpredictable ways when rotations are applied.

\begin{figure}[t]
    \centering
    \includegraphics[width=0.45\textwidth]{figures/cartoon.pdf}
    \caption{Illustration of the inner product between a collection of observed jet constituents (red dots with brightness indicating $p_\mathrm{T}$) and learned filter maps $\Phi_{1,2}(\eta, \phi)$.
    The observed jets $J_1$ and $J_2$ are completely identical up to a rotation about the jet axis.
    Yet an EFN would require two filters $\Phi_1(\eta, \phi)$ and $\Phi_2(\eta, \phi)$ in order to detect the same feature in both jets.
    }
    \label{fig:cartoon}
\end{figure}

We can ameliorate this issue by imposing some structure on the space of functions allowed for $\Phi$.
One possibility would be to constrain the form of $\Phi$ so that it is azimuthally symmetric about the jet axis.
In this case, a given jet would give the same filter response regardless of its orientation; this is an example of \textit{invariance}.
Although in most applications we do in fact want a network response which is invariant to such a symmetry, this would be a severe restriction of the functions that could be learned by the network.
Instead of forcing the network to immediately produce an invariant response, we can instead preserve the structure of the symmetry by enforcing \textit{equivariance} at each layer.
A particular function $f : X\rightarrow Y$ is equivariant w.r.t. a group $G$ acting on $X$ if there is some representation $\alpha$ of $G$ on $Y$ s.t. for all $g\in G$ and $x \in X$
\begin{equation}
\label{eq:equivariance}
    f(g\cdot x) = \alpha(g) f(x) \,.
\end{equation}
In other words, equivariance of $f$ implies that the group transformation $x\rightarrow g\cdot x$ can be effected either before evaluating the function (by acting $g$ on $\vec{x}$), or afterwards (by acting $\alpha(g)$ on the result).

To obtain an equivariant filter, we could require $\Phi(r, \theta)$ to be periodic in $\theta$, e.g. $\Phi(r, \theta) = \rho(r) e^{i m \theta}$ for integer $m$ and arbitrary radial profile $\rho$.
It is easy to verify the equivariance of this construction, as $\Phi(R_\Delta \vec{x}) = e^{i m \Delta} \Phi(\vec{x})$.
However, this is also too restrictive a form to be useful, as salient features in jet substructure are not necessarily periodic.

Instead, we consider the more general case in which $\Phi$ is a \textit{steerable function}.
A steerable function is one which may be written as a linear combination of equivariant basis functions ${\phi_i}$, so that
\begin{equation}
\label{eq:linear-combination}
    \Phi(g \cdot \vec{x}) = \sum_{i=1}^{n} \alpha_i(g) \phi_i(\vec{x}) = \boldsymbol{\alpha}^\mathrm{T}(g) \mathbf{\Phi}(\vec{x}) \,.
\end{equation}
This expression does not exactly conform to the equivariance condition~\ref{eq:equivariance}.
However we may associate the $n$-dimensional vector space $\mathcal{F} = \mathrm{span}\{\phi_1, \dots, \phi_n\}$ with a subspace of the original Hilbert space, $\mathrm{F} \subseteq \mathcal{H}$.
Then in terms of the column vector $\mathbf{\Phi}$, we have
\begin{equation}
\label{eq:interpolation}
    \mathbf{\Phi}(g \cdot \vec{x}) = \mathbf{A}(g) \mathbf{\Phi}(\vec{x}) \,,
\end{equation}
where $\mathbf{A}(g) = \mathrm{diag}\{\alpha_1(g), \dots, \alpha_n(g)\}$, which provides a manifestly equivariant form.
This is known as the interpolation equation, and depending on the group $G$ and the choice of basis functions $\{\phi_i\}$, $\mathbf{A}$ need not be diagonal.
Hence, this type of equivariant representation generalizes Eq.~\ref{eq:linear-combination}; the function space, rather than the function itself must be closed under transformations.
Function spaces whose elements satisfy Eq.~\ref{eq:interpolation} are called \textit{equivariant function spaces}.

Functions in an equivariant space $\mathcal{F}$ may be much more generic and expressive than any particular equivariant function.
Returning to the application of particle flow networks, we find that in order to learn sufficiently interesting filters, while preserving equivariance, we must represent intermediate results in the vector space $\mathcal{F}$.
For convenience, we define the projection of $\mathcal{J}_S$ onto $\mathbf{\Phi} \in \mathcal{F}$ as
\begin{align}
    \left<\mathbf{\Phi}, \mathcal{J}_S \right> := \begin{pmatrix}
    \left< \phi_1 , \mathcal{J}_S \right> \\
    \vdots \\
    \left< \phi_n , \mathcal{J}_S \right>
    \end{pmatrix} \in \mathbb{C}^n \,.
\end{align}

A rotation $R_\Delta$ of the jet $\mathcal{J}_S$ is equivalent to a change of coordinates so that $R_\Delta \mathcal{J}_S(\vec{x}) = \mathcal{J}_S(R^{-1}_\Delta \vec{x})$.
From the definition\footnote{fix this argument} of the inner product on $\mathcal{H}$,
\begin{align}
    \left<\phi_i, R_\Delta \mathcal{J}_S \right> &= \int_X \phi_i(\vec{x}) \mathcal{J}_S(R_\Delta^{-1} \vec{x}) \\
    &= \int_X \phi_i(R_\Delta \vec{x})\mathcal{J}_S(\vec{x}) \\
    &= \left< R^{-1}_\Delta \phi_i, \mathcal{J}_S \right>
\end{align}
From the interpolation Eq.~\ref{eq:interpolation},
\begin{align}
    \left<R_\Delta^{-1} \phi_j, \mathcal{J}_S \right> = A^*_{ij}(\Delta)\left<\phi_i, \mathcal{J}_S\right> \,,
\end{align}
and therefore the equivariance of the projection onto $\mathbf{\Phi}$ is established in terms of the interpolation matrix $\mathbf{A}$:
\begin{align}
    \left< \mathbf{\Phi}, \mathcal{J}_S \right> = \mathbf{A}^*(\Delta)\left< \mathbf{\Phi}, \mathcal{J}_S \right>
\end{align}

\section{SO(2) Equivariance}
In the preceeding section, we established a general formulation for designing projective filters by restricting the learnable filters to an equivariant function space and utilizing the interpolation matrix.
We now consider a concrete example for the rotation group $SO(2)$.
It can be shown~\cite{cheng1998theory} that the most general $SO(2)$-equivariant space of functions $\phi : \mathbb{R}^2 \rightarrow \mathbb{C}$ are spanned by functions of the form $\rho(r)e^{i m \theta}$, where $r$ and $\theta$ are polar coordinates on $\mathbb{R}^2$, in our case corresponding to $r=\sqrt{\eta^2 + \phi^2}$ and $\theta = \tan^{-1}(\phi/\eta)$.
Therefore the simplest possibility is to define our filter as $\mathbf{\Phi} = (\phi_{M}, \phi_{M-1}, \dots, \phi_{-M})^\mathrm{T}$ where
\begin{equation}
    \phi_m(r, \theta) = \rho_m(r) e^{i m \theta} \,.
\end{equation}
Here $M$ is a hyperparameter that represents an angular frequency cutoff for features.
Since the functions $\rho_m : \mathbb{R}^{\geq 0} \rightarrow \mathbb{C}$ may be of arbitrary form, we can parameterize them as trainable neural networks.
Because of the orthogonality of the complex exponentials, the interpolation matrix is diagonal and given simply by
\begin{equation}
\label{eq:interpolation-so2}
\mathbf{A}(\Delta) = \mathrm{diag}\{e^{i M \Delta}, \dots, e^{-i M \Delta}\} \,.
\end{equation}
That is, under rotation $R_\Delta$, each of the components of $\mathbf{\Phi}$ (and hence $\left \mathbf{\Phi}, \mathcal{J}_S \right>$) are simply modulated by a phase with frequency $m$.

\subsection{Convolution}
Now that we have obtained a trainable equivariant intermediate representation, $\vec{y} = \left<\mathbf{\Phi}, \mathcal{J}_S\right> \in \mathbb{C}^{2M+1}$, we will generally want to proceed to deeper representations by passing to additional network layers.
In order to preserve equivariance, the subsequent layer $f : \mathbb{C}^{2M+1} \rightarrow X'$ will need to satisfy the equivariance condition
\begin{equation}
    f(\mathbf{A}(\Delta) \vec{y}) = \alpha(\Delta) f(\vec{y})
\end{equation}
for some representation $\alpha$ of $SO(2)$ on $X'$.

The most general group-equivariant operation is that of generalized convolution (or correlation)\cite{}, which is defined for locally compact groups.
Using the interpolation equation, we may consider the projection of an arbitrarily rotated jet to be a function $g : SO(2)\rightarrow\mathbb{C}^m$:
\begin{equation}
    g(\Delta) = \left<\mathbf{\Phi}, R^{-1}_\Delta \mathcal{J} \right> = \mathbf{A}(\Delta)\left<\mathbf{\Phi}, \mathcal{J}_S \right> \,.
\end{equation}
The generalized convolution for signals on $SO(2)\cong U(1)$ looks essentially the same as usual 1D convolution, except that all functions on $SO(2)$ have cyclic boundary conditions:
\begin{align}
    [f \star g]_{U(1)}(\Delta) &= \int_{-\pi}^{\pi} d\theta f(\theta)g(\Delta - \theta) \\
    &= \left(\int_{-\pi}^{\pi} d\theta f(\theta) \mathbf{A}(\Delta - \phi)\right)\left<\mathbf{\Phi}, \mathcal{J}_S\right> \,.
\end{align}
We have already obtained the analytic form of $\mathbf{A}$; however, depending on the convolution kernel $f$, we may or may not be able to effect this integral.
In practice, convolution networks use a discretized approximation of the convolution, so that we may simply sample $f$ and $\mathbf{A}$ at intervals:
\begin{equation}
    [f \star g]_{U(1)}(\Delta_j) = \sum_i f_i \mathbf{A}(\Delta_j - \theta_i)
\end{equation}
where $f_i = f(\theta_i)$ are the learnable parameters of the covolutional filter kernel.

{\color{blue} editorial? should i drop the convolutional story or should I add in the results for the convolutional networks?}

\subsection{Invariance}
General equivariant representations will be more powerful than invariant ones.
Once an invariant representation is formed, all subsequent representations computed by the network must necessarily be invariant as well.
But since rotational invariance is desirable for most jet tagging objectives, an obvious (and simple) choice is to compute an invariant at some point before the output of the network.

In particular, the unitarity of $\mathbf{A}(\Delta)$ in Eq.~\ref{eq:interpolation-so2} ensures that the squared modulus $|\left< \mathbf{\Phi}, \mathcal{J}_S \right>|^2 := \sum_m |\left<\phi_m | \mathcal{J}_S \right>|^2$ is invariant under rotation.
It also has the benefit of being real-valued, making it easier to adapt to conventional neural network layers.

Using Parseval's theorem, it is possible to rewrite this as:
\begin{equation}
\sum_m |\left<\phi_m, \mathcal{J}_S \right>|^2 = \frac{1}{2 \pi}\int_{-\pi}^{\pi} |\left< \Phi , R_\Delta^{-1}\mathcal{J}_S \right>|^2 \,,
\end{equation}
Therefore we can interpret this invariant as the \textit{continuous average} projection of $\mathcal{J}_S$ onto the filter, averaged over all possible rotations.
In our experiments, we found that simply computing this norm was nearly as effective as including additional convolutional layers.


\begin{equation}
    \left<\mathbf{\Phi}, \mathcal{J}_S\right> \overset{\vec{x}\rightarrow g\cdot\vec{x}}{\longrightarrow}
    \left<\mathbf{\Phi}, g^{-1} \mathcal{J}_S \right> = \mathbf{A}(g)\left< \mathbf{\Phi}, \mathcal{J}_S \right> \,.
\end{equation}

Bra-ket notation lends itself well to this dualistic interpretation; $\left| \Phi \right> = \sum_i a_i \left|\phi_i\right>$ represents a column-vector in any basis $\{\left|\phi_i\right>\}$ of $\mathcal{F}$, while we understand $\left<\vec{x} | \Phi \right> = \sum_i \left<\vec{x} | \phi_i\right> = \sum_i \phi_i(\vec{x})$.
Then in terms of elements of the function space $\mathcal{F}$, we recover the equivariance condition
\begin{align}
    \Phi(g \cdot \vec{x}) &= \sum_i \alpha_i(g)\left< \vec{x} | \phi_i \right> = \left<\vec{x} | \mathbf{A}(g) | \Phi \right> \\
    \implies& \left|g^{-1} \Phi \right> = \mathbf{A}(g) \left| \Phi \right>
\end{align}
where $\mathbf{A}(g) = \mathrm{diag}\{\alpha_1(g), \dots, \alpha_n(g)\}$.
This is known as the interpolation equation, and depending on the group $G$ and the choice of functions $\{ \phi_i \}$, $\mathbf{A}$ does not necessarily need to be diagonal.
Hence, this type of equivariant representation generalizes Eq.~\ref{eq:linear-combination}, and function spaces whose elements satisfy Eq.~\ref{eq:interpolation} are called \textit{equivariant function spaces}.

Returning to the application of flow projections, we suppose that $\Phi$ is restricted to some equivariant function space $\Phi \in \mathcal{F}=\mathrm{span}\{\phi_1,\dots,\phi_n\}$.
Using the completeness of $\{\phi_i\}$ on $\mathcal{F}$, we may then write
\begin{align}
    \left< \Phi | \mathcal{J}_S \right> &= \sum_i \left< \Phi | \phi_i \right>\left<\phi_i | \mathcal{J}_S \right> \\
    &= \sum_i \left< \Phi | \phi_i \right> \sum_{(\vec{q},\vec{x}) \in S} Z(\vec{q}_j) \phi_i^*(\vec{x}_j) \\
    &= \sum_i c_S \left< \Phi | \phi_i \right> \,,
\end{align}
where $c_S \in \mathbb{C}$.
Effecting a rotation of the jet $\mathcal{J}_S$ by an angle $\Delta$, we have
\begin{equation}
    \left<\Phi|R_\Delta \mathcal{J}_S \right> = \sum_i c_S \left< \Phi | A(\Delta) | \phi_i \right> Phi_i A_ik phi_kj
\end{equation}


we see that when $\mathbf{\Phi}$ is constructed as above, we should rewrite
\begin{align}
    \left< \Phi | R_\Delta \mathcal{J}_S \right> &= \left<R_\Delta^{-1} \Phi | \mathcal{J}_S \right> \\
    &= \sum_i \left<R_\Delta^{-1}\Phi | \phi_i \right>\left< \phi_i | \mathcal{J}_S \right> \\
    &= \left< \Phi | \mathbf{A}(g^{-1}) | \phi_i \right>\left<\phi_i | \mathcal{J}_S \right>
\end{align}

To this end, we propose expressing $\Phi$ in terms of some basis of filters, each of which is equivariant.
$\Phi(\vec{x}) = \sum_m a_m \phi_m(\vec{x})$

\begin{align}
    \left< \Phi, R_\Delta \mathcal{J}_S \right> &= \int_X \Phi^*(\vec{x}) \mathcal{J}_S(R^{-1}_\Delta \vec{x}_i) \\
    &= \int_X \Phi^*(R_\Delta\vec{x}) \mathcal{J}_S(\vec{x}) \\
    &= \int_X \alpha(\Delta)\Phi^*(\vec{x})\mathcal{J}_S(\vec{x}) \\
    &= \alpha(\Delta) \left< \Phi, \mathcal{J}_S \right>
\end{align}

In order to endow our network with the desired structure, we will require the filters $\Phi$ to be an element of an \textit{equivariant function space} $\mathcal{F} = \mathrm{span}\{\phi_1, \phi_n\}$.
An $n$-dimensional equivariant function space has the property that for every $g\inG$, there is an $n \times n$ matrix $\mathbf{A}(g)$ such that
\begin{equation}
\label{eq:interpolation-old}
    \mathbf{\Phi}(g \cdot \vec{x}) = \mathbf{A}(g) \mathbf{\Phi}(\vec{x}) \,,
\end{equation}
where the bold-face notation $\mathbf{\Phi}$ denotes that we express the function $\Phi(\vec{x}) = \sum c_i \phi_i(\vec{x})$ as an $n$-dimensional column vector $(c_1, \dots, c_n)^\mathrm{T} \in \mathcal{F}$.
Eq.~\ref{eq:interpolation} is known as the \textit{interpolation equation}.

In the case of one-parameter transformation group, the interpolation matrix may be expressed as $A(g) = e^{g B}$ for an arbitrary matrix $B$.

For the specific case of Particle Flow networks, we shall consider equivariance w.r.t. the group of rotations about the jet axis, $SO(2)\cong U(1)$.
Since $SO(2)$ is a one-parameter group transformation, it is always possible to change to coordinates in which the group action resembles translation, with infinitesimal generator $L_t = \partial/\partial x$.

In this case the relevant coordinates are polar coordinates in the $(\eta, \phi)$ plane, defined by $r=\sqrt{\eta^2 + \phi^2}$, $\theta = \tan^{-1}(\phi/\eta))$.
In these coordinates, a rotation of $\Delta$ about the jet axis is given by $f(R_\Delta \vec{x}) = f(r, \theta + \Delta)$.

For any Lie group, it is possible to construct a space of equivariant functions satisfying Eq.~\ref{eq:equivariance} by requiring closure under the exponential map of the group's generators~\cite{cheng1998theory}.
In the case of $SO(2)$, the generator of rotation is $L = -y \partial/\partial x + x \partial/\partial y$, and any equivariant function space is of the form $h(r)$

for any such group; for $SO(2)$ any such space is spanned by functions of the form $f(r, \theta) = e^{i m \theta} R(r)$ for $m \in \mathbb{Z}$.
Therefore, we will take $\Phi$ to be of the form:
\begin{equation}
    \Phi(\vec{x}) = \Phi(r, \theta) = \sum_{m=-M}^{M} R_m(r) e^{i m \theta} =: \sum_{m=-M}^{M} \phi_m(r, \theta) \,,
\end{equation}
where $M$ is a hyperparameter associated with some maximum azimuthal frequency, and $R_m : \mathbb{R}^{\geq 0} \rightarrow \mathbb{C}$ are arbitrary functions.

Clearly, for any rotation $R_\Delta$, we have $\phi_m(R_\Delta\vec{x}) = \phi_m(r, \theta + \Delta) = e^{i m \Delta} \phi_m(r, \theta)$.
If we view $\Phi$ as a column vector with $2M+1$ components, i.e.:
\begin{equation}
    \vec{\Phi}(\vec{x}) = \left( \phi_{-M}(\vec{x}), \phi_{-(M-1)}(\vec{x}), \dots, \phi_M(\vec{x}))
\end{equation}

\section{Deep Sets as Projections}
\label{sec:deepsets}

The Deep Sets architecture~\cite{2017arXiv170306114Z} provides a permutation-invariant universal approximator for continuous functions $f\colon 2^X\rightarrow \mathcal{Y}$ which act on sets of features $S\in2^X$ where $S = \{s_i\}_{i\in \mathrm{Ind}_S}; s_i\in X$.
This is the case when $f$ is of the form
\begin{equation}
    f(S) = f(\{s_1,\dots,s_k\}) = F\left(\sum_{i=1}^{|S|} \Phi(s_i) \right) \,,
\end{equation}
for suitable latent space dimensionality $L$ and functions $\Phi\colon X\rightarrow\mathbb{R}^L$ and $F\colon\mathbb{R}^L\rightarrow\mathcal{Y}$.
In particular, when $F$ and $\Phi$ are modeled by neural networks, $f$ is a neural network which can be trained with the usual algorithms, and is manifestly invariant with respect to the input order of the features $x_i$.

Strictly speaking, the universal approximator property is only proved for functions $f$ on sets of fixed-size, i.e. on $X^k$.
The Deep Sets theorem relies on results for polynomial approximation, and applies only to functions $f$ which are continuous with respect to the Euclidean metric on $X^k\cong\mathbb{R}^{k\times d}$.
Weaker conditions have been obtained for universality, e.g.~\cite{NEURIPS2018_f5f8590c}.
Nevertheless, the architecture has seen successful application even when the underlying feature space $X$ represents quantities such as color, words, etc., which does not have a meaningful topology.

When the feature space is measurable\footnote{We suppress the specification of $\sigma$-algebras throughout the discussion, assuming only that $(X,\Sigma,\mu)$ forms a measure space for \textit{some} $\Sigma$.} with measure $\mu$, we may naturally consider the inner product between functions $f,g$ in the Hilbert space $H=L^2(X,\mu)$:
\begin{equation}
    \langle f, g \rangle = \int_{X} f^* g d\mu \,.
\end{equation}


We shall also define the improper integral of the Dirac delta function on this space as {\color{blue}error: this requires binary operation $+$ with inverses, so we probably require $X$ to be a locally-compact group}
\begin{equation}
\label{eq:deepsets-proj}
    \int_X \delta(x-y)f(x)d\mu(x) := f(y) \,.
\end{equation}

Now we propose that rather than considering $S$ to be a set of features $s\in X$, we instead view it as a generalized function $\mathcal{J}_S\colon X\rightarrow \mathbb{C}$.
We make this association via the functional defined as ${\mathcal{J}_{(\cdot)}\colon 2^X\times X}\rightarrow\mathbb{C}$ where
\begin{equation}
    \mathcal{J}_S(x) \mapsto \sum_{i=1}^{|S|} \delta(x - s_i) \,.
\end{equation}
This object can be viewed as a ``function'' representing the spatial density of the observed elements $s\in S$.

Given these definitions, we can evaluate the projection of the object $\mathcal{J}_S$ onto the function $\Phi\in H$:
\begin{equation}
\label{eq:projection}
    \langle \Phi, \mathcal{J}_S \rangle = \int_X \Phi^*(x) \mathcal{J}_S(x) d\mu(x) = \sum_{i=1}^{|S|} \Phi^*(s) \,.
\end{equation}
We immediately recognize this quantity as the inner term of \cref{eq:deepsets}, for the case of real-valued $\Phi$.
Therefore, when acting on spatial (or otherwise locally compact) coordinates $X$, the Deep Sets architecture may be understood as a function $F$ applied to the projection of the spatial distribution of observations onto learned filters $\Phi\in H$:
\begin{equation}
    f(S) = F\left(\langle \Phi, \mathcal{J}_S \rangle \right) \,.
\end{equation}


\section{Energy Flow Networks}
\label{sec:energyflow}

Energy Flow Networks~\cite{Komiske:2018cqr} are a special case of the Deep Sets architecture, in which coordinate space $X\subset\mathbb{R}^3$ represents the momentum 4-vectors of massless particles measured within an energetic hadronic jet.
Each particle can represented by the coordinates $(z,\eta,\phi)$, where $z$ is an energy scale (taken to be the transverse momentum $p_\mathrm{T}$) and $(\eta,\phi)\in\mathbb{R}^2$ are the angular coordinates of the particle relative to the jet axis.
Specifically, the EFN corresponds to the case where $\Phi(z,\eta,\phi)$ is linear in $z$, by construction:
\begin{equation}
\label{eq:energyflow}
    f_\mathrm{efn} = F\left( \sum_{i=1}^{|S|} z_i \Phi(\eta_i, \phi_i) \right) \,.
\end{equation}
In the notation of \cref{eq:deepsets-proj}, we may write
\begin{equation}
    f_\mathrm{efn} = F\left( \langle \Phi, \mathcal{E}_S \rangle \right)
\end{equation}
where the functional $\mathcal{E}_{(\cdot)}\colon \mathbb{R}^2\rightarrow\mathbb{C}$ defined by
\begin{equation}
    \mathcal{E}_S(x) \mapsto \sum_{i=1}^{|S|} z_i \delta(\eta - \eta_i)\delta(\phi - \phi_i)
\end{equation}
has a natural physical interpretation as the spatial distribution of energy (or momentum) within jet.

While a network of the form in \cref{eq:energyflow} clearly has much less expressive power than \cref{eq:deepsets}, this particular case is of special interest in high energy physics applications.
It is shown in~\cite{Komiske:2018cqr} that for suitable $F$ and $\Phi$, \cref{eq:energyflow} is a universal approximator for a specific class of functions on the particles of a jet, known as IRC-safe observables.

Essentially, IRC-safety means that the observable is not measurably different for any two events which differ only by the collinear splitting of an energetic particle, or the infrared emission of additional low-energy particles.
IRC-safety is an important invariant which applies to many perturbative calculations in QCD theory, allowing for these theoretical values to be compared to experimental measurements of non-perturbative phenomena.
As experimentalists rely increasingly on deep learning methods in analyses, the EFN is an important step towards learning powerful models with appropriately regularized representations.

\section{Dual Space Regression}
\label{sec:tensordual}

In \cref{sec:energyflow}, we see an example of how linearity of the mapping $\Phi$ with respect to a subspace of the features $X$ can have useful properties.
We will generalize this concept further by supposing that the input feature space is $W=X\times V$, where $X$ is required to have locally-compact topology as described in \cref{sec:deepsets}.
For the present purposes, we assume only that ${V\subseteq\mathbb{C}^K}$ for some dimension $K$.
We say that each element of a set $S\in2^W = \{(x_i,v_i)\}$ has \textit{coordinate features} $x_i\in X$ and \textit{qualitative features} $v_i\in V$.

Then we may define a new functional $\mathcal{G}_{(\cdot)}\colon 2^W\times X\rightarrow V$ as
\begin{equation}
\label{eq:functional}
    \mathcal{G}_S(x) \mapsto \sum_{(x_i,v_i)\in S} v_i \delta(x-x_i) \,,
\end{equation}
In particular, the EFN corresponds to the case where the $v_i\in\mathbb{R}$ are particle energies.

Now we may write the projection operation between a filter $\Phi$ and functional $\mathcal{G}$:
\begin{equation}
\label{eq:proj-dual}
    \langle \Phi, \mathcal{G}_S \rangle = \int_X\Phi^* \mathcal{G}_S d\mu = \sum_{(x_i,v_i)\in S} \Phi^*(x_i) v_i \,.
\end{equation}
Here we take the map $\Phi\colon X\rightarrow V$ to be an element of some Hilbert space.
Then by the Riesz representation theorem, $\langle \Phi, \cdot\rangle$ uniquely represents an element of the dual space $\mathcal{H}^*$.

Strictly speaking, neither $S\in 2^W$ nor the generalized function $\mathcal{G}_S$ are members of any Hilbert space; however, they do have well-defined inner product.
Therefore, by analogy with the Riesz representation theorem, we propose that the map $\langle \Phi, \mathcal{G}_{(\cdot}) \rangle \colon 2^W\rightarrow \mathcal{C}$ may be thought of as an element of the ``dual space'' to the set of inputs $2^W$.

We now consider the case where the qualitative feature space $V$ has some kind of tensor structure associated with the coordinates $X$.
That is, suppose that $G$ is a group that acts {\color{blue}(isometrically?)} on $X$ with group action $\cdot:G\times X\rightarrow X$ and that $G$ has a representation $L\colon G\rightarrow GL(V)$ on $V$.
Then if each element of the input space $2^W$ transforms with group action $(x_i, v_i) \mapsto (g^{-1} x_i, L_g v_i)$, we say that $V$ has a tensor structure (specified by the representation $L$) with respect to the coordinates $X$.
In the notation of the functional in \cref{eq:functional}, a transformation $g\in G$ acts on a tensor-valued $\mathcal{G}_S$ via:
\begin{equation}
    \mathcal{G}_{S}(x) \overset{g}{\rightarrow} L_g \mathcal{G}_S(g^{-1} x)
\end{equation}

The simplest case is when $V$ is a \textit{scalar}: then $L_g$ is the identity on $V$ for all $g\in G$.
So for example, if $x_i$ are 2D coordinates, and $V$ represents say, the colors, or masses of some particle at that location, we expect that under a global coordinate rotation, these features will remain the same even though their location changes, illustrated in [figure].

Another interesting case is when $V$ transforms as a \textit{vector} under $G$.
Taking again the example of 2D rotations, we consider $v_i$ represents the 2D (vector) velocity of a particle at point $x_i$.
Then under a rotation, the 

When the function $\Phi$ is parameterized, trainable neural network layers can be constructed which act on inputs in the form of \cref{eq:proj-dual}.
Such networks can be thought of as learning elements of $\mathcal{H}^*$ which optimally project from potentially high-dimensional input spaces to a single number in $\mathbb{C}$.
These layers can therefore be used as a form of ``pooling'' which is naturally invariant with respect to the ordering of the input set.
Moreover, they provide an interpretable artifact, $\Phi$, which can be thought of as a ``density'' map of qualitative features $v$ as a function of coordinate features $x$.

\section{other}

Given this projective interpretation, we may generalize a bit further by considering a network of the form \cref{eq:deepsets-proj} to be a map $F\colon \mathbb{C}^L\rightarrow\mathcal{Y}$ composed with a learnable collection of dual space elements.
That is, the left-bilinear action $\langle \Phi, \cdot\rangle$ may be considered an element of the dual space of $2^X$.

This interpretation can be generalized further to the case where the feature space is a product of spaces $X=W\times Z$, where $W$ is locally compact and $Z$ is invariant with respect to transformations in $W$.
For example $W$ could represent spatial coordinates of objects in $\mathbb{R}^3$, while $Z$ represents arbitrary features of those objects, such as color, size, etc.
In this case, we write $\mathbf{x} = \{(w_i, z_i)\}$, and define $\mathcal{J}_\mathbf{x}\colon W\rightarrow V$:
\begin{equation}
\label{eq:Jmap}
    \mathcal{J}_\mathbf{x}(x) := \sum_{(x_i,z_i)\in\mathbf{x}} g(z_i) \delta(x - x_i) \,,
\end{equation}
where $g\colon Z\rightarrow V$ is an arbitrary function mapping the non-spatial coordinates to some latent representation space $V$.
Now $\mathcal{J}_\mathbf{x}$ maps from the spatial coordinates of $W$ to a latent representation in $V$.
We also promote the functions $\Phi$ so that $\Phi*\colon W\rightarrow V^*$ takes values in the dual space to $V$.
In this work we shall always consider $V\subseteq\mathbb{C}^K$, so that conjugation $\Phi\rightarrow \Phi^*$ corresponds to the Hermitian transpose.
However, it may be interesting to consider other cases, such as tensor-valued maps.

Taking the definition of \cref{eq:Jmap}, the projections now read
\begin{equation}
    \langle \Phi, \mathcal{J}_\mathbf{x} \rangle = \sum_{(w_i,z_i)\in\mathbf{x}} g(z_i)\Phi^*(w_i) \,,
\end{equation}
and we may define a new network architecture analogous to \cref{eq:deepsets}:
\begin{equation}
\label{eq:factorsets}
    f_\mathrm{proj}(\mathbf{x}) = F\left( \langle \Phi, \mathcal{J}_\mathbf{x} \rangle\right) = F\left( \sum_{(w_i,z_i)\in\mathbf{x}} g(z_i) \Phi^*(w_i) \right)
\end{equation}

We note that \cref{eq:factorsets} is strictly a special case of \cref{eq:deepsets}; essentially all that has changed is that the function mapping from feature space to latent space has been factorized.
In the following sections, we shall demonstrate that factorizing the covariant and invariant features has distinct benefits in certain applications.

As a corollary to the EFN result, we note that there is a sort of generalized ``IRC-safety'' that applies to networks of the form in~\cref{eq:factorsets}.
Specifically, it is invariant under the exchange:
\begin{equation}
\begin{aligned}
&\{\dots, (w_i,g(z_i)), \dots\} \leftrightarrow \\ 
&\{\dots, (w_i+\delta, a g(z_i)), (w_i+\delta', (1-a)g(z_i)), \dots\} \,,
\end{aligned}
\end{equation}
for sufficiently small $\delta,\delta' \in W$ and $a\in[0,1]$.
In other words, elements may be added or removed from the set $\mathbf{x}$, as long as the total quantity of $g(z_i)$ in a localized area remains unchanged.
This follows in analogy with the EFN result from the fact that $\Phi$ is continuous and $\langle \Phi, \mathcal{J}_\mathbf{x} \rangle$ is linear in $g(z_i)$.

\section{Plancherel Networks}
\label{sec:plancherel}

Let $G$ be a {\color{blue} something?} group which acts on $\mathcal{J}$.
The Plancherel network is of the form:
\begin{equation}
    f_\mathrm{plan}(\mathbf{x}) = F\left( \int_G \left|\langle \Phi, g\cdot\mathcal{J}_\mathbf{x} \rangle\right|^2 d\mu(g) \right) \,,
\end{equation}
where $G$ is {\color{blue} some} group and $\mu$ is the Haar measure on $G$.
In other words, this network computes the square modulus of the projection $\langle\Phi,\mathcal{J}\rangle$ averaged over all possible transformations $g\in G$.

\section{Harmonic Flow Networks}

We start with the special case of the Energy Flow Network.
Both EFNs and PFN propose to solve the issue of jet representation in neural networks by leveraging the universal approximation theorem known as the Deep Sets theorem~\cite{2017arXiv170306114Z}.
A network $f:X\rightarrow\mathcal{Y}$ applied to a jet $\mathcal{J}\in X$ is represented by:
\begin{equation}
    f(\mathcal{J}) = F\left( \sum_i z_i \Phi(\eta_i,\phi_i) \right)
\end{equation}

We note that if we choose to view a jet $J$ as the observed energy density as a function of the spatial coordinates $\eta$ and $\phi$, we may write
\begin{equation}
    \mathcal{J}(\eta, \phi) = \sum_i z_i \delta(\eta - \eta_i)\delta(\phi - \phi_i) \,.
\end{equation}
Then for a filter $\Phi$ which is some arbitrary function in $L^2(\mathbb{R}^2)$, we may compute the projection of the jet onto the filter using the ordinary inner product
\begin{equation}
    \left<\Phi, \mathcal{J}\right> = \int d\eta d\phi\ \Phi^*(\eta, \phi) \mathcal{J}(\eta, \phi) = \sum_i z_i \Phi^*(\eta_i, \phi_i) \,,
\end{equation}
which we may immediately recognize as the inner term of Eq.~\ref{eq:energyflow}, for the case of real-valued $\Phi$.
In other words, we may understand energyflow networks as a neural network $F$ applied to the projection of the jet's spatial energy distribution upon a learned filter $\Phi$.

We will first define the specific case of a Harmonic Energy Flow network in order to demonstrate its advantages, while preserving the IRC-safe and permutation invariant features of the EFN.
We will then generalize the concept to show that Harmonic Flow networks can be used on a variety problems.

We begin by noting that an arbitrary $L^2$ function $g:\mathbb{R}^2 \rightarrow \mathbb{C}$ may be expressed in polar coordinates in the form~\cite{}
\begin{align}
    g(r, \theta) &= \sum_{m=0}^{\infty}\sum_{n=m}^{\infty} A_{nm} \Psi_{nm}(r,\theta) \,;\\
    \Psi_{nm}(r,\theta) &= \sqrt{\frac{2 n!}{(n+m)!}}e^{-r^2/2}r^m L^m_n(r^2) e^{i m \theta} \,,
\end{align}
where $L^m_n$ are the generalized Laguerre polynomials.
This is the so-called Laguerre-Gaussian series, and it may come in handy later.
For now, it is sufficient to note in particular that we may therefore choose to write an arbitrary filter $\Phi$ in the form
\begin{align}
    \Phi(r, \theta) &= \sum_m \Phi_m(r,\theta) \,; \\
    \Phi_m :&= R_m(r) e^{i m \theta} \,,
\end{align}
where $R_m(r)$ depends only on the radius.
This is essentially a multipole expansion of $\Phi$.

Assuming $\Phi$ to take this form, we may now write
\begin{equation}
    \left<\Phi,\mathcal{J}\right> = \sum_m \left<\Phi_m, \mathcal{J}\right> \,.
\end{equation}
While $\left<\Phi, \mathcal{J}\right>$ has no particular equivariance with respect to rotations about the jet axis, the individual $\left<\Phi_m, \mathcal{J}\right>$ do.
Under a rotation $R_\Delta$ about the jet axis,
\begin{equation}
\label{eq:equivariant-projection}
    \left<\Phi_m, \mathcal{J}\right> \rightarrow \left<\Phi_m, R_\Delta \mathcal{J}\right> = e^{i m \Delta} \left<\Phi_m, \mathcal{J}\right> \,.
\end{equation}
So the projections of the jet onto the individual angular frequency modes are linearly equivariant with respect to rotations about the jet axis; in particular, they just pick up a phase.

This immediately suggests replacing the argument to $F$ in Eq.~\ref{eq:energyflow} with a vector of $M+1$ projections onto the definite-frequency modes of the filter $\Phi_m$:
\begin{equation}
\label{eq:energyflow-equivariant}
    f(\mathcal{J}) = F\left(\left[\left<\Phi_0, \mathcal{J}\right>, \dots, \left<\Phi_M, \mathcal{J}\right> \right]^\mathrm{T}\right) \,.
\end{equation}
This forces the overall network to preserve equivariance w.r.t. rotations at least until the first layer of the subnetwork $F$.
We call this the equivariant EFN, and note that this formulation is strictly more general than the EFN of Eq.~\ref{eq:energyflow}, since $F$ is always free to learn to simply sum the individual components of the input vector; in the case of a neural network, this would correspond to a neuron connected to the input with all weights equal to one.
Therefore, Eq.\ref{eq:energyflow-equivariant} possesses the same universal-approximation guarantees as Eq.~\ref{eq:energyflow}, and since each of the $\Phi_m$ are continuous, this form also preserves the IRC-safety.

While we shall show that preserving equivariance to a deeper level of representation within the overall network $f$, it is unlikely that the outer subnetwork $F$ will be able to continue to learn equivariant representations without additional structure.
It turns out that the natural way to ensure this is by promoting $F$ to be a $U(1)$ group-convolutional architecture; we treat this topic in detail in \cite{}[OTHER PAPER].

Since the task at hand (quark vs. gluon tagging) is manifestly rotationally invariant, we shall consider imposing rotational invariance, a special case of equivariance.
This can be easily accomplished by simply taking the magnitude of the projection onto each $m$ mode, i.e. by replacing $\left<\Phi_m,\mathcal{J}\right>$ with $\left|\left<\Phi_m,\mathcal{J}\right>\right|^2$ in Eq.~\ref{eq:energyflow-equivariant}.
However, in so doing, we no longer can guarantee that every observable learned by Eq.~\ref{eq:energyflow} can be learned in this form since.
This is obvious, for instance, if $\left<\Phi,\mathcal{J}\right>$ happens to be negative.

It turns out that simply summing the projections onto the various $m$ modes in quadrature is a very powerful way to ensure equivariance while maintaining the universally expressive and IRC-safe features of the EFN~\ref{eq:energyflow}.
This seems counterintuitive, as it is strictly less expressive than Eq.~\ref{eq:energyflow-equivariant}.
However, using the equivariance property of Eq.~\ref{eq:equivariant-projection}, we note that the projection of a rotated jet, $R_\Delta$ onto the filter $\Phi$ is
\begin{equation}
    \langle\Phi, R_\Delta\mathcal{J}\rangle = \sum_m \langle \Phi_m, \mathcal{J}\rangle e^{i m \Delta} \,,
\end{equation}
where $\langle\Phi_m,\mathcal{J}\rangle$ is simply a complex number.
Therefore, by Parseval's theorem we have the fundamental result:
\begin{equation}
\label{eq:parseval}
    \sum_m\left|\langle \Phi_m, \mathcal{J}\rangle\right|^2 = \frac{1}{2\pi}\int_{-\pi}^\pi\hspace{-0.5em}d\Delta \left|\langle\Phi, R_\Delta \mathcal{J}\rangle\right|^2 \,.
\end{equation}

We therefore propose the formulation:
\begin{align}
    f(\mathcal{J}) &= F\left(\sum_m \left|\left<\Phi_m, \mathcal{J}\right>\right|^2 \right) = F\left(\frac{1}{2\pi}\int_{-\pi}^\pi\hspace{-0.5em}d\Delta\left|\langle\Phi,R_\Delta\mathcal{J}\rangle\right|^2\right) \,; \\
    \Phi_m(r,\theta) &= R_m(r) e^{i m \theta} \,.
\end{align}
which we refer to as the invariant form of the Harmonic Flow network.
It is so called because it is invariant with respect to rotations about the origin.

\section{Experiments}
\label{sec:experiments}

\begin{table*}[ht]
    \centering
    \begin{tabular}{l r c c}
    \vspace{0.25 em}
    & Model & AUC & $1/\epsilon_b |_{\epsilon_s = 50\%}$ \\
    \hline
    \hline
    \multicolumn{1}{l}{\textit{IRC-safe}} & & & \\
    & EFN & 0.8824 & 28.6 \\
    & {\color{blue} s-EFN (Ours)} & 0.8895 & 30.5 \\
    \hline
    \multicolumn{1}{l}{\textit{w/o PID}} & & & \\
    & ResNeXt-50 & 0.8960 & 30.9 \\
    & P-CNN & 0.8915 & 31.0 \\
    & PFN & 0.8911 & 30.8 \\
    & {\color{blue} s-PFN (Ours)} & 0.8958 & 32.4 \\
    & ParticleNet-Lite & 0.8993 & 32.8 \\
    & ParticleNet & 0.9014 & 33.7 \\
    \hline
    \multicolumn{1}{l}{\textit{w/ PID}} & & & \\
    & P-CNN & 0.9002 & 34.7 \\
    & PFN-Ex & 0.9005 & 34.7 \\
    & {\color{blue} s-PFN-Ex (Ours)} & 0.9058 & \textbf{37.6} \\
    & ParticleNet-Lite & 0.9079 & 37.1 \\
    & ParticleNet & 0.9116 & 39.8 \\
    \hline
    \end{tabular}
    \caption{
    Comparison of similar network architectures on the quark-gluon tagging dataset of Ref.~\cite{Komiske:2018cqr}.
    ResNeXt is a deep 2D CNN adapted for jet images in~\cite{Qu:2019gqs}, while P-CNN~\cite{CMS-DP-2017-049} is a 1D convolution operating on ordered lists of particles, also implemented in~\cite{Qu:2019gqs}.
    Our s-EFN model is slightly more sensitive than the EFN, while maintaining IRC-safety.
    The s-PFN model achiaves comparable performance to ParticleNet without adding more inputs, by no longer guaranteeing IRC-safety.
    The s-PFN with experimental PID information performs comparably to ParticleNet.
    }
    \label{tab:my_label}
\end{table*}


\section{Conclusion}
\label{sec:conclusion}

\bibliographystyle{unsrt}  

\section{Introduction}
A common problem in high energy physics is the desire to use state-of-the art machine learning methods to solve real-world problems encountered in physics experiments.
Often, the leading edge of machine learning is driven by a standard suite of problems which are motivated by the practical interests of tech industries.
For example, convolutional neural networks (CNNs) for image processing and Recurrent Neural Networks (RNNs) for natural language processing are some of the most well-studied and developed areas in deep learning.
In an effort to reap the benefits of these advances, physicists have been systematically studying the efficacy of these existing methods when applied to typical problems in our field such triggering, event reconstruction, and object tagging.\cite{Cheng:2017rdo,Egan:2017ojy,Luo:2017ncs,deOliveira:2015xxd,Cogan:2014oua,Guest:2018yhq,CMS-DP-2017-049,Shlomi2020GraphNN}

However, the structure and underlying processes of the data for which these architectures were developed are often fundamentally different from those present in physics experiments.
In order to bridge this gap, physicists have often resorted (with a few exceptions\cite{Louppe:2017ipp,pmlr-v119-bogatskiy20a,Moreno:2019bmu,Komiske:2018cqr}) to remedial measures by re-structuring or re-formatting experimental data from its natural representation to conform with the format required by a specific existing model architecture.

We consider the problem of jet-tagging\cite{Asquith:2018igt,Gallicchio:2011xq,Larkoski:2017jix}, in which the goal is to discriminate between various types of hadronic decays based on the detailed structure of energy and tracking measurements from the detector.
Hadronic jets are the most abundant products of proton collisions at the LHC and are formed when energetic quarks or gluons fragment recursively into lower-energy quarks and gluons until stable particles form, a process known as hadronization.
These particles are then observed by the particle detector.
Experimentally, jets are defined by a specific clustering algorithm\cite{Cacciari:2008gp} which groups together energy deposits localized within a certain angular scale set by a chosen radius parameter.
Jets can be initiated directly from hard QCD processes as well as rare particle decays\cite{Plehn:2009rk,Abdesselam:2010pt,Butterworth:2008iy,Shimmin:2016vlc}.
The exceedingly large QCD cross sections often lead to a scenario where the rate of background jets overwhelms the signal of interest.

In this work, we will consider two common tasks: quark/gluon (q/g) identification\cite{Komiske:2018vkc,Gallicchio:2011xq,Dasgupta:2007wa,Gallicchio:2012ez,Larkoski:2014pca}, and top-quark tagging\cite{Banfi:2006hf,Larkoski:2017jix,Asquith:2018igt}.
In q/g identification, the goal is to determine whether a jet originated from a final-state quark or gluon, which subsequently fragmented into the pattern of radiation observed by the detector.
Top quarks decay before hadronization can occur, and frequently result in three collinear jets, which may overlap significantly in the detector when produced at high momentum.
The goal of top-tagging is to determine whether or not a large-radius jet pattern originated from a top quark.

In data, a jet is most naturally represented by a variable-sized collection of four-vectors representing momentum as measured by calorimeter deposits and/or charged particle trajectories.
These vectors may also be annotated with additional information, such as charge or particle type.
The precision with which these constituents are measured depends on the detector properties.
For example, charged particle trajectories yield very precise directional measurements but do not include neutral particles, while calorimeter deposits have coarser directional precision but better energy precision.

While the field of HEP has used certain machine learning methods, such as Multilevel Perceptrons\cite{Kolanoski:1995zn,Bass:1996ez,Milek:1999xm} and Boosted Decision Trees\cite{Roe:2004na}, for decades, these models are not directly amenable to the data in question.
Instead, physicists have constructed a large number of mathematical and heuristic \textit{jet substructure observables} which quantify various properties on a per-jet level.\cite{Thaler:2010tr,Larkoski:2014pca,Komiske:2017aww,Moult:2016cvt,Larkoski:2014gra,Dasgupta:2013ihk,Larkoski:2014wba}
These observables may then be used in standard multivariate analysis techniques.

Developments in the field of deep learning have led to a reexamination of the problem of jet tagging, as new architectures have emerged which can operate on constituent-level data.
Alternative proposals include casting the particle data into a 2D pixel image in order to apply to a CNNs\cite{deOliveira:2015xxd,Cogan:2014oua}, passing the particle's features one at a time as a sequence into an RNN\cite{Egan:2017ojy,Cheng:2017rdo,Louppe:2017ipp}, or embedding groups of nearby constituents into a graph for use in a Graph Neural Networks\cite{Qu:2019gqs}.
While RNNs in particular have seen successful application in experiments, more recently the field has been moving to permutation-invariant set based networks\cite{2017arXiv170306114Z,Komiske:2018cqr} which are easier to train and achieve comparable performance.

Until now, the question of equivariance properties has not received much attention from the field.
The CNN-based jet image approach, which respects an approximate 2D translational symmetry, was for a long time the only example of an equivariant architecture that had been studied.
While translation equivariance has proved highly useful in image processing tasks, it is unclear whether it is physically meaningful in the context of jet substructure.
Very recently, it has been proposed to consider a more physically-relevant class of equivariance, specifically, the Lorentz Group Network\cite{pmlr-v119-bogatskiy20a}. The LGN architecture operates at the constituent level and is fully equivariant with respect to arbitrary lorentz transformations.
However, the LGN has not yet proved to work as well as existing methods, possibly due to the exceedingly large memory structures required to implement sufficiently complex networks.

In this paper, we turn our attention to the curious gap between the PFN and CNNs.
PFNs, which are the current state-of-the-art in most experimental applications, generally perform as well as any other methods (with the possible exception of the newly-proposed GNNs), but they possess no particular type of equivariance.
On the other hand, the CNN approach which does exhibit equivariance, is generally outperformed by the other methods mentioned, possibly due to the sparsity/discretization or due to the equivariance being of the wrong type.

It turns out there is a specific connection between the PFN and the CNN architectures.
They can be viewed from a common mathematical perspective, illustrated schematically in Fig.~\ref{fig:jetland}.
This mathematical connection, which is elaborated in Sec.~\ref{sec:particle-projection} and \ref{sec:rotational-convolution}, is essentially that with a subtle modification, PFNs can be viewed as the geometrical operation of projection.
This projection operation then can be easily promoted to a convolution, which we call Particle Convolution.
The Particle Convolution allows us to build networks which feature equivariance with respect to a much larger class of symmetry groups, while operating directly on the constituent-level data.
We show that the jet image method is a special case of Particle Convolution in the case of a discrete shift operator and binned coordinates, while the EFN is a special case in which the operator is nullary.

In Sec.~\ref{sec:equivariance} we will discuss the concept of equivariance, and motivate the particular case of rotational equivariance for the problem of jet tagging.
We then consider in Sec.~\ref{sec:particle-projection} a formal connection between the permutation-invariant set-based models and the notion of Hilbert space projection.
In Sec.~\ref{sec:rotational-convolution}, we demonstrate how to promote these projective operations to convolutions possessing equivariance properties by construction, using the particular example of rotation.
In Sec.~\ref{sec:steerable-convolutions}, we consider some of the technical challenges in implementing Particle Convolution, and detail a more efficient solution based on the notion of steerable functions.
In Sec.~\ref{sec:experiments} we provide details for experiments conducted on the two benchmark tasks, q/g tagging and top-tagging, and in Sec.~\ref{sec:conclusion} we present results and conclusions.

\section{Equivariance}
\label{sec:equivariance}

In this section, we begin with a formal definition of the mathematical concept of equivariance.
Then, before proceeding to the technical details of how this definition can be applied, we present an intuitive argument for how equivariance can benefit machine learning models in the specific application of jet tagging.

A map $f : X\rightarrow Y$ is \textit{equivariant} with respect to a group $G$ acting on $X$ if for every $g\in G$, there exists some $\Pi_g : Y \rightarrow Y$ such that:
\begin{equation}
    f(g \cdot x) = \Pi_g f(x)\,,\ \forall x \in X \,.
\end{equation}
In other words, given an equivariant function $f$, it is possible to determine the result of $f(g \cdot x)$ by applying either $g$ to the input, or $\Pi_g$ to the output of the function.
Note that invariance is a special case of equivariance: the function $f$ is said to be invariant when $\Pi_g$ is equal to the identity for every $g$.

The most commonly known equivariant neural architecture is the Convolutional Neural Network (CNN).
These networks are equivariant with respect to discrete translational shifts in one or more dimensions.
In particular, two dimensional CNNs excel at computer vision tasks.
The intuitive reason behind this is that CNNs can learn generic features, such as textures or edges, and is able to match any of its learned features at any location on an image.
If a particular feature is shifted to another location in the image, the CNN's response to that feature will be shifted a corresponding amount.
In contrast, a simple fully-connected neural network would need to re-learn a new instance of each feature 


Many recent developments in the field of machine learning have focused on the analysis and design of equivariance properties of neural networks.
In many cases, equivariance can be generalized to aribrary homogeneous spaces via group convolution, and these architectures often lead to improved performance when a relevant symmetry of the data can be exploited\cite{cohen2019general,Weiler2018LearningSF,2018arXiv180208219T}.

In the context of jet tagging, we propose to consider the specific case of rotation.
The \textit{rotational} Particle Convolution Network (rPCN), detailed in Sec.~\ref{sec:rotational-convolution}, evaluates a different kind of convolution with respect to rotation about the jet axis in the $\eta - \phi$ plane.
This case is physically relevant and also happens to be mathematically simple.
The rationale for this particular equivariance stems from the physical processes governing jet formation, which are approximately invariant under rotation about the jet axis.
Therefore, common features may emerge in radiation patterns which differ only in arbitrary rotation about the jet axis.

To address this, some works \cite{deOliveira:2015xxd} have proposed removing the rotational degree of freedom via a pre-processing step that imposes a standardized reference orientation for all jets.
This technique can in some cases improve performance of non-equivariant models, although in other cases it can have a detrimental effect.
In any case, it effectively manages to render the entire model invariant under \textit{global} rotations.

In many applications we do in fact desire an invariant network response; for example, when determining whether a jet originated from a boosted $Z$ boson decay, we expect the same answer regardless of the random orientation of the parent particle's decay axis.
However, rather than constructing an invariant input to eliminate the rotational variation in the data, we can instead attempt to preserve the structure of the underlying symmetry within the network itself by enforcing equivariance at each layer.
This allows a deep network to build rich representations to examine and compare the angular structure of various features.

In the following sections, we will describe one way in which a rotationally-equivariant network can be constructed.
This architecture, which we call a Particle Convolution Network, can be generalized to effect equivariance with respect to additional types symmetry as well.


\begin{figure}[ht]
    \centering
    \includegraphics[width=0.5\textwidth]{figures/clocks.pdf}
    \caption{Illustration of invariant vs. equivariant feature matching.
    Colors correspond different learned filter channels.
    In this example, the input stimuli are a superposition of two radiation patterns, each of which has an arbitrary angular rotation.
    The invariant approach is able to match inputs differing by a global rotation; however, it must learn different features for each relative orientation.
    A network with equivariant filters can match each pattern regardless of relative orientation.
    }
    \label{fig:rotation}
\end{figure}

Figure~\ref{fig:rotation} illustrates this point.
In this example, we consider two independent features that might appear in a jet: a diffuse blob of low-energy particles, and a pair of high-$p_\mathrm{T}$ subjets.
If the diffuse particles and the hard subjets are produced at different stages during parton showering, their relative angular orientations about the jet axis could be largely random and uncorrelated.
Rotational pre-processing does effectively limit the amount the network needs to learn in the case where two jets differ only by a global rotation.
However, such a network still needs to learn a completely different set of features to recognize various relative orientations between the hard and soft particles.
On the other hand, an equivariant network might learn a single feature representing ``two subjets'' and one feature for ``diffuse lobe'', and by construction understands that these two features independently have a rotational degree of freedom.
Having learned such features, the equivariant model would easily detect either element of the substructure at arbitrary angles, and can pass this information to deeper layers within the network.
Subsequent layers could then proceed to execute computations which reason about either the absolute or relative angular position of the features.

This is precisely the intuition behind the rotational convolution: features in the jet may be matched by evaluating the projection of the jet onto a series of filters representing learned features.
These projections are sampled along a range of rotational orientations, resulting in a periodic ``waveform'' representing each individual filter's response as a function of angle.
At this point, the network has formed a representation consisting of a discrete 1D waveform with multiple channels corresponding to the different filters.

\section{Particle Projection}
\label{sec:particle-projection}

In order to define the Particle Convolution, we begin by examining the structure of the ``Deep Sets''-based Energy Flow and Particle Flow networks (EFN and PFN).
These networks operate on a jet $S$ represented by a set $S = \{s_1, \dots, s_{|S|}\}$ of observations $s_i$.
We shall assume the observations $s_i = (\vec{x}_i, q_i)$ are composed of a 2D coordinate $\vec{x}_i$ representing the direction of particles relative to the jet axis in the $\eta - \phi$ plane\footnote{To simplify notation, we will always assume jets have been centered in the $\eta - \phi$ plane, so that the coordinate $\vec{x} = (\Delta \eta, \Delta \phi) = (\eta, \phi)$.}, and some non-coordinate quantities $q_i$.
The non-coordinate observable $q_i=(e_i, \dots)$ is composed of at least an energy $e_i$ and may be supplemented by additional observables such as charge, mass, etc.

The EFN depends only on the particle coordinate $\vec{x}_i$ and energy $e_i$, and is defined as:
\begin{equation}
\label{eq:efn-def}
    \mathrm{EFN}(S) = F\left( \sum_{k=1}^{|S|} e_k \Phi(\vec{x}_k) \right) \,,
\end{equation}
where $F$ and $\Phi$ are arbitrary continuous functions, generally represented by neural networks.

Since the learnable function $\Phi$ is defined over (approximately) Euclidean spatial coordinates, we might consider it as an element of the Hilbert space $L^2(\mathbb{R}^2)$.
If we likewise consider the values of the non-coordinate observables of the jet $S$ to be represented as a function on $\mathbb{R}^2$, we can ``expand'' the jet in the continuous basis $\left| \vec{x} \right>$.
That is, for each collection $S$, we define an associated representation $\left| Z \right>_S$ such that
\begin{equation}
    \left<\vec{x} | Z \right>_S = \mathcal{Z}_S(\vec{x}) := \sum_{k=1}^{|S|} Z(q_k) \delta(\vec{x} - \vec{x}_k) \,,
\end{equation}
where $Z$ denotes some arbitrary function of the non-coordinate observable(s) $q$.
That is, in the coordinate basis, $\left| Z \right>_S$ corresponds to a function $\mathcal{Z}_S : \mathbb{R}^2 \rightarrow \mathbb{R}$ which represents the spatial distribution of the quantity $Z(q)$ within the jet $S$.
This function, $\mathcal{Z}_S$, is of course not continuous, and also highly sparse.
Strictly speaking, it is not an element of $L^2(\mathbb{R}^2)$, but rather a generalized function.
Nonetheless, we shall show that this analytic expression of the jet and its observables provides a useful mathematical tool.

For instance, if we consider an arbitrary continuous function $\Phi : \mathbb{R}^2 \rightarrow \mathbb{R}$, we may exploit the definition of the delta function in order to evaluate the Hilbert space inner product
\begin{align}
\label{eq:proj-z}
    \left< \Phi, Z\right>_S &= \int d^2\vec{x}\ \Phi(\vec{x}) \mathcal{Z}_S(\vec{x}) \\
    &= \int d^2\vec{x}\ \Phi(\vec{x}) \sum_k Z(q_k) \delta(\vec{x}-\vec{x}_k) \\
    &= \sum_k Z(q_k) \Phi(\vec{x}_k) \,.
\end{align}

In the special case $Z = \mathrm{E}$ where $\mathrm{E}(q_k) = e_k$ simply returns the energy of the particle $k$, we have:
\begin{equation}
\label{eq:proj-efn}
    \left< \Phi, \mathrm{E} \right>_S = \sum_k e_k \Phi(\vec{x}_k) \,.
\end{equation}
This last equality can immediately be recognized as the inner term of Eq.~\ref{eq:efn-def}.
In other words, we may consider the EFN to be an arbitrary function $F$ acting on the \textit{projection} of a jet's empirical energy distribution $\left| \mathrm{E} \right>_S$ onto a learned filter $\left| \Phi \right>$:
\begin{equation}
    \mathrm{EFN}(S) = F\big( \left< \Phi, \mathrm{E} \right>_S \big) \,.
\end{equation}

In the following section, we will exploit this interpretation in order to define general equivariant convolution operators.
But first, we will consider the PFN from this perspective.
The PFN as defined in~\cite{Komiske:2018cqr} is given by
\begin{equation}
    \mathrm{PFN}(S) = F\left( \sum_k \Phi(\vec{x}_k, q_k) \right) \,.
\end{equation}
Comparing with Eq.~\ref{eq:efn-def}, it is clear that the EFN is a special case of this PFN, where $\Phi$ is linear in the particle energy $e_k$.
It is this linearity which guarantees IRC-safety in the EFN case.
However, it is the fact that $\Phi$ depends only on $\vec{x}$ that allows the projection Eq.~\ref{eq:proj-z} to be related to the EFN.
This is because the inner product is defined in terms of an integral with respect to a meaningful topological measure -- in this case, two dimensional Euclidean space.
As we shall see, this is important for defining convolution with respect to locally compact groups, such as rotation and translation.

Therefore, in order to extend the concept of projection to the more general case, we define a modified version of the PFN directly in terms of the projection operator of Eq.~\ref{eq:proj-z}:
\begin{align}
    \mathrm{PFN}'(S) &= F\left( \left< \Phi, Z \right>_S \right) \\
    &= F\left( \sum_k Z(q_k) \Phi(\vec{x}_k) \right) \,.
\end{align}
Here, $F$, $Z$, and $\Phi$ are all arbitrary continuous functions that could be implemented via neural networks.
This network is equivalent to a PFN where the per-particle function is required to be separable into terms depending on the coordinate and non-coordinate observables.
In principle, this represents a strictly less general model than the original PFN.
However, in this form the $\mathrm{PFN}'$ readily admits generalization via convolution, which we shall find can result in much more effective models.

\subsection{Equivariant Projections}
\label{subsec:equivariant-projections}
Before proceeding to define the Particle Convolution, we first make some observations about the potential for equivariance in the EFN and $\mathrm{PFN}'$ architectures, which are based on projection.
In particular, it is straightforward to show that whenever $\Phi$ possesses a specific form of equivariance, so does its projection.
Suppose $\Phi$ is equivariant with respect to $G$ so that $\Phi(g\cdot x) = \Pi_g \Phi(x)$ for all $x \in X$ and $g \in G$. Then, if $\left|Z\right>_S \rightarrow g\left|Z\right>_S$, we have:
\begin{align}
    \left<\Phi | Z \right>_S \rightarrow &\left<\Phi | g | Z\right>_S \\
     = &\int dx\ \Phi(g \cdot x) \mathcal{Z}_S(x) \\
     = &\int dx\ \Pi_g \Phi(x) \mathcal{Z}_S(x) \\
     = &\ \Pi_g \left< \Phi | Z \right>_S \,.
\end{align}
So we can see that the projection is also equivariant to transformations $g$ applied to the particles of $\left|Z\right>_S$.

In practice this is often of limited utility.
For example, if $\Phi$ is constant on the $\eta - \phi$ plane, it is invariant \textit{w.r.t.} translations, and clearly so is the projection.
If $\Phi$ has some specific periodic behavior under rotations, so that in polar coordinates $\Phi(r, \theta + \delta) = \Phi(r, \theta)$ for some $\delta$, then the projection will also have this periodicity.
However, such filters are able to express only limited pattern-matching ability, and may not lead to sufficiently complex representations for the task at hand.

In the following section, we will see that by working with convolutions, arbitrary filters can be learned while retaining equivariance.

\section{Particle Convolution}
\label{sec:rotational-convolution}

Having re-cast the core operation of Particle Flow Networks in terms of the geometric concept of projection, we can immediately generalize the network by promoting projection to convolution.
In this section, we demonstrate this concretely with the example of rotational convolution.

The rotational Particle Convolution between a jet $\left|Z\right>_S$ and filter $\left<\Phi\right|$ is defined as the function:
\begin{align}
    h(\Delta;\ \Phi, Z, S) :=&\ [\Phi \star Z ]_S(\Delta) \\
    =& \left< \Phi | R_\Delta | Z \right>_S \\
    =& \int d^2\vec{x}\ \Phi(R_\Delta \vec{x}) \mathcal{Z}_S(\vec{x}) \\
    =& \sum_k Z(q_k) \Phi(R_\Delta \vec{x}_k) \,,
\end{align}
where $\Delta$ is the angle about the jet axis, and $R_\Delta \in SO(2)$ is the corresponding rotation operator.
We give the convolution the handle $h(\Delta)$, omitting the independent arguments $(\Phi,Z,S)$ when convenient, to emphasize that the result of the convolution is a \textit{function} of angle.

\begin{figure*}
    \centering
    \includegraphics[width=\textwidth]{figures/pconv.pdf}
    \caption{Illustration of the particle convolution of an example jet, represented as a point cloud, and two different filters.}
    \label{fig:pconv}
\end{figure*}

Let us examine the equivariance of this operation with respect to a rotation of the input particles $\left| Z \right>_S$:
\begin{align}
    |Z\rangle_S \rightarrow &R_\delta | Z \rangle_S \implies \\
    h(\Delta) \rightarrow &\langle \Phi | R_\Delta \big[\ R_\delta | Z\rangle_S \big] \\
    =&\langle \Phi | R_{\Delta + \delta} | Z \rangle_S \\
    =&h(\Delta + \delta) \\
    =&T_{-\delta}\ h(\Delta) \,,
\end{align}
where $T$ is the coordinate shift operator acting on the function $h$, defined by $T_y f(x) = f(x-y)$.

Since this is true for every $R_\delta \in SO(2)$ and every $\Delta$, we have established the equivariance of the convolution $h(\Delta)$ with respect to the group $SO(2)$ acting on the jet $\left| Z \right>_S$.
Moreover, we see that rotation of the input particles corresponds to a shift of the output convolution.

This convolution operation comprises the first layer of the rPCN.
The network parameters of this layer are encoded by the functions $\Phi$ and $Z$, which could for instance be themselves neural networks.
Generally, we would like to pass the result of this first layer to additional layers in a deep neural network.
However, it is not clear what to do with a continuous function.

In practice, the convolution can be \textit{sampled} at $n$ discrete points $\Delta_i = 2\pi i/n$.
In this case, we can represent the convolution $h$ as a tensor with index $i$ representing the sampled points of the convolutional ``waveform'':
\begin{equation}
\label{eq:conv-sample}
    h_i = \left<\Phi | R_{\Delta_i} | Z \right>_S = \sum_k Z_S(q_k) \Phi(R_{\Delta_i} \vec{x}_i) \,.
\end{equation}
The sampled convolution is now, strictly speaking, equivariant with respect to the subgroup of discrete rotations by $2\pi i/n$:
\begin{align}
    \left|Z\right>_S \rightarrow &R_{\frac{2\pi j}{n}} \left| Z \right>_S \implies \\
    h_i \rightarrow &\left< \Phi \right| R_{\frac{2 \pi (i+j)}{n}} \left| Z \right>_S = h_{i+j} \,.
\end{align}
Hence, a discrete rotation of the particles in $S$ corresponds to a cyclic permutation of the indices $h_{i+j}$.

By sampling a larger number of points $n$, the network can approximate continuous equivariance.
Conversely, by limiting the number of samples in accordance with the Shannon-Nyquist theorem, the network can be designed so as to ignore high frequency information which might be considered noise.

Having obtained a shift-equivariant tensor $h_i$, it is straightforward to build deeper equivariant representations by apply the standard 1D CNN layer.
The discrete 1D convolution layer is shift equivariant; however, care must be taken to also enforce the cyclic boundary conditions.
This can be done by appropriate padding of the inputs: for a convolution with kernel size $k$, the tensor $h = (h_1, \dots, h_n)$ should be extended as follows:
\begin{equation}
    h^+ = ( h_{(n-k+1)/2}, \dots, h_{n} \Big| h_1, \dots, h_n \Big| h_1, \dots, h_{(k-1)/2} ) \,.
\end{equation}
Note that in contrast to conventional CNNs, the number of samples in the convolution output is not reduced, but rather stays the same due to the periodic boundary condition.
However, it is possible to reduce the tensor along its sample axis via downsampling.

After processing by any number of additional convolutional layers, an invariant representation may be formed by a global pooling operation.
For example, it is clear that taking the maximum or average value of $h_i$ along the sample axis yields an invariant quantity.
After the pooling operation, each filter channel yields a single quantity describing an invariant feature of the input jet.
Once the invariant is formed, any additional functions applied will also result in invariants.
For example, the collection of invariant filter responses may be passed to densely connected layers, and finally to whatever output is suitable for the task objective.

\section {Steerable Convolutions}
\label{sec:steerable-convolutions}

In Sec.~\ref{sec:rotational-convolution}, we defined the rotational Particle Convolution operation.
This convolution can in principle be sampled by applying a series of rotations to the particle coordinates $R_{\Delta_i}\vec{x}_k$, and re-evaluating $\Phi$ in the projection of Eq.~\ref{eq:conv-sample}.
In practice, this can be problematic as a sizeable neural network may be required to model $\Phi$.
For example, in some experiments described in Sec.~\ref{sec:experiments}, we consider architectures in which $\Phi$ must be applied to up to 150 particles, and re-sampled at up to 21 orientations.
Therefore, even at a moderate batch size of 64 jets, a single network layer of 128 units results in a 32-bit tensor occupying nearly one GiB in memory.
This means that models based on direct sampling of convolutions cannot scale well due to memory limitations of current GPUs, and are also very time-intensive to train and optimize.

In this section, we show how the convolution may be implemented more efficiently using \textit{steerable functions}\cite{cheng1998theory,Weiler2018LearningSF}.
A function is said to be steerable if it can be expressed as a linear combination of equivariant functions\cite{cheng1998theory}.
By structuring a PCN such that the learnable filters $\Phi$ are expressed in an appropriate equivariant basis, it is possible to efficiently sample convolutions at arbitrary points by evaluating $\Phi$ only once.

This works by imposing some particular structure on our learnable filters.
Again, we will demonstrate with the specific case of rotation.
Consider for example, a filter $\psi$ of the form:
\begin{equation}
\label{eq:phi-m-def}
    \phi_m(\vec{x}) = \rho_m(r) e^{i m \theta} \,,
\end{equation}
where $m$ is an arbitrary integer, $\rho_m$ is a arbitrary function, and $r$ and $\theta$ are polar coordinates about the jet axis in the $\eta - \phi$ plane.
For any filter of this form, we can see that the rotation operator acts as:
\begin{equation}
\label{eq:phi-m-equivariance}
    \phi_m(R^{-1}_{\Delta} \vec{x}) = \rho_m(r) e^{i m (\theta - \Delta)} = e^{-i m \Delta} \phi_m(\vec{x}) \,.
\end{equation}

We can easily sample the convolution with this function at any point $\Delta$ in terms of the un-rotated projection:
\begin{align}
    h(\Delta; \phi_m, Z, S) &= \left<\phi_m | R_{\Delta} | Z \right>_S \\
    &= \sum_k \phi_m(R_\Delta \vec{x}_k) Z(q_k) \\
    &= \sum_k e^{i m \Delta} \phi_m(\vec{x}_k) Z(q_k) \\
    &= e^{i m \Delta} \left<\phi_m | Z \right>_S \,.
\end{align}
In other words, we need only compute $h(0)$ once and then use $h(\Delta) = e^{im\Delta} h(0)$.
This is, of course, an example of an equivariant projection as described in Sec.~\ref{subsec:equivariant-projections}.

Unfortunately, the filter $\phi_m$ has a definite angular frequency $m$.
As mentioned in Sec.~\ref{subsec:equivariant-projections}, the pattern-matching ability of such a filter is quite limited.
However, this suggests we could exploit a similar type of equivariance by imposing a specific structure related to Eq.~\ref{eq:phi-m-def} to construct a more general $\Phi$.


\begin{figure*}[h!]
    \centering
    \includegraphics[width=\textwidth]{figures/series.pdf} 
    \caption{
    Visualization of the series representation of learned filter instances $\Phi(\vec{x}) = \sum_m \phi_m(r, \theta)$ for various values of the cutoff hyperparameter $M$.
    }
    \label{fig:series}
\end{figure*}

In particular, let us re-define $\Phi$ in terms of a series of functions $\phi_m$:
\begin{equation}
\label{eq:phi-series}
    \Phi(\vec{x}) = \Phi(r, \theta) = \sum_{m=-M}^{M} \phi_m(r, \theta) = \sum_m \rho_m(r) e^{i m \theta} \,.
\end{equation}
Here, the hyperparameter $M$ represents a cutoff for angular frequencies captured by the filter.
The $2M+1$ radial functions $\rho_m$ are now the arbitrary learnable components of $\Phi$, which could be represented by ordinary neural networks, radial basis functions, \textit{etc}.
The relationship between $\Phi$ and this hyperparameter is depicted in Fig~\ref{fig:series}.

Note that in Eq.~\ref{eq:phi-series}, $\Phi$ is generally complex-valued, and so are the functions $\rho_m$.
For our purposes, we shall constrain it to be real-valued by imposing the condition $\rho_{-m} = \rho_{m}^\star$.
Alternatively, we could define $\Phi$ in terms of sines and cosines; however, the present form allows for simpler mathematical manipulation.

Each of the component functions $\phi_m$ are rotationally equivariant via Eq.~\ref{eq:phi-m-equivariance}.
However the filter $\Phi$, which combines arbitrarily many frequency modes, does not.
The trick is to instead associate $\Phi$ with a column vector
\begin{equation}
\Tilde{\Phi} = (\phi_{m}, ..., \phi_{-m})^\mathrm{T} \,,
\end{equation}
related via:
\begin{equation}
    \Phi(\vec{x}) = \mathbf{1}^\mathrm{T} \Tilde{\Phi}(\vec{x}) = \sum_m \phi_m(\vec{x}) \,.
\end{equation}
Now under a rotation $R_\Delta$, we have:
\begin{align}
    \Tilde{\Phi}(R^{-1}_\Delta \vec{x}) &= \big(\phi_M(R^{-1}_\Delta \vec{x}), \dots, \phi_{-M}(R^{-1}_\Delta \vec{x})\big)^\mathrm{T} \\
    &= \big( e^{-iM\Delta}\phi_M(\vec{x}), \dots, e^{iM\Delta}\phi_{-M}(\vec{x}) \big)^\mathrm{T} \\
    &= \mathbf{A}(\Delta) \Tilde{\Phi}(\vec{x}) \,,
\end{align}
where the square matrix
\begin{equation}
    \mathbf{A}(\Delta) = \mathrm{diag}\{e^{-iM\Delta}, \dots, e^{iM\Delta} \} \,.
\end{equation}
In other words, while the function $\Phi$ is not equivariant, the vector of functions $\Tilde{\Phi}$ \textit{is}, transforming as $\Tilde{\Phi} \rightarrow \mathbf{A}(\Delta)\Tilde{\Phi}$ under a rotation of coordinates by $\Delta$.

Therefore, in order to perform a convolution with $\Phi$, the network should compute the projection:
\begin{equation}
    \langle \Tilde{\Phi} | Z \rangle_S := \big(\left< \phi_M | Z\right>_S, \dots, \left< \phi_{-M} | Z\right>_S\big)^\mathrm{T} \,,
\end{equation}
once for the un-rotated case, and sample additional rotations at points $\Delta_i$ via
\begin{equation}
\label{eq:h-interpolation}
    \Tilde{h}_i = \langle\Tilde{\Phi} | R(\Delta_i) | Z \rangle_S = \mathbf{A}(\Delta_i) \langle\Tilde{\Phi}|Z\rangle_S \,.
\end{equation}
Finally, the sampled vectors $\Tilde{h}_i$ can be ``collapsed'' into a single numerical value by summing over the $m$-components:
\begin{equation}
    h_i = \mathbf{1}^\mathrm{T} \Tilde{h}_i
\end{equation}
Note that given a fixed set of sample points $\Delta_i$, the operators $\mathbf{A}(\Delta_i)$ are simply constant, diagonal matrices, and the matrix multiplication of Eq.~\ref{eq:h-interpolation} is a trivial operation for GPUs.

The network based on steerable convolutions is depicted schematically in Figure~\ref{fig:flowchart}.

\begin{figure}[ht]
    \centering
    \includegraphics[width=0.5\textwidth]{figures/flowchart.pdf}
    \caption{Schematic representation of the rotational Particle Convolution layer as described in Sec.~\ref{sec:steerable-convolution}.
    The functions $z_c$ and $\rho_{cm}$ can be specifically chosen (e.g. when setting $z_c = p_\mathrm{T}$ for IRC-safety), or can be implemented as neural networks.
    }
    \label{fig:flowchart}
\end{figure}

Having re-formulated our filter $\Phi$ in terms of the equivariant functions $\phi_m$, we can revisit the question of efficiency and scalability for practical networks.
Note that in either case, for a filter $\Phi$ to capture information about given angular frequency $M$, we require $2M+1$ samples for $h_i$.
In Eq.~\ref{eq:conv-sample}, we do this by re-evaluating the projection with $\Phi$ $2M+1$ times.
In Eq.~\ref{eq:h-interpolation}, we evaluate the projection only once, but there are $2M+1$ functions $\phi_m$ which must be evaluated.
Therefore, it may seem that nothing has been gained.

However, in practice, the functions $\rho_m$ can be collectively implemented by a single neural network which shares most of its weights.
As long as the number of samples required, $2M+1$, is generally less than the size of hidden layers for typical dense networks, then the network should have a substantially smaller memory footprint.
This is because the intermediate tensors for the hidden layers need only be computed for the un-rotated case, as opposed to being calculated at every orientation as in Eq.~\ref{eq:conv-sample}.
This allows for smaller networks which are faster to train.

Moreover, the functions $\rho_m$ only need to learn a radial profile, rather than a response on the full 2-dimensional $\eta - \phi$ plane.
Therefore, it is reasonable to expect that simpler, smaller networks could be able to achieve equivalent results.
Lastly, since $\Phi$ in Eq.~\ref{eq:phi-series} is expressed in terms of functions with a definite cutoff frequency $M$, this acts as low-pass filter for angular structure.
Without this cutoff, higher frequency information may be aliased to lower modes during sampling, which could lead to a network sensitive to undesirable artifacts.

We find empirically that in either the direct sampling or steerable case, rPCNs of comparable angular resolution and depth achieve similar performance.
Therefore, when optimizing hyperparameters for the experiments in section~\ref{sec:experiments}, the more efficient steerable convolution architecture is used; however, it is possible that a more exhaustive optimization would find better results via the direct sampling approach.


\section{Experiments}
\label{sec:experiments}
We conduct experiments in training the rotational Particle Convolution Network for two benchmark problems in jet physics: quark/gluon identification, and top tagging.
The details of these datasets are described in Sec.~\ref{subsec:datasets}.
For both tasks, we train a general PCN as well as the IRC-safe variant.
The coordinates $\vec{x}$ are centered about the jet axis in the $\eta - \phi$ plane, as described in~\cite{Komiske:2018cqr}, and the scalar feature for each particle is the particle's transverse momentum, \textit{i.e.} $q_k = p_{\mathrm{T},k}$, which is invariant w.r.t. the coordinate-centering operation.
In the case of the q/g tagging task, we also test the performance of models when supplied with per-particle identification information as well, $q_k = (p_{\mathrm{T},k}, \mathrm{PID}_k)$.

\subsection{Datasets}
\label{subsec:datasets}

For reference, we consider two typical benchmark problems: q/g identification and top-quark tagging, using the same datasets as in Refs.~\cite{Komiske:2018cqr,Qu:2019gqs}.
The q/g dataset, described in Ref.\cite{Komiske:2018cqr} is comprised of $Z(\nu\nu) + q$ events (signal) and $Z(\nu\nu)+g$ events (background).
The events are simulated using \textsc{Pythia}8\cite{Sjostrand:2014zea}, which performs parton showering and hadronization.
No detector simulation is performed for this sample.
Particles (excluding neutrinos) are clustered using the anti-$k_\mathrm{T}$ algorithm\cite{Cacciari:2008gp} with radius parameter $R=0.4$.
The 4-momenta and particle ID of particles within the leading jet are saved for jets with $p_\mathrm{T}\in[500, 550]$ GeV and $|\eta|<2$.
The dataset is split into 1.2 million training events and 400k each of validation and test events.
In our experiments, we truncate events with greater than 68 particles by discarding those with the lowest $p_\mathrm{T}$.

For the top-tagging benchmark, we use the dataset available at Ref.\cite{kasieczka_gregor_2019_2603256}.
The events are generated with \textsc{Pythia}8 and passed into Delphes\cite{deFavereau:2013fsa} fast detector simulation, without pileup.
The signal process is $t\bar{t}$ and the background is inclusive QCD.
Jets are reconstructed using Delphes EFlow module, with anti-$k_\mathrm{T}$ radius parameter $R=0.8$, and are required to have $p_\mathrm{T} \in [550, 650]$ GeV and $|\eta|<2$.
For each jet, we record the 4-momenta for the leading 140 particle-flow constituents.



\subsection{Architectures}
In previous sections, the discussion has focused on projection and convolution of a jet with a single filter $\Phi$.
In practice, we allow a neural network to learn a moderate number of independent filter channels $\Phi_c(r, \theta) = \sum_m \rho_{cm}(r) e^{i m \theta}$ and feature embeddings $Z_q$.
By convolving each feature embedding with each filter, the network produces a convolution waveform $\langle \Phi_c | R(\Delta) | Z_j\rangle$ with $C \times J$ feature channels, which can be passed on as multi-feature input to a standard 1D CNN architecture.
Hence, to define the architecture of a particular rPCN, we must to specify the following:
\begin{itemize}
    \item The number of filter channels, $C$;
    \item The number of feature embeddings, $J$;
    \item The maximum angular frequency mode, $M$;
    \item The feature embedding function(s), $Z_j(q)$;
    \item The radial activation function(s), $\rho_{cm}(r)$;
    \item The remaining 1D convolutional network structure.
\end{itemize}

Except in the IRC-safe case, we implement both $Z_j$ and $\rho_{cm}$ as densely-connected neural networks.
Therefore, $Z_j$ and $\rho_{cm}$ are specified by the number of hidden layers and their units, as well as their nonlinearies.
We apply a ReLU nonlinearity at the output of $Z_c$, so that it can more naturally act as an ``attention'' mechanism which can learn to ignore particles based on their scalar features $q_k$ in the region where the ReLU response is zero.
The complex-valued network $\rho_{cm}$ is implemented by a single network which learns real and imaginary parts separately, $\rho_{cm} = \alpha_{cm} + i \beta_{cm}$.
We impose the constraint $\rho_{-m} = \rho_m^*$ so that $\Phi_c$ are real by construction; therefore it is sufficient to learn $\alpha_{cm}$ for $m\geq 0$ and $\beta_{cm}$ for $m\geq 1$.

When particle ID is included as part of the scalar input, we follow the ``experimentally plausible'' labeling scheme of \cite{Komiske:2018cqr}, where particles are categorized into one of eight types: photons, neutral hadrons, and positively- and negatively-charged muons, electrons, and hadrons.
These labels are input to the network via a trainable 3-dimensional embedding layer, which are concatenated with the particle's $p_\mathrm{T}$ and passed as a triplet into the $Z_c$ network.

In the IRC-safe case, the function $Z_c$ is simply replaced with the transverse momentum of each particle.
Since the majority of trainable parameters in the Particle Convolution network tend to be from the $Z_c$ and $\rho_{cm}$ sub-networks, the IRC-safe networks will usually have substantially fewer parameters than their non-IRC-safe counterparts.

In our experiments, the remainder of the network can be specified by:
\begin{itemize}
    \item The nonlinearity following the Particle Convolution layer;
    \item The number of 1D convolution layers, their kernel size, and stride;
    \item The global pooling operation (maximum or average);
    \item Any remaining hidden layers, units, and nonlinearity;
    \item A final 2-unit softmax output layer representing the categorical prediction.
\end{itemize}

We found that a resnet-like~\cite{2015arXiv151203385H} convolutional architecture worked best.
We pass the particle convolution output without nonlinearity as a ``skip connection'' across residual blocks consisting of ReLU activation and batch normalization preceding two 1D CNN layers.

\subsection{Training}
In all experiments, the data are split into train, validation, and test samples in a 6:1:1 ratio.
Networks are implemented using Tensorflow v2~\cite{tensorflow2015-whitepaper} and Keras~\cite{chollet2015keras}.
The loss function used for training is the categorical crossentropy, optimized via Adam\cite{2014arXiv1412.6980K} with learning rate $10^{-4}$ and batch size $128$.
The validation loss is used to determine training convergence; training is stopped when the validation loss has not improved for $16$ consecutive epochs.
The model epoch with the lowest validation loss is evaluated on the test set to report unbiased AUC and rejection metrics.

After training networks by randomly sampling a wide range of hyperparameter configurations, we selected the best-performing architecture for each task based on the validation loss.
These configurations were then frozen, and the networks were retrained ten times.
The values reported in tables~\ref{tab:comparison-qg} and~\ref{tab:comparison-top} are the test scores for the specific network whose validation score was nearest to the median value on each task.
The errors quoted are the standard deviation of the test scores over the series of ten retrainings.

\subsection{Inspection}
It is possible to visualize the learned functions $\Phi_c(\vec{x})$ in the $\eta - \phi$ plane.
Moreover, we can visualize the learned behavior of the nonlinear feature embedding, $Z_j(p_\mathrm{T})$.
Each convolution represents a learned spatial pattern as well as a gated response based on momentum.
Interestingly, it seems that the network often spontaneously learns a sort of binning in $p_\mathrm{T}$, as shown in Figure~\ref{fig:z-response}.
It also tends to ignore many of the lowest-$p_\mathrm{T}$ particles, suggesting those particles could potentially be pruned from the input to further speed up the network.

\begin{figure}
    \centering
    \includegraphics[width=0.5\textwidth]{figures/z-vs-pt.pdf}
    \caption{Colored lines representing functions $Z(p_\mathrm{T})$ learned by a network on the q/g tagging task. The shaded histogram indicates the distribution of particle $p_\mathrm{T}$ over the training dataset.}
    \label{fig:z-response}
\end{figure}

\section{Results \& Conclusions}
\label{sec:conclusion}

The rPCN performance on the two benchmark tasks are given in Tables~\ref{tab:comparison-qg} and \ref{tab:comparison-top}.
These results show that by promoting the permutation-invariant architectures EnergyFlow and ParticleFlow to include rotational convolution significantly improves their performance.
In particular, the rPCN significantly improves the state-of-the-art for IRC-safe taggers.
When IRC safety can be dispensed with, the rPCN can achieve performance comparable to the current state-of-the-art set by graph-based models such as ParticleNet.

From a practical perspective, the rPCN is much more similar to deep-sets based architectures already in production use by LHC experiments, as compared to the relatively novel graph architectures.
The rPCN model also requires minimal preprocessing, is compatible with arbitrary-length collections of particles, and is permutation-equivariant.
PCNs in general are conceptually analogous to the jet image approach in the limit of infinitely small pixels and continuous filter kenels, although the rPCN targets rotational rather than translational equivariance as a physically-motivated inductive bias.

Moreover, by careful design of the filter representation $\Phi$, the PCN approach can admit further generalization by constructing convolutional operations that can endow additional equivariance properties to jet tagging models.
For example, in forthcoming work we study whether a PCN which is equivariant \textit{w.r.t.} to scaling in the $\eta - \phi$ plane could potentially be useful to capture similar substructure properties across a wide range of jet $p_\mathrm{T}$.

\section{Acknowledgements}
The author would like to especially thank Paul Tipton for many iterations of feedback throughout this project.
He thanks Aishik Ghosh, Daniel Whiteson, Dan Guest, and Ema Smith for helpful conversations and comments on this manuscript.
Training experiments were made possible by the Grace GPU cluster operated by Yale Center for Research Computing. 
This work was supported by grant DE-SC0017660 funded by the U.S. Department of Energy, Office of Science.



\begin{table*}[ht]
    \centering
    \begin{tabular}{l r c c}
    \vspace{0.25 em}
    & Model & AUC & $1/\epsilon_b |_{\epsilon_s = 50\%}$ \\
    \hline
    \hline
    \multicolumn{1}{l}{\textit{IRC-safe}} & & & \\
    & EFN & 0.8824 & $28.6 \pm 0.3$ \\
    & rPCN \textbf{(Ours)} & \textbf{0.8944} & $\mathbf{32.5 \pm 0.4}$ \\
    \hline
    \multicolumn{1}{l}{\textit{w/o PID}} & & & \\
    & PFN & 0.8911 & 30.8 $\pm 0.4$ \\
    & ResNeXt-50 & 0.8960 & 30.9 \\
    & P-CNN & 0.8915 & 31.0 \\
    & ParticleNet-Lite & 0.8993 & 32.8 \\
    & ParticleNet & \textbf{0.9014} & 33.7 \\
    & rPCN \textbf{(Ours)} & 0.8997 & $\mathbf{34.2 \pm 0.4}$ \\
    \hline
    \multicolumn{1}{l}{\textit{w/ PID}} & & & \\
    & P-CNN & 0.9002 & 34.7 \\
    & PFN-Ex & 0.9005 & 34.7 \\
    & ParticleNet-Lite & 0.9079 & 37.1 \\
    & rPCN +Ex \textbf{(Ours)} & 0.9081 & $38.6 \pm 0.5$ \\
    & ParticleNet & \textbf{0.9116} & $\mathbf{39.8 \pm 0.2}$ \\
    \hline
    \end{tabular}
    \caption{
    Comparison of known network architectures on the quark-gluon tagging dataset of Ref.~\cite{Komiske:2018cqr}.
    ResNeXt is a deep 2D CNN adapted for jet images in~\cite{Qu:2019gqs}, while P-CNN~\cite{CMS-DP-2017-049} is a 1D convolution operating on ordered lists of particles, also implemented in~\cite{Qu:2019gqs}.
    Our IRC-safe rPCN model is substantially more sensitive than the EFN.
    When including a nonlinear $p_\mathrm{T}$ function, the rPCN model achieves comparable performance to ParticleNet.
    }
    \label{tab:comparison-qg}
\end{table*}


\begin{table*}[ht]
    \centering
    \begin{tabular}{l  r c c c}
    \vspace{0.25 em}
    & Model & AUC & $1/\epsilon_b |_{\epsilon_s = 50\%}$ & $1/\epsilon_b|_{\epsilon_s = 30\%}$ \\
    \hline
    \hline
    \multicolumn{1}{l}{\textit{IRC-safe}} & & & & \\
    & EFPs & 0.9803 & 184 & 384 \\
    & EFN & 0.9789 & 181 & 619 \\
    & rPCN \textbf{(Ours)} &  $\mathbf{0.9821}$ & $\mathbf{257 \pm 6}$ & $\mathbf{1038 \pm 41}$ \\
    \hline
    \multicolumn{1}{l}{\textit{w/o PID}} & & & & \\
    & P-CNN & 0.9803 & 201 & 759 \\
    & PFN & 0.9819 & 247 & 888 \\
    & ResNeXt-50 & 0.9837 & 302 & 1147 \\
    & ParticleNet-Lite & 0.9844 & 325 $\pm 5$ & 1262 $\pm 49$ \\
    & rPCN \textbf{(Ours)} & 0.9845 & 364 $\pm 9$ & $\mathbf{1642 \pm 93}$ \\ 
    & ParticleNet & $\mathbf{0.9858}$ & $\mathbf{397 \pm 7}$ & $1615 \pm 93$ \\
    \hline
    \end{tabular}
    \caption{
    Comparison of similar network architectures on the top tagging dataset of Ref.~\cite{Komiske:2018cqr}.
    ResNeXt is a deep 2D CNN adapted for jet images in~\cite{Qu:2019gqs}, while P-CNN~\cite{CMS-DP-2017-049} is a 1D convolution operating on ordered lists of particles, also implemented in~\cite{Qu:2019gqs}.
    As with the q/g task, the rPCN sets a new state-of-the-art for IRC-safe tagging, while achieving comparable performance to the graph-based ParticleNet model.
    }
    \label{tab:comparison-top}
\end{table*}




\newpage

\bibliographystyle{unsrt}  
