# Group-invariant max filtering

Jameson Cahill111Department of Mathematics and Statistics, University of North
Carolina Wilmington, Wilmington, NC Joseph W. Iverson222Department of
Mathematics, Iowa State University, Ames, IA Dustin G. Mixon333Department of
Mathematics, The Ohio State University, Columbus, OH 444Translational Data
Analytics Institute, The Ohio State University, Columbus, OH Daniel
Packer33footnotemark: 3

###### Abstract

Given a real inner product space $V$ and a group $G$ of linear isometries, we
construct a family of $G$-invariant real-valued functions on $V$ that we call
max filters. In the case where $V=\mathbb{R}^{d}$ and $G$ is finite, a
suitable max filter bank separates orbits, and is even bilipschitz in the
quotient metric. In the case where $V=L^{2}(\mathbb{R}^{d})$ and $G$ is the
group of translation operators, a max filter exhibits stability to
diffeomorphic distortion like that of the scattering transform introduced by
Mallat. We establish that max filters are well suited for various
classification tasks, both in theory and in practice.

##  1 Introduction

Modern machine learning has been extraordinarily successful in domains where
large volumes of labeled data are available [48, 65]. Indeed, highly
expressive models can generalize once they fit an appropriately large training
set. Unfortunately, many important domains are plagued by a scarcity of data
or by expensive labels (or both). One way to bridge this gap is by augmenting
the given dataset with the help of a large family of innocuous distortions. In
many cases, the distortions correspond to the action of a group, meaning the
ground truth exhibits known symmetries. Augmenting the training set by
applying the group action encourages the model to learn these symmetries.
While this approach has been successful [66, 28, 48, 26], it is extremely
inefficient to train a large, symmetry-agnostic model to find a highly
symmetric function. One wonders:

Why not use a model that already accounts for known symmetries?

This motivates invariant machine learning (e.g., [76, 78, 41, 69, 6]), where
the model is invariant to underlying symmetries in the data. To illustrate,
suppose an object is represented by a point $x$ in a set $V$, but there is a
group $G$ acting on $V$ such that the same object is also represented by
$gx\in V$ for every $g\in G$. This ambiguity emerges, for example, when using
a matrix to represent a point cloud or a graph, since the representation
depends on the labeling of the points or vertices. If we apply a $G$-invariant
feature map $\Phi\colon V\to F$, then the learning task can be performed in
the feature domain $F$ without having to worry about symmetries in the
problem. Furthermore, if $\Phi$ separates the $G$-orbits in $V$, then no
information is lost by passing to the feature domain. In practice, $V$ and $F$
tend to be vector spaces out of convenience, and $G$ is frequently a linear
group.

While our interest in invariants stems from modern machine learning, maps like
$\Phi$ have been studied since Cayley established invariant theory in the
nineteenth century [25]. Here, we take $V=\mathbb{C}^{d}$ and
$G\leq\operatorname{GL}(V)$, and the maps of interest consist of the
$G$-invariant polynomials $\mathbb{C}[V]^{G}$. In 1890, Hilbert [42] proved
that $\mathbb{C}[V]^{G}$ is finitely generated as a $\mathbb{C}$-algebra in
the special case where $G$ is the image of a representation of
$\operatorname{SL}(\mathbb{C}^{k})$, meaning one may take the feature domain
$F$ to be finite dimensional (one dimension for each generator). Since $G$ is
not a compact subset of $\mathbb{C}^{d\times d}$ in such cases, there may
exist distinct $G$-orbits whose closures intersect, meaning no continuous
$G$-invariant function can separate them; this subtlety plays an important
role in Mumford’s more general geometric invariant theory [55]. In general,
the generating set of $\mathbb{C}[V]^{G}$ is often extraordinarily large [36],
making it impractical for machine learning applications. To alleviate this
issue, there has been some work to construct separating sets of polynomials
[32, 33, 34], i.e., sets that separate as well as $\mathbb{C}[V]^{G}$ does
without necessarily generating all of $\mathbb{C}[V]^{G}$. For every reductive
group $G$, there exists a separating set of $2d+1$ invariant polynomials [36,
19], but the complexity of evaluating these polynomials is still quite large.
Furthermore, these polynomials tend to have high degree, and so they are
numerically unstable in practice. In practice, one also desires a quantitative
notion of separating so that distant orbits are not sent to nearby points in
the feature space, and this behavior is not always afforded by a separating
set of polynomials [19]. Despite these shortcomings, polynomial invariants are
popular in the data science literature due in part to their rich algebraic
theory, e.g., [8, 58, 19, 11, 6].

In this paper, we focus on the case where $V$ is a real inner product space
and $G$ is a group of linear isometries of $V$. We introduce a family of non-
polynomial invariants that we call max filters. In Section 2, we define max
filters, we identify some basic properties, and we highlight a few familiar
examples. In Section 3, we use ideas from [37] to establish that $2d$ generic
max filters separate all $G$-orbits when $G$ is finite (see Corollary 13), and
then we describe various settings in which max filtering is computationally
efficient. In Section 4, we show that when $G$ is finite, a sufficiently large
random max filter bank is bilipschitz with high probability; see Theorem 18.
This is the first known construction of invariant maps for a general class of
groups that enjoy a lower Lipschitz bound, meaning they separate orbits in a
quantitative sense. In the same section, we later show that when
$V=L^{2}(\mathbb{R}^{d})$ and $G$ is the group of translations, certain max
filters exhibit stability to diffeomorphic distortion akin to what Mallat
established for his scattering transform in [52]; see Theorem 22. In Section
5, we explain how to select max filters for classification in a couple of
different settings, we determine the subgradient of max filters to enable
training, and we characterize how random max filters behave for the symmetric
group. In Section 6, we use max filtering to process real-world datasets.
Specifically, we visualize the shape space of voting districts, we use
electrocardiogram data to classify whether patients had a heart attack, and we
classify a multitude of textures. Surprisingly, we find that even in cases
where the data do not appear to exhibit symmetries in a group $G$, max
filtering with respect to $G$ can still reveal salient features. We conclude
in Section 7 with a discussion of opportunities for follow-on work.

##  2 Preliminaries

Given a real inner product space $V$ and a group $G$ of linear isometries of
$V$, consider the quotient space $V/G$ consisting of the $G$-orbits
$[x]:=G\cdot x$ for $x\in V$. This quotient space is equipped with the metric

$d([x],[y]):=\inf_{p\in[x],q\in[y]}\|p-q\|.$

(Indeed, $d$ satisfies the triangle inequality since $G$ is a group of
isometries of $V$.) This paper is concerned with the following function.

######  Definition 1.

The max filtering map $\langle\langle\cdot,\cdot\rangle\rangle\colon V/G\times
V/G\to\mathbb{R}$ is defined by

$\langle\langle[x],[y]\rangle\rangle:=\sup_{p\in[x],q\in[y]}\langle
p,q\rangle.$

Sometimes, “max filtering map” refers to the related function
$\langle\langle[\cdot],[\cdot]\rangle\rangle\colon V\times V\to\mathbb{R}$
instead. The intended domain should be clear from context.

We note that since $G$ consists of linear isometries, it is closed under
adjoints, and so the max filtering map can be alternatively expressed as

$\langle\langle[x],[y]\rangle\rangle=\sup_{g\in G}\langle
gx,y\rangle=\sup_{g\in G}\langle x,gy\rangle.$

Furthermore, if $G$ is topologically closed (e.g., finite), then the supremum
can be replaced with a maximum:

$\langle\langle[x],[y]\rangle\rangle=\max_{g\in G}\langle
gx,y\rangle=\max_{g\in G}\langle x,gy\rangle.$

The max filtering map satisfies several important properties, summarized
below.

######  Lemma 2.

Suppose $x,y,z\in V$. Then each of the following holds:

  * (a)

$\langle\langle[x],[x]\rangle\rangle=\|x\|^{2}$.

  * (b)

$\langle\langle[x],[y]\rangle\rangle=\langle\langle[y],[x]\rangle\rangle$.

  * (c)

$\langle\langle[x],[ry]\rangle\rangle=r\langle\langle[x],[y]\rangle\rangle$
for every $r\geq 0$.

  * (d)

$\langle\langle[x],[\cdot]\rangle\rangle\colon V\to\mathbb{R}$ is convex.

  * (e)

$\langle\langle[x],[y+z]\rangle\rangle\leq\langle\langle[x],[y]\rangle\rangle+\langle\langle[x],[z]\rangle\rangle$.

  * (f)

$d([x],[y])^{2}=\|x\|^{2}-2\langle\langle[x],[y]\rangle\rangle+\|y\|^{2}$.

  * (g)

$\langle\langle[x],\cdot\rangle\rangle\colon V/G\to\mathbb{R}$ is
$\|x\|$-Lipschitz.

###### Proof.

First, (a) and (b) are immediate, as is the $r=0$ case of (c). For $r>0$,
observe that $q\in[ry]$ precisely when $q^{\prime}:=r^{-1}q\in[y]$, since each
member of $G$ is a linear isometry. Thus,

$\langle\langle[x],[ry]\rangle\rangle=\sup_{p\in[x],q\in[ry]}\langle
p,q\rangle=\sup_{p\in[x],q^{\prime}\in[y]}\langle
p,rq^{\prime}\rangle=r\cdot\sup_{p\in[x],q^{\prime}\in[y]}\langle
p,q^{\prime}\rangle=r\langle\langle[x],[y]\rangle\rangle.$

Next, (d) follows from the identity
$\langle\langle[x],[y]\rangle\rangle=\sup_{p\in[x]}\langle p,y\rangle$, which
expresses $\langle\langle[x],[\cdot]\rangle\rangle$ as a pointwise supremum of
convex functions. For (e), we apply (d) and (c):

$\langle\langle[x],[y+z]\rangle\rangle=\langle\langle[x],[\tfrac{1}{2}\cdot
2y+\tfrac{1}{2}\cdot
2z]\rangle\rangle\leq\tfrac{1}{2}\langle\langle[x],[2y]\rangle\rangle+\tfrac{1}{2}\langle\langle[x],[2z]\rangle\rangle=\langle\langle[x],[y]\rangle\rangle+\langle\langle[x],[z]\rangle\rangle.$

Next, (f) is immediate. For (g), select any $y,z\in V$. We may assume
$\langle\langle[x],[y]\rangle\rangle\geq\langle\langle[x],[z]\rangle\rangle$
without loss of generality. Select any $p\in[y]$, $q\in[z]$, and $\epsilon>0$,
and take $g\in G$ such that $\langle
x,gp\rangle>\langle\langle[x],[y]\rangle\rangle-\epsilon$. Then

$\displaystyle|\langle\langle[x],[y]\rangle\rangle-\langle\langle[x],[z]\rangle\rangle|$
$\displaystyle=\langle\langle[x],[y]\rangle\rangle-\langle\langle[x],[z]\rangle\rangle$
$\displaystyle<\langle x,gp\rangle+\epsilon-\langle x,gq\rangle=\langle
x,g(p-q)\rangle+\epsilon\leq\|x\|\|p-q\|+\epsilon.$

Since $p$, $q$, and $\epsilon$ were arbitrary, the result follows. ∎

######  Definition 3.

Given a template $z\in V$, we refer to
$\langle\langle[z],\cdot\rangle\rangle\colon V/G\to\mathbb{R}$ as the
corresponding max filter. Given a (possibly infinite) sequence
$\\{z_{i}\\}_{i\in I}$ of templates in $V$, the corresponding max filter bank
is $\Phi\colon V/G\to\mathbb{R}^{I}$ defined by
$\Phi([x]):=\\{\langle\langle[z_{i}],[x]\rangle\rangle\\}_{i\in I}$.

In what follows, we identify a few familiar examples of max filters.

######  Example 4 (Norms).

There are several norms that can be thought of as a max filter with some
template. For example, consider $V=\mathbb{R}^{n}$. Then taking
$G=\operatorname{O}(n)$ and any unit-norm template $z$ gives

$\langle\langle[z],[x]\rangle\rangle=\max_{g\in\operatorname{O}(n)}\langle
gz,x\rangle=\|x\|.$

Similarly, the infinity norm is obtained by taking $G$ to be the group of
signed permutation matrices and $z$ to be a standard basis element, while the
$1$-norm comes from taking $G$ to be the group of diagonal orthogonal matrices
and $z$ to be the all-ones vector. We can also recover various matrix norms
when $V=\mathbb{R}^{m\times n}$. For example, taking
$G\cong\operatorname{O}(m)\times\operatorname{O}(n)$ to be the group of linear
operators of the form $X\mapsto Q_{1}XQ_{2}^{-1}$ for
$Q_{1}\in\operatorname{O}(m)$ and $Q_{2}\in\operatorname{O}(n)$, then max
filtering with any rank-$1$ matrix of unit Frobenius norm gives the spectral
norm.

######  Example 5 (Power spectrum).

Consider the case where $V=L^{2}(\mathbb{R}/\mathbb{Z})$ and $G$ is the group
of circular translation operators $T_{a}$ defined by $T_{a}g(t):=g(t-a)$ for
$a\in\mathbb{R}$. (Here and throughout, functions in $L^{2}$ will be real
valued by default.) Given a template $z_{k}$ of the form $z_{k}(t):=\cos(2\pi
kt)$ for some $k\in\mathbb{N}$, it holds that

$\displaystyle\langle\langle[z_{k}],[f]\rangle\rangle=\max_{a\in[0,1)}\langle
T_{a}z_{k},f\rangle$ $\displaystyle=\max_{a\in[0,1)}\int_{0}^{1}\cos(2\pi
k(t-a))f(t)dt$ $\displaystyle=\max_{a\in[0,1)}\operatorname{Re}\bigg{(}e^{2\pi
ika}\int_{0}^{1}f(t)e^{-2\pi ikt}dt\bigg{)}=|\hat{f}(k)|.$

A similar choice of templates recovers the power spectrum over finite abelian
groups.

######  Example 6 (Unitary groups).

While we generally assume $V$ is a real inner produce space, our theory also
applies in the complex setting. For example, consider the case where
$V=\mathbb{C}^{n}$ and $G\leq\operatorname{U}(n)$. Then $V$ is a
$2n$-dimensional real inner product space with

$\langle x,y\rangle:=\operatorname{Re}(x^{*}y),$

where $x^{*}$ denotes the conjugate transpose of $x$. As such,
$\operatorname{U}(n)\leq\operatorname{O}(V)$ since $g\in\operatorname{U}(n)$
implies

$\langle
gx,gy\rangle=\operatorname{Re}((gx)^{*}(gy))=\operatorname{Re}(x^{*}g^{*}gy)=\operatorname{Re}(x^{*}y)=\langle
x,y\rangle.$

Thus, $G\leq\operatorname{O}(V)$.

######  Example 7 (Phase retrieval).

Suppose $V=\mathbb{C}^{r}$ and
$G=\\{c\cdot\operatorname{id}:|c|=1\\}\leq\operatorname{U}(r)$. Then

$\langle\langle[z],[x]\rangle\rangle=\max_{|c|=1}\operatorname{Re}((cz)^{*}x)=|z^{*}x|.$

The max filter bank corresponding to $\\{z_{i}\\}_{i=1}^{n}$ in $V$ is given
by $\Phi([x])=\\{|z_{i}^{*}x|\\}_{i=1}^{n}$. The inverse problem of recovering
$[x]$ from $\Phi([x])$ is known as complex phase retrieval [4, 9, 29, 70], and
over the last decade, several algorithms were developed to solve this inverse
problem [22, 31, 20, 72, 21, 27]. In the related setting where
$V=\mathbb{R}^{d}$ and $G=\\{\pm\operatorname{id}\\}\leq\operatorname{O}(d)$,
the analogous inverse problem is known as real phase retrieval [4].

######  Example 8 (Matched filtering).

In classical radar, the primary task is to locate a target. Here, a
transmitter emits a pulse $p\in L^{2}(\mathbb{R})$, which then bounces off the
target and is received at the transmitter’s location with a known direction of
arrival. The return signal $q$ is a noisy version of $T_{a}p$ for some $a>0$,
where $T_{a}$ denotes the translation-by-$a$ operator defined by
$T_{a}f(t):=f(t-a)$. Considering the transmitter-to-target distance is $a/2$
times the speed of light, the objective is to estimate $a$, which can be
accomplished with matched filtering: simply find $a$ for which $\langle
T_{a}p,q\rangle$ is largest. This is essentially a max filter with
$V=L^{2}(\mathbb{R})$ and $G$ being the group of translation operators, though
for this estimation problem, the object of interest is the maximizer $a$, not
the maximum value $\langle\langle[p],[q]\rangle\rangle$. Meanwhile, the
maximum value is used for the detection problem of distinguishing noise from
noisy versions of translates of $p$. (This accounts for half of the etymology
of max filtering.)

######  Example 9 (Max pooling).

In a convolutional neural network, it is common for a convolutional layer to
be followed by a max pooling layer. Here, the convolutional layer convolves
the input image with several localized templates, and then the max pooling
layer downsamples each of the resulting convolutions by partitioning the scene
into patches and recording the maximum value in each patch. In the extreme
case where the max pooling layer takes the entire scene to be a single patch
to maximize over, these layers implement a max filter bank in which $V$ is the
image space and $G$ is the group of translation operators. (This accounts for
the other half of the etymology of max filtering.)

##  3 The complexity of separating orbits

For practical reasons, we are interested in orbit-separating invariants
$\Phi\colon V\to\mathbb{R}^{n}$ of low complexity, which we take to mean two
different things simultaneously:

  * (i)

$n$ is small (i.e., the map has low sample complexity), and

  * (ii)

one may evaluate $\Phi$ efficiently (i.e., the map has low computational
complexity).

While these notions of complexity are related, they impact the learning task
in different ways. In what follows, we study both notions of complexity in the
context of max filtering.

###  3.1 Generic templates separate orbits

In this subsection, we focus on the case in which $V=\mathbb{R}^{d}$ and
$G\leq\operatorname{O}(d)$ is semialgebraic, which we will define shortly.
Every polynomial function $p\colon\mathbb{R}^{n}\to\mathbb{R}$ determines a
basic semialgebraic set

$\\{x\in\mathbb{R}^{n}:p(x)\geq 0\\}.$

By closing under finite unions, finite intersections, and complementation, the
basic semialgebraic sets of $\mathbb{R}^{n}$ generate an algebra of sets known
as the semialgebraic sets in $\mathbb{R}^{n}$. A semialgebraic subgroup
$G\leq\operatorname{GL}(d)$ is a subgroup of $\operatorname{GL}(d)$ that is
also a semialgebraic set in $\mathbb{R}^{d\times d}$, e.g.,
$\operatorname{O}(d)$. A semialgebraic function is a function
$f\colon\mathbb{R}^{s}\to\mathbb{R}^{t}$ for which the graph
$\\{(x,f(x)):x\in\mathbb{R}^{s}\\}$ is a semialgebraic set in
$\mathbb{R}^{s+t}$.

######  Lemma 10.

For every semialgebraic subgroup $G\leq\operatorname{O}(d)$, the corresponding
max filtering map
$\langle\langle[\cdot],[\cdot]\rangle\rangle\colon\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}$
is semialgebraic.

###### Proof.

The graph of the max filtering map can be expressed in first-order logic:

$\\{(x,y,t)\in\mathbb{R}^{d}\times\mathbb{R}^{d}\times\mathbb{R}:(\forall g\in
G,t\geq\langle x,gy\rangle)\wedge(\forall\epsilon\in\mathbb{R},\exists g\in
G,\epsilon>0\Rightarrow t-\epsilon<\langle x,gy\rangle)\\}.$

(To be precise, one should replace our quantifiers over $G$ with the
polynomial conditions that define the semialgebraic set $G$ to obtain a
condition in first-order logic.) It follows from Proposition 2.2.4 in [12]
that the graph is semialgebraic. ∎

Every semialgebraic set $A$ can be decomposed as a disjoint union
$A=\bigcup_{i}A_{i}$ of finitely many semialgebraic sets $A_{i}$, each of
which is homeomorphic to an open hypercube $(0,1)^{d_{i}}$ (where $(0,1)^{0}$
is a point). The dimension of $A$ can be defined in terms of this
decomposition as $\operatorname{dim}(A):=\max_{i}d_{i}$. (It does not depend
on the decomposition.)

######  Definition 11.

Given a semialgebraic subgroup $G\leq\operatorname{O}(d)$, we say the
corresponding max filtering map
$\langle\langle[\cdot],[\cdot]\rangle\rangle\colon\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}$
is $k$-strongly separating if for every $x,y\in\mathbb{R}^{d}$ with
$[x]\neq[y]$, it holds that

$\operatorname{dim}\big{\\{}z\in\mathbb{R}^{d}:\langle\langle[z],[x]\rangle\rangle=\langle\langle[z],[y]\rangle\rangle\big{\\}}\leq
d-k.$

As an example, consider the case where $G=\operatorname{O}(d)$. Then
$\langle\langle[z],[x]\rangle\rangle=\langle\langle[z],[y]\rangle\rangle$
holds precisely when $\|z\|\|x\|=\|z\|\|y\|$, i.e., $z=0$ or $[x]=[y]$. Thus,
the max filtering map is $d$-strongly separating in this case.

######  Theorem 12.

Consider any semialgebraic subgroup $G\leq\operatorname{O}(d)$ with
$k$-strongly separating max filtering map
$\langle\langle[\cdot],[\cdot]\rangle\rangle\colon\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}$
for some $k\in\mathbb{N}$. For generic $z_{1},\ldots,z_{n}\in\mathbb{R}^{d}$,
the max filter bank
$x\mapsto\\{\langle\langle[z_{i}],[x]\rangle\rangle\\}_{i=1}^{n}$ separates
$G$-orbits in $\mathbb{R}^{d}$ provided $n\geq 2d/k$.

In the $d$-strongly separating case where $G=\operatorname{O}(d)$, Theorem 12
implies that $n=2$ generic templates suffice to separate orbits. (Of course,
any single nonzero template suffices in this case.) Theorem 12 is an
improvement to Theorem 1.9 in [37], which gives the condition $n\geq 2d+1$. We
obtain the improvement $n\geq 2d/k$ by leveraging a more detailed notion of
strongly separating, as well as the positive homogeneity of max filtering. The
proof makes use of a lift-and-project technique that first appeared in [4] and
was subsequently applied in [29, 73, 19, 62].

###### Proof of Theorem 12.

Fix $n\geq 2d$, and let $\mathcal{Z}\subseteq(\mathbb{R}^{d})^{n}$ denote the
set of $\\{z_{i}\\}_{i=1}^{n}$ for which the max filter bank
$x\mapsto\\{\langle\langle[z_{i}],[x]\rangle\rangle\\}_{i=1}^{n}$ fails to
separate $G$-orbits in $\mathbb{R}^{d}$. We will show that $\mathcal{Z}$ is
semialgebraic with dimension $\leq dn-1$, from which the result follows. To do
so, observe that $\\{z_{i}\\}_{i=1}^{n}\in\mathcal{Z}$ precisely when there
exists a witness, namely, $(x,y)\in\mathbb{R}^{d}\times\mathbb{R}^{d}$ with
$[x]\neq[y]$ such that
$\langle\langle[z_{i}],[x]\rangle\rangle=\langle\langle[z_{i}],[y]\rangle\rangle$
for every $i\in\\{1,\ldots,n\\}$. In fact, we may assume that the witness
$(x,y)$ satisfies $\|x\|^{2}+\|y\|^{2}=1$ without loss of generality since the
set of witnesses for $\\{z_{i}\\}_{i=1}^{n}$ avoids $(0,0)$ and is closed
under positive scalar multiplication by Lemma 2(c). This suggests the
following lift of $\mathcal{Z}$:

$\displaystyle\mathcal{L}:=\Big{\\{}~{}(\\{z_{i}\\}_{i=1}^{n},(x,y))\in(\mathbb{R}^{d})^{n}\times(\mathbb{R}^{d})^{2}~{}:$
$\displaystyle~{}[x]\neq[y],~{}\|x\|^{2}+\|y\|^{2}=1,$
$\displaystyle~{}\langle\langle[z_{i}],[x]\rangle\rangle=\langle\langle[z_{i}],[y]\rangle\rangle~{}~{}\forall
i\in\\{1,\ldots,n\\}~{}\Big{\\}}.$

Since $G$ is semialgebraic, we have that $[x]\neq[y]$ is a semialgebraic
condition. Furthermore,
$\langle\langle[z_{i}],[x]\rangle\rangle=\langle\langle[z_{i}],[y]\rangle\rangle$
is a semialgebraic condition for each $i$ by Lemma 10. It follows that
$\mathcal{L}$ is semialgebraic. Next, we define the projection maps
$\pi_{1}\colon(\\{z_{i}\\}_{i=1}^{n},(x,y))\mapsto\\{z_{i}\\}_{i=1}^{n}$ and
$\pi_{2}\colon(\\{z_{i}\\}_{i=1}^{n},(x,y))\mapsto(x,y)$. Then
$\mathcal{Z}=\pi_{1}(\mathcal{L})$ is semialgebraic by Tarski–Seidenberg
(Proposition 2.2.1 in [12]). To bound the dimension of $\mathcal{Z}$, we first
observe that

$\pi_{2}^{-1}(x,y)=\big{\\{}z\in\mathbb{R}^{d}:\langle\langle[z],[x]\rangle\rangle=\langle\langle[z],[y]\rangle\rangle\big{\\}}^{n}\times\\{(x,y)\\},$

and so $\operatorname{dim}(\pi_{2}^{-1}(x,y))\leq n(d-k)$, since the max
filtering map is $k$-strongly separating by assumption. We use the fact that
$\pi_{2}(\mathcal{L})$ is contained in the unit sphere in
$(\mathbb{R}^{d})^{2}$ together with Lemma 1.10 in [37] to obtain

$\operatorname{dim}(\mathcal{Z})\leq\operatorname{dim}(\mathcal{L})\leq\operatorname{dim}(\pi_{2}(\mathcal{L}))+\max_{x,y\in\mathbb{R}^{d}}\operatorname{dim}(\pi_{2}^{-1}(x,y))\leq(2d-1)+n(d-k)\leq
dn-1,$

where the last step is equivalent to the assumption $n\geq 2d/k$. ∎

######  Corollary 13.

Consider any finite subgroup $G\leq\operatorname{O}(d)$. For generic
$z_{1},\ldots,z_{n}\in\mathbb{R}^{d}$, the max filter bank
$x\mapsto\\{\langle\langle[z_{i}],[x]\rangle\rangle\\}_{i=1}^{n}$ separates
$G$-orbits in $\mathbb{R}^{d}$ provided $n\geq 2d$.

Corollary 13 follows immediately from Theorem 12 and the following lemma:

######  Lemma 14.

For every finite subgroup $G\leq\operatorname{O}(d)$, the corresponding max
filtering map
$\langle\langle[\cdot],[\cdot]\rangle\rangle\colon\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}$
is $1$-strongly separating.

###### Proof.

Consider any $x,y\in\mathbb{R}^{d}$ with $[x]\neq[y]$. Then
$\langle\langle[z],[x]\rangle\rangle=\langle\langle[z],[y]\rangle\rangle$ only
if there exists $g\in G$ such that $\langle z,x\rangle=\langle z,gy\rangle$,
i.e., $z\in\operatorname{span}\\{x-gy\\}^{\perp}$. Thus,

$\big{\\{}z\in\mathbb{R}^{d}:\langle\langle[z],[x]\rangle\rangle=\langle\langle[z],[y]\rangle\rangle\big{\\}}\subseteq\bigcup_{g\in
G}\operatorname{span}\\{x-gy\\}^{\perp}.$

Since the max filtering map is semialgebraic by Lemma 10, it follows that the
left-hand set is also semialgebraic. Since $G$ is finite, the right-hand set
is semialgebraic with dimension $d-1$, and the result follows. ∎

We would like to know if a version of Corollary 13 holds for all semialgebraic
groups, but we do not have a proof of strongly separating for infinite groups
in general. This motivates the following problem:

######  Problem 15.

  * (a)

For which semialgebraic groups is the max filtering map $k$-strongly
separating?

  * (b)

How many templates are needed to separate orbits for a given group?

We identify a couple of interesting instances of Problem 15. First, we
consider the case of complex phase retrieval (as in Example 7), where
$V=\mathbb{C}^{r}$ and $G$ is the center of $\operatorname{U}(r)$. It is known
that $n=4r-4=2\operatorname{dim}(V)-4$ generic templates separate orbits for
every $r$, and this is the optimal threshold for infinitely many $r$ [29], but
there also exist $11$ templates in $\mathbb{C}^{4}$ that separate orbits, for
example [70].

As another example, consider the case where $V=\mathbb{R}^{d}$ and $G\cong
S_{d}$ is the group of $d\times d$ permutation matrices. Then Corollary 13
gives that $2d$ generic templates separate orbits. However, it is
straightforward to see that the templates $z_{j}:=\sum_{i=1}^{j}e_{i}$ for
$j\in\\{1,\ldots,d\\}$ also separate orbits, where $e_{i}$ denotes the $i$th
standard basis element. Indeed, take $\operatorname{sort}(x)$ to have weakly
decreasing entries. Then the first entry equals
$\langle\langle[z_{1}],[x]\rangle\rangle$, while for each $j>1$, the $j$th
entry equals
$\langle\langle[z_{j}],[x]\rangle\rangle-\langle\langle[z_{j-1}],[x]\rangle\rangle$.
As such, this max filter bank determines $\operatorname{sort}(x)$, which is a
separating invariant of $V/G$. Considering Theorem 12, one might suspect that
the max filtering map is $2$-strongly separating in this case, but this is not
so. Indeed, the cone $C$ of sorted vectors in $\mathbb{R}^{d}$ has dimension
$d$, and so there exists a subspace $H$ of co-dimension $1$ that intersects
the interior of $C$. Select any $x$ in the interior of $C$ and any unit vector
$v\in H^{\perp}$, and then take $y=x+\epsilon v$ for $\epsilon>0$ sufficiently
small so that $y\in C$. Then

$\big{\\{}z\in\mathbb{R}^{d}:\langle\langle[z],[x]\rangle\rangle=\langle\langle[z],[y]\rangle\rangle\big{\\}}\supseteq
H\cap C,$

which is a semialgebraic set of dimension $d-1$, and so the claim follows.

###  3.2 Low-complexity max filtering

In this subsection, we focus on the case in which $V\cong\mathbb{R}^{d}$.
Naively, one may compute the max filtering map
$(x,y)\mapsto\langle\langle[x],[y]\rangle\rangle$ over a finite group
$G\leq\operatorname{O}(d)$ of order $m$ by computing $\langle x,gy\rangle$ for
every $g\in G$ and then returning the maximum. This approach costs $O(md)$
operations. Of course, this is not possible when $G$ is infinite, and it is
prohibitive when $G$ is finite but large. Interestingly, many of the groups
that we encounter in practice admit a faster implementation. In particular,
for many quotient spaces $V/G$, the quotient metric $d\colon V/G\times
V/G\to\mathbb{R}$ is easy to compute, and by Lemma 2(f), the max filtering map
is equally easy to compute in these cases:

$\langle\langle[x],[y]\rangle\rangle=\frac{1}{2}\Big{(}d([x],[y])^{2}-\|x\|^{2}-\|y\|^{2}\Big{)}.$

In this subsection, we highlight a few examples of such quotient spaces before
considering the harder setting of graphs.

####  3.2.1 Point clouds

Consider $V=\mathbb{R}^{k\times n}$ with Frobenius inner product. We can
represent a point cloud of $n$ points in $\mathbb{R}^{k}$ as a member of $V$
by arbitrarily labeling the points with column indices. In this setting, we
identify members of $V$ that reside in a common $G$-orbit with $G\cong S_{n}$
permuting the columns. The resulting quotient metric is known as the
$2$-Wasserstein distance:

$d([X],[Y])=\min_{P\in\Pi(n)}\|X-YP\|_{F},$

where $\Pi(n)$ denotes the set of $n\times n$ permutation matrices and
$\|\cdot\|_{F}$ denotes the Frobenius norm. The corresponding max filtering
map is then given by

$\langle\langle[X],[Y]\rangle\rangle=\max_{P\in\Pi(n)}\langle
X,YP\rangle=\max_{P\in\Pi(n)}\operatorname{tr}(X^{\top}YP)=\max_{S\in\operatorname{conv}\Pi(n)}\operatorname{tr}(X^{\top}YS),$

where $\operatorname{conv}\Pi(n)$ denotes the convex hull of $\Pi(n)$, namely,
the doubly stochastic matrices. By this formulation, the max filtering map can
be computed in polynomial time by linear programming. In the special case
where $k=1$, the max filtering map has an even faster implementation:

$\langle\langle[X],[Y]\rangle\rangle=\langle\operatorname{sort}(X),\operatorname{sort}(Y)\rangle,$

which can be computed in linearithmic time.

####  3.2.2 Circular translations

Consider the case where $V\cong\mathbb{R}^{n}$ is the space of vectors with
entries indexed by the cyclic group $C_{n}:=\mathbb{Z}/n\mathbb{Z}$, and
$G\cong C_{n}$ is the group of circular translations $T_{a}$ defined by
$T_{a}f(x):=f(x-a)$. Then the max filtering map is given by

$\langle\langle[f],[g]\rangle\rangle=\max_{a\in C_{n}}\langle
f,T_{a}g\rangle=\max_{a\in C_{n}}\sum_{x\in C_{n}}f(x)g(x-a)=\max_{a\in
C_{n}}(f\star Rg)(a),$

where $\star$ denotes the circular convolution and $R$ denotes the reversal
operator. Thus, the max filtering map can be computed in linearithmic time
with the help of the fast Fourier transform.

####  3.2.3 Shape analysis

In geometric morphometics [53], it is common for data to take the form of a
sequence of $n$ landmarks in $\mathbb{R}^{k}$ (where $k$ is typically $2$ or
$3$) with a global rotation ambiguity. This corresponds to taking
$V=\mathbb{R}^{k\times n}$ and $G\cong\operatorname{O}(k)$ acting on the left,
and so the max filtering map is given by

$\langle\langle[X],[Y]\rangle\rangle=\max_{R\in\operatorname{O}(k)}\langle
X,RY\rangle=\max_{R\in\operatorname{O}(k)}\operatorname{tr}(RYX^{\top})=\|YX^{\top}\|_{*},$

where $\|\cdot\|_{*}$ denotes the nuclear norm. As such, the max filtering map
can be computed in polynomial time with the aid of the singular value
decomposition.

####  3.2.4 Separation hierarchy for weighted graphs

Here, we focus on the case in which $V$ is the vector space of real symmetric
$n\times n$ matrices with zero diagonal and $G\cong S_{n}$ is the group of
linear isometries of the form $A\mapsto PAP^{-1}$, where $P$ is a permutation
matrix. We think of $V/G$ as the space of weighted graphs on $n$ vertices (up
to isomorphism). One popular approach for separating graphs uses message-
passing graph neural networks, but the separation power of such networks is
limited by the so-called Weisfeler–Lehman test [75, 77, 54]. For example,
message-passing graph neural networks fail to distinguish $C_{3}\cup C_{3}$
from $C_{6}$. See [64, 43] for surveys of this rapidly growing literature.

As an alternative, we consider a max filtering approach. Given two adjacency
matrices $A_{1}$ and $A_{2}$, Lemma 2(f) implies that the corresponding graphs
are isomorphic if and only if

$\|A_{1}\|_{F}^{2}=\langle\langle[A_{1}],[A_{2}]\rangle\rangle=\|A_{2}\|_{F}^{2}.$

As such, max filtering is graph isomorphism–hard in this setting.
Interestingly, there exist $A\in V$ for which the map
$X\mapsto\langle\langle[A],[X]\rangle\rangle$ can be computed in linearithmic
time, and furthermore, these easy-to-compute max filters help with separating
orbits. To see this, we follow [2], which uses color coding to facilitate
computation by dynamic programming.

######  Definition 16.

A tuple $\\{f_{i}\\}_{i=1}^{N}$ with $f_{i}\colon[n]\to[k]$ for each $i\in[N]$
is an $(n,k)$-color coding if $n\geq k$ and for every $S\subseteq[n]$ of
cardinality $k$, there exists $i\in[N]$ such that $f_{i}(S)=[k]$.

######  Lemma 17.

Given $n,k\in\mathbb{N}$ with $n\geq k$, there exists an $(n,k)$-color coding
of size $\lceil ke^{k}\log n\rceil$.

###### Proof.

We show that $N$ random colorings form a color coding with positive
probability. For each $i\in[N]$ and $S\in\binom{[n]}{k}$, we have
$\mathbb{P}\\{f_{i}(S)=[k]\\}=k!/k^{k}$, and so the union bound gives

$\mathbb{P}\Big{\\{}\\{f_{i}\\}_{i=1}^{N}\text{ is not $(n,k)$-color
coding}\Big{\\}}\leq\sum_{S\in\binom{[n]}{k}}\prod_{i=1}^{N}\mathbb{P}\\{f_{i}(S)\neq[k]\\}=\binom{n}{k}\bigg{(}1-\frac{k!}{k^{k}}\bigg{)}^{N}.$

It suffices to select $N$ so that the right-hand side is strictly smaller than
$1$. The result follows by applying the bounds $\binom{n}{k}\leq n^{k}$,
$k!\geq(k/e)^{k}$, and $(1-1/t)^{t}<1/e$ for $t>1$:

$\tbinom{n}{k}(1-\tfrac{k!}{k^{k}})^{N}\leq
n^{k}(1-e^{-k})^{e^{k}e^{-k}N}<e^{k\log n-e^{-k}N},$

which is at most $1$ when $N\geq ke^{k}\log n$. ∎

Data: Weighted tree with vertex set $[k]$ (in post-order traversal order) and
with adjacency matrix $A\in\mathbb{R}^{k\times k}$, and $(n,k)$-color coding
$\\{f_{i}\\}_{i=1}^{N}$

Input: Weighted graph with adjacency matrix $B\in\mathbb{R}^{n\times n}$

Result: Max filter $\langle\langle[\tilde{A}],[B]\rangle\rangle$, where
$\tilde{A}:=[\begin{smallmatrix}A&0\\\
0&0\end{smallmatrix}]\in\mathbb{R}^{n\times n}$

for  _$i\in[N]$ and $\pi\in S_{k}$_ do

Initialize $\ell(v)\leftarrow 0$ for all $v\in[n]$

for  _$u\in[k-1]$_ do

Take the unique vertex $u^{\prime}\in\\{u+1,\ldots,k\\}$ adjacent to the leaf
$u$ in $H-[u-1]$

for  _$v^{\prime}\in(\pi\circ f_{i})^{-1}(u^{\prime})$_ do

$\ell(v^{\prime})\leftarrow\ell(v^{\prime})+\max_{v\in(\pi\circ
f_{i})^{-1}(u)}(\ell(v)+A_{u,u^{\prime}}B_{v,v^{\prime}})$

end for

end for

Put $s(i,\pi):=\max_{v\in(\pi\circ f_{i})^{-1}(k)}\ell(v)$

end for

Output $\max_{i\in[N],\pi\in S_{k}}s(i,\pi)$

Algorithm 1 Max filtering with a weighted tree template by color coding

Algorithm 1 computes the max filter with a small weighted tree using a color
coding and dynamic programming. Lemma 17 implies that Algorithm 1 has runtime
$e^{O(k\log k)}n^{2}\log n$, which is linearithmic in the size of the data
when $k$ is fixed. Notice that max filtering with the path on $k=4$ vertices
already separates the graphs $C_{3}\cup C_{3}$ and $C_{6}$. Furthermore, using
techniques from [2], one can modify Algorithm 1 to max filter with any
template graph $H$ on $k$ vertices, though the runtime becomes $e^{O(k\log
k)}n^{t+1}\log n$, where $t$ is the treewidth of $H$. Letting
$\mathcal{H}(k,t)$ denote the set of weighted graphs on at most $k$ vertices
with treewidth at most $t$, we have the following hierarchy:

$\begin{array}[]{ccccccccc}\mathcal{H}(n,1)&\subseteq&\cdots&\subseteq&\mathcal{H}(n,n-1)\\\
\rotatebox[origin={c}]{90.0}{$\subseteq$}\\\ \vdots&&\iddots\\\
\rotatebox[origin={c}]{90.0}{$\subseteq$}\\\ \mathcal{H}(2,1)\end{array}$ (1)

Corollary 13 gives that $n(n-1)$ generic templates from $\mathcal{H}(n,n-1)$
separate all isomorphism classes of weighted graphs on $n$ vertices. It would
be interesting to study the separation power of templates of logarithmic order
and bounded treewidth.

##  4 Stability of max filtering

###  4.1 Bilipschitz max filter banks

Upper and lower Lipschitz bounds are used to quantify the stability of a
mapping between metric spaces, but it is generally difficult to estimate such
bounds; see [9, 7, 17, 45, 5] for examples from phase retrieval and [79, 19,
18, 6] for other examples. In this subsection, we prove the following:

######  Theorem 18.

Fix a finite group $G\leq\operatorname{O}(d)$ of order $m$ and select

$n\geq
12m^{2}d\log(\tfrac{2}{\delta}+1),\qquad\delta:=(\tfrac{\pi}{128m^{4}}\cdot\tfrac{1}{2d+3\log(4m^{2})})^{1/2}.$

Draw independent random vectors
$z_{1},\ldots,z_{n}\sim\mathsf{Unif}(S^{d-1})$. With probability $\geq
1-e^{-n/(12m^{2})}$, it holds that the max filter bank
$\Phi\colon\mathbb{R}^{d}/G\to\mathbb{R}^{n}$ with templates
$\\{z_{i}\\}_{i=1}^{n}$ has lower Lipschitz bound $\delta$ and upper Lipschitz
bound $n^{1/2}$.

This result distinguishes max filtering from separating polynomial invariants,
which do not necessarily enjoy upper or lower Lipschitz bounds [19]. In
Theorem 18, we may take the embedding dimension to be $n=\Theta^{*}(m^{2}d)$
with bilipschitz bounds $\Theta^{*}(\frac{1}{m^{2}d^{1/2}})$ and
$\Theta^{*}(md^{1/2})$, where $\Theta^{*}(\cdot)$ suppresses logarithmic
factors. For comparison, we consider a couple of cases that have already been
studied in the literature. First, the case where
$G=\\{\pm\operatorname{id}\\}$ reduces to the setting of real phase retrieval
(as in Example 7), where it is known that there exist $n=\Theta(d)$ templates
that deliver lower- and upper-Lipschitz bounds $\frac{1}{4}$ and $4$, say; see
equation (17) in [9]. Notably, these bounds do not get worse as $d$ gets
large. It would be interesting if a version of Theorem 18 held for infinite
groups, but we do not expect it to hold for infinite-dimensional inner product
spaces. Case in point, for $V=\ell^{2}$ with $G=\\{\pm\operatorname{id}\\}$,
it was shown in [17] that for every choice of templates, the map is not
bilipschitz.

Another interesting phenomenon from finite-dimensional phase retrieval is that
separating implies bilipschitz; see Lemma 16 and Theorem 18 in [9] and
Proposition 1.4 in [17]. This suggests the following:

######  Problem 19.

Is every separating max filter bank
$\Phi\colon\mathbb{R}^{d}/G\to\mathbb{R}^{n}$ bilipschitz?

If the answer to Problem 19 is “yes,” then Corollary 13 implies that $2d$
generic templates produce a bilipschitz max filter bank
$\Phi\colon\mathbb{R}^{d}/G\to\mathbb{R}^{2d}$ whenever
$G\leq\operatorname{O}(d)$ is finite.

Theorem 18 follows immediately from Lemmas 20 and 21 below. Our proof uses the
following notion that was introduced in [1]. We say
$\\{z_{i}\\}_{i=1}^{n}\in(\mathbb{R}^{d})^{n}$ exhibits
$(k,\delta)$-projective uniformity if

$s_{k}\\{|\langle z_{i},x\rangle|\\}_{i=1}^{n}\geq\delta\|x\|$

for every $x\in\mathbb{R}^{d}$, where $s_{k}\colon\mathbb{R}^{n}\to\mathbb{R}$
returns the $k$th smallest entry of the input. In what follows, we denote
$\|\\{z_{i}\\}_{i=1}^{n}\|_{F}:=(\sum_{i=1}^{n}\|z_{i}\|^{2})^{1/2}$.

######  Lemma 20.

Fix a finite subgroup $G\leq\operatorname{O}(d)$ and suppose
$\\{z_{i}\\}_{i=1}^{n}\in(\mathbb{R}^{d})^{n}$ exhibits
$(\lceil\frac{n}{|G|^{2}}\rceil,\delta)$-projective uniformity. Then the max
filter bank $\Phi\colon\mathbb{R}^{d}/G\to\mathbb{R}^{n}$ with templates
$\\{z_{i}\\}_{i=1}^{n}$ has lower Lipschitz bound $\delta$ and upper Lipschitz
bound $\|\\{z_{i}\\}_{i=1}^{n}\|_{F}$.

###### Proof.

The upper Lipschitz bound follows from Lemma 2(g):

$\|\Phi([x])-\Phi([y])\|^{2}=\sum_{i=1}^{n}|\langle\langle[z_{i}],[x]\rangle\rangle-\langle\langle[z_{i}],[y]\rangle\rangle|^{2}\leq\|\\{z_{i}\\}_{i=1}^{n}\|_{F}^{2}\cdot
d([x],[y])^{2}.$

For the lower Lipschitz bound, fix $x,y\in\mathbb{R}^{d}$ with $[x]\neq[y]$,
and then for each $i\in\\{1,\ldots,n\\}$, select $g_{i},h_{i}\in G$ such that
$\langle\langle[z_{i}],[x]\rangle\rangle=\langle z_{i},g_{i}x\rangle$ and
$\langle\langle[z_{i}],[y]\rangle\rangle=\langle z_{i},h_{i}y\rangle$. Then

$\|\Phi([x])-\Phi([y])\|^{2}=\sum_{i=1}^{n}\langle
z_{i},g_{i}x-h_{i}y\rangle^{2}\geq\bigg{(}\sum_{i=1}^{n}\langle
z_{i},\tfrac{g_{i}x-h_{i}y}{\|g_{i}x-h_{i}y\|}\rangle^{2}\bigg{)}\cdot
d([x],[y])^{2},$ (2)

where the inequality follows from the bound $\|g_{i}x-h_{i}y\|\geq
d([x],[y])$. Next, consider the map $p\colon i\mapsto(g_{i},h_{i})$, and
select $(g,h)\in G^{2}$ with the largest preimage. By pigeonhole, we have
$|p^{-1}(g,h)|\geq\lceil\frac{n}{|G|^{2}}\rceil=:k$, and so

$\sum_{i=1}^{n}\langle
z_{i},\tfrac{g_{i}x-h_{i}y}{\|g_{i}x-h_{i}y\|}\rangle^{2}\geq\max_{i\in
p^{-1}(g,h)}\langle z_{i},\tfrac{gx-hy}{\|gx-hy\|}\rangle^{2}\geq
s_{k}\\{\langle z_{i},\tfrac{gx-hy}{\|gx-
hy\|}\rangle^{2}\\}_{i=1}^{n}\geq\delta^{2}.$

Combining with (2) gives the result. ∎

The following lemma gives that random templates exhibit projective uniformity.

######  Lemma 21 (cf. Lemma 6.9 in [1]).

Select $p\in(0,1)$ and take

$\delta:=(\tfrac{\pi}{128}\cdot\tfrac{p^{2}}{2d+3\log(4/p)})^{1/2},\qquad
n\geq\tfrac{12d}{p}\log(\tfrac{2}{\delta}+1).$ (3)

Draw independent random vectors
$z_{1},\ldots,z_{n}\sim\mathsf{Unif}(S^{d-1})$. Then $\\{z_{i}\\}_{i=1}^{n}$
exhibits $(\lceil pn\rceil,\delta)$-projective uniformity with probability
$\geq 1-e^{-pn/12}$.

###### Proof.

Put $k:=\lceil pn\rceil$, let $\mathcal{E}$ denote the failure event that
$\\{z_{i}\\}_{i=1}^{n}$ does not have $(k,\delta)$-projective uniformity, and
let $N_{\delta}$ denote a $\delta$-net of $S^{d-1}$ of minimum size. Note that
if $v$ is within $\delta$ of $x$, then for every $z_{i}$, it holds that

$|\langle z_{i},v\rangle|\leq|\langle z_{i},x\rangle|+|\langle
z_{i},v-x\rangle|\leq|\langle z_{i},x\rangle|+\|v-x\|\leq|\langle
z_{i},x\rangle|+\delta.$

Thus, we may pass to the $\delta$-net to get

$\displaystyle\mathbb{P}(\mathcal{E})$
$\displaystyle=\mathbb{P}\Big{\\{}\text{ $\exists x\in S^{d-1}$ s.t.\
$s_{k}\\{|\langle z_{i},x\rangle|\\}_{i=1}^{n}<\delta$ }\Big{\\}}$
$\displaystyle\leq\mathbb{P}\Big{\\{}\text{ $\exists v\in N_{\delta}$ s.t.\
$s_{k}\\{|\langle z_{i},v\rangle|\\}_{i=1}^{n}<2\delta$ }\Big{\\}}$
$\displaystyle\leq|N_{\delta}|\cdot\mathbb{P}\Big{\\{}s_{k}\\{|\langle
z_{i},e_{1}\rangle|\\}_{i=1}^{n}<2\delta\Big{\\}}=|N_{\delta}|\cdot\mathbb{P}\Big{\\{}\sum_{i=1}^{n}\mathbf{1}_{\\{|\langle
z_{i},e_{1}\rangle|<2\delta\\}}\geq k\Big{\\}},$

where the second inequality applies the union bound and the rotation
invariance of the distribution $\mathsf{Unif}(S^{d-1})$. A standard volume
comparison argument gives $|N_{\delta}|\leq(\frac{2}{\delta}+1)^{d}$. The
final probability concerns a sum of independent Bernoulli variables with some
success probability $q=q(d,\delta)$, which can be estimated using the
multiplicative Chernoff bound:

$\mathbb{P}\Big{\\{}\sum_{i=1}^{n}\mathbf{1}_{\\{|\langle
z_{i},e_{1}\rangle|<2\delta\\}}\geq
k\Big{\\}}\leq\mathbb{P}\Big{\\{}\sum_{i=1}^{n}\mathbf{1}_{\\{|\langle
z_{i},e_{1}\rangle|<2\delta\\}}\geq
pn\Big{\\}}\leq\exp(-\tfrac{(p-q)^{2}}{p+q}\cdot n),$

provided $p>q$. Next, we verify that $q(d,\delta)\leq\frac{p}{2}$. Denoting
$g\sim\mathsf{N}(0,I_{d})$, we have

$\displaystyle
q:=\mathbb{P}\\{|\langle\tfrac{g}{\|g\|},e_{1}\rangle|<2\delta\\}$
$\displaystyle\leq\inf_{t>0}\Big{(}\mathbb{P}\\{|\langle
g,e_{1}\rangle|<2\delta t\\}+\mathbb{P}\\{\|g\|^{2}>t^{2}\\}\Big{)}$
$\displaystyle\leq\inf_{t\geq\sqrt{2d}}\Big{(}\sqrt{\tfrac{2}{\pi}}\cdot
2\delta t+e^{-(t^{2}-2d)/3}\Big{)},$

where the final inequality uses the facts that $|\langle g,e_{1}\rangle|$ has
half-normal distribution and $\|g\|^{2}$ has chi-squared distribution with $d$
degrees of freedom. We select $t:=(2d+3\log(\frac{4}{p}))^{1/2}$ so that the
second term equals $\frac{p}{4}$, and then our choice (3) for $\delta$ ensures
that the first term equals $\frac{p}{4}$. Overall, we have

$\displaystyle\mathbb{P}(\mathcal{E})$
$\displaystyle\leq|N_{\delta}|\cdot\mathbb{P}\Big{\\{}\sum_{i=1}^{n}\mathbf{1}_{\\{|\langle
z_{i},e_{1}\rangle|<2\delta\\}}\geq k\Big{\\}}$
$\displaystyle\leq(\tfrac{2}{\delta}+1)^{d}\cdot\exp(-\tfrac{(p-q)^{2}}{p+q}\cdot
n)\leq\exp(d\log(\tfrac{2}{\delta}+1)-\tfrac{pn}{6})\leq e^{-pn/12},$

where the last step applied our assumption that
$n\geq\frac{12d}{p}\log(\frac{2}{\delta}+1)$. ∎

###  4.2 Mallat-type stability to diffeomorphic distortion

In this subsection, we focus on the case in which $V=L^{2}(\mathbb{R}^{d})$
and $G$ is the group of translation operators $T_{a}$ defined by
$T_{a}f(x):=f(x-a)$ for $a\in\mathbb{R}^{d}$. Given a template $h\in
L^{2}(\mathbb{R}^{d})$, the corresponding max filter is

$\langle\langle[h],[f]\rangle\rangle=\sup_{a\in\mathbb{R}^{d}}\langle
h,T_{a}f\rangle=\sup_{a\in\mathbb{R}^{d}}\int_{\mathbb{R}^{d}}h(x)f(x-a)dx=\sup_{a\in\mathbb{R}^{d}}(Rh\star
f)(a),$

where $R$ denotes the reversal operator defined by $Rh(x):=h(-x)$ and $\star$
denotes convolution. (Of course, the supremum of a member of
$L^{2}(\mathbb{R}^{d})$ is not well defined, but $Rh\star f$ is continuous
since $Rh,f\in L^{2}(\mathbb{R}^{d})$.)

Our motivation for this setting stems from image analysis, in which case
$d=2$. For a familiar example, consider the task of classifying handwritten
digits. Intuitively, each class is translation invariant, and so it makes
sense to treat images as members of $V/G$. In addition, images that are slight
elastic distortions of each other should be sent to nearby points in the
feature domain. The fact that image classification is invariant to such
distortions has been used to augment the MNIST training set and boost
classification performance [66]. Instead of using data augmentation to learn
distortion-invariant features, it is desirable to restrict to feature maps
that already exhibit distortion invariance. (Indeed, such feature maps would
require fewer parameters to train.) This compelled Mallat to introduce his
scattering transform [52], which has since played an important role in the
theory of invariant machine learning [15, 16, 71, 39, 57]. Mallat used the
following formalism to analyze the stability of the scattering transform to
distortion.

Given a diffeomorphism $g\in C^{1}(\mathbb{R}^{d})$, we consider the
corresponding distortion operator $L_{g}$ defined by
$L_{g}f(x):=f(g^{-1}(x))$. It will be convenient to interact with the vector
field $\tau:=\operatorname{id}-g^{-1}\in C^{1}(\mathbb{R}^{d})$, since
$L_{g}f(x)=f(x-\tau(x))$. For example, if $\tau(x)=a$ for every
$x\in\mathbb{R}^{d}$, then $L_{g}$ is translation by $a$. In what follows,
$J\tau(x)\in\mathbb{R}^{d\times d}$ denotes the Jacobian matrix of $\tau$ at
$x$.

######  Theorem 22.

Take any continuously differentiable $h\in L^{2}(\mathbb{R}^{d})$ for which

$(1+\|x\|_{2})^{(d+1)/2}\cdot
h(x)\qquad\text{and}\qquad(1+\|x\|_{2})^{(d+3)/2}\cdot\nabla h(x)$ (4)

are bounded. There exists $C(h)>0$ such that for every $f\in
L^{2}(\mathbb{R}^{d})$ and every diffeomorphism $g\in C^{1}(\mathbb{R}^{d})$
for which $\tau:=\operatorname{id}-g^{-1}$ satisfies
$\sup_{x\in\mathbb{R}^{d}}\|J\tau(x)\|_{2\to 2}\leq\frac{1}{2}$, it holds that

$|\langle\langle[h],[f]\rangle\rangle-\langle\langle[h],[L_{g}f]\rangle\rangle|\leq
C(h)\cdot\|f\|_{L^{2}(\mathbb{R}^{d})}\cdot\sup_{x\in\mathbb{R}^{d}}\|J\tau(x)\|_{2\to
2}.$

This matches Mallat’s bound [52] on the stability of the scattering transform
to diffeomorphic distortion, except our bound has no Hessian term. The proof
of Theorem 22 follows almost immediately from the following modification of
Lemma E.1 in [52], which bounds the commutator between the filter and the
distortion by the magnitude of the distortion:

######  Lemma 23.

Take $h\in L^{2}(\mathbb{R}^{d})$ as in Theorem 22, and consider the linear
operator $Z_{h}$ defined by $Z_{h}f:=h\star f$. There exists $C(h)>0$ such
that for every diffeomorphism $g\in C^{1}(\mathbb{R}^{d})$ for which
$\tau:=\operatorname{id}-g^{-1}$ satisfies
$\sup_{x\in\mathbb{R}^{d}}\|J\tau(x)\|_{2\to 2}\leq\frac{1}{2}$, it holds that

$\|L_{g}Z_{h}-Z_{h}L_{g}\|_{L^{2}(\mathbb{R}^{d})\to
L^{\infty}(\mathbb{R}^{d})}\leq
C(h)\cdot\sup_{x\in\mathbb{R}^{d}}\|J\tau(x)\|_{2\to 2}.$

Assuming Lemma 23 for the moment, we can prove Theorem 22.

###### Proof of Theorem 22.

The change of variables $a=g^{-1}(a^{\prime})$ gives

$\displaystyle|\langle\langle[h],[f]\rangle\rangle-\langle\langle[h],[L_{g}f]\rangle\rangle|$
$\displaystyle=\Big{|}\sup_{a\in\mathbb{R}^{d}}(Z_{Rh}f)(a)-\sup_{b\in\mathbb{R}^{d}}(Z_{Rh}L_{g}f)(b)\Big{|}$
$\displaystyle=\Big{|}\sup_{a^{\prime}\in\mathbb{R}^{d}}(L_{g}Z_{Rh}f)(a^{\prime})-\sup_{b\in\mathbb{R}^{d}}(Z_{Rh}L_{g}f)(b)\Big{|}$
$\displaystyle\leq\|L_{g}Z_{Rh}f-Z_{Rh}L_{g}f\|_{L^{\infty}(\mathbb{R}^{d})}$
$\displaystyle\leq\|L_{g}Z_{Rh}-Z_{Rh}L_{g}\|_{L^{2}(\mathbb{R}^{d})\to
L^{\infty}(\mathbb{R}^{d})}\cdot\|f\|_{L^{2}(\mathbb{R}^{d})},$

and so the result follows from Lemma 23. ∎

The rest of this section proves Lemma 23. Our proof follows some of the main
ideas in the proof of Lemma E.1 in [52].

###### Proof of Lemma 23.

Denote $K:=Z_{h}-L_{g}Z_{h}L_{g}^{-1}$. Then $L_{g}Z_{h}-Z_{h}L_{g}=-KL_{g}$,
and so

$\|L_{g}Z_{h}-Z_{h}L_{g}\|_{L^{2}(\mathbb{R}^{d})\to
L^{\infty}(\mathbb{R}^{d})}=\|KL_{g}\|_{L^{2}(\mathbb{R}^{d})\to
L^{\infty}(\mathbb{R}^{d})}\leq\|K\|_{L^{2}(\mathbb{R}^{d})\to
L^{\infty}(\mathbb{R}^{d})}\|L_{g}\|_{L^{2}(\mathbb{R}^{d})\to
L^{2}(\mathbb{R}^{d})}.$

We first bound the second factor. For $f\in L^{2}(\mathbb{R}^{d})$, a change
of variables gives

$\|L_{g}f\|_{L^{2}(\mathbb{R}^{d})}^{2}=\int_{\mathbb{R}^{d}}f(g^{-1}(x))^{2}dx=\int_{\mathbb{R}^{d}}f(u)^{2}|\operatorname{det}(Jg(u))|du\leq\sup_{u\in\mathbb{R}^{d}}|\operatorname{det}(Jg(u))|\cdot\|f\|_{L^{2}(\mathbb{R}^{d})}^{2}.$

For $x\in\mathbb{R}^{d}$, the fact that $\|J\tau(x)\|_{2\to 2}\leq\frac{1}{2}$
implies

$\operatorname{det}(Jg^{-1}(x))=\operatorname{det}(I_{d}-J\tau(x))\geq(1-\|J\tau(x)\|_{2\to
2})^{d}\geq 2^{-d},$

and so combining with the above estimate gives

$\|L_{g}\|_{L^{2}(\mathbb{R}^{d})\to
L^{2}(\mathbb{R}^{d})}\leq\sup_{u\in\mathbb{R}^{d}}|\operatorname{det}(Jg(u))|=\sup_{x\in\mathbb{R}^{d}}|\operatorname{det}(Jg^{-1}(x))|^{-1}\leq
2^{d}.$ (5)

It remains to bound $\|K\|_{L^{2}(\mathbb{R}^{d})\to
L^{\infty}(\mathbb{R}^{d})}$. To this end, one may verify that $K$ can be
expressed as $Kf(x)=\int_{\mathbb{R}^{d}}k(x,u)f(u)du$, where the kernel $k$
is defined by

$k(x,u):=h(x-u)-\operatorname{det}(I_{d}-J\tau(u))\cdot
h(x-u-\tau(x)+\tau(u)).$

We will bound the $L^{2}$ norms of every $k(x,\cdot)$ and $k(\cdot,u)$, and
then appeal to Young’s inequality for integral operators to bound
$\|K\|_{L^{2}(\mathbb{R}^{d})\to L^{\infty}(\mathbb{R}^{d})}$. We decompose
$k=k_{1}+k_{2}+k_{3}$, where

$\displaystyle k_{1}(x,u)$
$\displaystyle:=h(x-u)-h\big{(}(I_{d}-J\tau(u))(x-u)\big{)},$ $\displaystyle
k_{2}(x,u)$
$\displaystyle:=\big{(}1-\operatorname{det}(I_{d}-J\tau(u))\big{)}\cdot
h\big{(}(I_{d}-J\tau(u))(x-u)\big{)},$ $\displaystyle k_{3}(x,u)$
$\displaystyle:=\operatorname{det}(I_{d}-J\tau(u))\cdot\Big{(}h\big{(}(I_{d}-J\tau(u))(x-u)\big{)}-h(x-u-\tau(x)+\tau(u))\Big{)}.$

First, we analyze $k_{1}$. Letting $p_{1}\colon[0,1]\to\mathbb{R}^{d}$ denote
the parameterized line segment of constant velocity from
$(I_{d}-J\tau(u))(x-u)$ to $x-u$, we have

$k_{1}(x,u)=\int_{0}^{1}\nabla h(p_{1}(t))\cdot J\tau(u)(x-u)~{}dt,$

and so

$\displaystyle|k_{1}(x,u)|$ $\displaystyle\leq\sup_{t\in[0,1]}\|\nabla
h(p_{1}(t))\|_{2}\cdot\|J\tau(u)(x-u)\|_{2}$
$\displaystyle\leq\sup_{t\in[0,1]}\|\nabla
h(p_{1}(t))\|_{2}\cdot\sup_{z\in\mathbb{R}^{d}}\|J\tau(z)\|_{2\to
2}\cdot\|x-u\|_{2}.$ (6)

To bound the first factor, let $C_{\infty}(h)>0$ denote a simultaneous bound
on the absolute value and $2$-norm of (4). To use this, we bound
$\inf_{t\in[0,1]}\|p_{1}(t)\|_{2}$ from below:

$\displaystyle\|p_{1}(t)\|_{2}=\|x-u+tJ\tau(u)(x-u)\|_{2}$
$\displaystyle\geq\|x-u\|_{2}-t\|J\tau(u)(x-u)\|_{2}$
$\displaystyle\geq(1-t\sup_{z\in\mathbb{R}^{d}}\|J\tau(u)\|_{2\to
2})\cdot\|x-u\|_{2}\geq\tfrac{1}{2}\|x-u\|_{2}.$

Then

$\sup_{t\in[0,1]}\|\nabla
h(p_{1}(t))\|_{2}\leq\sup_{t\in[0,1]}\frac{C_{\infty}(h)}{(1+\|p_{1}(t)\|_{2})^{(d+3)/2}}\leq\frac{C_{\infty}(h)}{(1+\frac{1}{2}\|x-u\|_{2})^{(d+3)/2}},$

which allows us to further bound (6):

$|k_{1}(x,u)|\leq\frac{C_{\infty}(h)}{(1+\frac{1}{2}\|x-u\|_{2})^{(d+3)/2}}\cdot\sup_{z\in\mathbb{R}^{d}}\|J\tau(z)\|_{2\to
2}\cdot\|x-u\|_{2}.$ (7)

Next, we analyze $k_{2}$. Since $\|J\tau(x)\|_{2\to 2}\leq\frac{1}{2}$ by
assumption, Bernoulli’s inequality gives

$1-\operatorname{det}(I_{d}-J\tau(x))\leq 1-(1-\|J\tau(x)\|_{2\to 2})^{d}\leq
d\|J\tau(x)\|_{2\to 2}\leq d\sup_{z\in\mathbb{R}^{d}}\|J\tau(z)\|_{2\to 2}.$

Also, the convexity bound $(1+t)^{d}\leq 1+(2^{d}-1)t$ for $t\in[0,1]$ implies

$\displaystyle 1-\operatorname{det}(I_{d}-J\tau(x))$ $\displaystyle\geq
1-\|I_{d}-J\tau(x)\|_{2\to 2}^{d}$ $\displaystyle\geq 1-(1+\|J\tau(x)\|_{2\to
2})^{d}\geq-2^{d}\sup_{z\in\mathbb{R}^{d}}\|J\tau(z)\|_{2\to 2}.$

Furthermore, we have

$\|(I_{d}-J\tau(x))(x-u)\|_{2}\geq\|x-u\|_{2}-\|J\tau(u)\|_{2\to
2}\|x-u\|_{2}\geq\tfrac{1}{2}\|x-u\|_{2},$

and so

$\displaystyle|k_{2}(x,u)|$
$\displaystyle\leq\big{|}1-\operatorname{det}(I_{d}-J\tau(u))\big{|}\cdot\big{|}h\big{(}(I_{d}-J\tau(u))(x-u)\big{)}\big{|}$
$\displaystyle\leq 2^{d}\sup_{z\in\mathbb{R}^{d}}\|J\tau(z)\|_{2\to
2}\cdot\frac{C_{\infty}(h)}{(1+\frac{1}{2}\|x-u\|_{2})^{(d+1)/2}}.$ (8)

Finally, we analyze $k_{3}$. Put

$r:=x-u-\tau(x)+\tau(u),\qquad s:=\tau(x)-\tau(u)-J\tau(u)(x-u).$

Then letting $p_{2}\colon[0,1]\to\mathbb{R}^{d}$ denote the parameterized line
segment of constant velocity from $r$ to $r+s$, we have

$k_{3}(x,u)=\operatorname{det}(I_{d}-J\tau(u))\cdot\big{(}h(r+s)-h(r)\big{)}=\operatorname{det}(I_{d}-J\tau(u))\int_{0}^{1}\nabla
h(p_{2}(t))\cdot s~{}dt,$

and so

$|k_{3}(x,u)|\leq|\operatorname{det}(I_{d}-J\tau(u))|\cdot\sup_{t\in[0,1]}\|\nabla
h(p_{2}(t))\|_{2}\cdot\|s\|_{2}.$ (9)

For the first factor of (9), we have
$|\operatorname{det}(I_{d}-J\tau(u))|\leq\|I_{d}-J\tau(u)\|_{2\to
2}^{d}\leq(3/2)^{d}$. To bound the second factor of (9), we use our bound
$C_{\infty}(h)>0$ on (4). To do so, we bound
$\inf_{t\in[0,1]}\|p_{2}(t)\|_{2}$ from below. First, we note that

$\tau(x)-\tau(u)=\int_{0}^{1}J\tau(p_{3}(t))(x-u)dt,$

where $p_{3}\colon[0,1]\to\mathbb{R}^{d}$ is the parameterized line segment of
constant velocity from $u$ to $x$. Thus,

$\|\tau(x)-\tau(u)\|_{2}\leq\int_{0}^{1}\|J\tau(p_{3}(t))\|_{2\to
2}\|x-u\|_{2}dt\leq\sup_{z\in\mathbb{R}^{d}}\|J\tau(z)\|_{2\to
2}\cdot\|x-u\|_{2},$ (10)

and so for each $t\in[0,1]$, we have

$\displaystyle\|p_{2}(t)\|_{2}=\|r+ts\|_{2}$
$\displaystyle=\|(x-u)-(1-t)(\tau(x)-\tau(u))-tJ\tau(u)(x-u)\|_{2}$
$\displaystyle\geq\|x-u\|_{2}-(1-t)\|\tau(x)-\tau(u)\|-t\|J\tau(u)\|_{2\to
2}\|x-u\|_{2}$
$\displaystyle\geq(1-\sup_{z\in\mathbb{R}^{d}}\|J\tau(z)\|_{2\to
2})\cdot\|x-u\|_{2}\geq\tfrac{1}{2}\|x-u\|_{2}.$

Overall, we have

$\sup_{t\in[0,1]}\|\nabla
h(p_{2}(t))\|_{2}\leq\sup_{t\in[0,1]}\frac{C_{\infty}(h)}{(1+\|p_{2}(t)\|_{2})^{(d+3)/2}}\leq\frac{C_{\infty}(h)}{(1+\frac{1}{2}\|x-u\|_{2})^{(d+3)/2}}.$

Finally, we apply (10) to bound the third factor of (9):

$\|s\|_{2}\leq\|\tau(x)-\tau(u)\|_{2}+\|J\tau(u)(x-u)\|_{2}\leq
2\cdot\sup_{z\in\mathbb{R}^{d}}\|J\tau(z)\|_{2\to 2}\cdot\|x-u\|_{2}.$

We combine these estimates to obtain the following bound on (9):

$|k_{3}(x,u)|\leq(3/2)^{d}\cdot\frac{C_{\infty}(h)}{(1+\frac{1}{2}\|x-u\|_{2})^{(d+3)/2}}\cdot
2\cdot\sup_{z\in\mathbb{R}^{d}}\|J\tau(z)\|_{2\to 2}\cdot\|x-u\|_{2}.$ (11)

Finally, (7), (8), and (11) together imply

$\displaystyle|k(x,u)|$
$\displaystyle\leq|k_{1}(x,u)|+|k_{2}(x,u)|+|k_{3}(x,u)|$ $\displaystyle\leq
C_{\infty}(h)\cdot\sup_{z\in\mathbb{R}^{d}}\|J\tau(z)\|_{2\to
2}\cdot\bigg{(}\frac{2^{d}}{(1+\frac{1}{2}\|x-u\|_{2})^{(d+1)/2}}+\frac{(2(\frac{3}{2})^{d}+1)\|x-u\|_{2}}{(1+\frac{1}{2}\|x-u\|_{2})^{(d+3)/2}}\bigg{)}.$

Importantly, this is a bounded function of $x-u$ that decays like
$\|x-u\|_{2}^{-(d+1)/2}$. By integrating the square, this simultaneously
bounds the $L^{2}$ norm of every $k(x,\cdot)$ and $k(\cdot,u)$ by a quantity
of the form $C_{0}(h)\cdot\sup_{z\in\mathbb{R}^{d}}\|J\tau(z)\|_{2\to 2}$. By
Young’s inequality for integral operators (see Theorem 0.3.1 in [67], for
example), it follows that

$\|K\|_{L^{2}(\mathbb{R}^{d})\to L^{\infty}(\mathbb{R}^{d})}\leq
C_{0}(h)\cdot\sup_{z\in\mathbb{R}^{d}}\|J\tau(z)\|_{2\to 2}.$

Combining with (5) then gives the result with $C(h):=2^{d}\cdot C_{0}(h)$. ∎

##  5 Template selection for classification

###  5.1 Classifying characteristic functions

In this subsection, we focus on the case in which $V=L^{2}(\mathbb{R}^{d})$
and $G$ is the group of translation operators. Suppose we have $k$ distinct
$G$-orbits of indicator functions of compact subsets of $\mathbb{R}^{d}$.
According to the following result, there is a simple classifier based on a
size-$k$ max filter bank that correctly classifies these orbits. (This cartoon
setting enjoys precursors in [24, 23].)

######  Theorem 24.

Given compact sets $S_{1},\ldots,S_{k}\subseteq\mathbb{R}^{d}$ of positive
measure satisfying

$[\mathbf{1}_{S_{i}}]\neq[\mathbf{1}_{S_{j}}]\qquad\text{whenever}\qquad i\neq
j,$

there exist templates $z_{1},\ldots,z_{k}\in L^{2}(\mathbb{R}^{d})$ satisfying

$\langle\langle[z_{i}],[\mathbf{1}_{S_{j}}]\rangle\rangle<\langle\langle[z_{i}],[\mathbf{1}_{S_{i}}]\rangle\rangle\qquad\text{whenever}\qquad
i\neq j.$

###### Proof.

By compactness, there exists $r>0$ such that every $S_{i}$ is contained in the
closed ball centered at the origin with radius $r$. For reasons that will
become apparent later, we take $B$ to be the closed ball centered at the
origin with radius $3r$, and we define
$z_{i}:=\mathbf{1}_{S_{i}}-\mathbf{1}_{B\setminus S_{i}}$. Then for every $i$
and $j$, it holds that

$\langle\langle[z_{i}],[\mathbf{1}_{S_{j}}]\rangle\rangle=\sup_{a\in\mathbb{R}^{d}}\langle\mathbf{1}_{S_{i}}-\mathbf{1}_{B\setminus
S_{i}},T_{a}\mathbf{1}_{S_{j}}\rangle\leq\sup_{a\in\mathbb{R}^{d}}\langle\mathbf{1}_{S_{i}},T_{a}\mathbf{1}_{S_{j}}\rangle=\langle\langle[\mathbf{1}_{S_{i}}],[\mathbf{1}_{S_{j}}]\rangle\rangle.$
(12)

In the special case where $j=i$, this implies

$|S_{i}|=\langle
z_{i},\mathbf{1}_{S_{i}}\rangle\leq\langle\langle[z_{i}],[\mathbf{1}_{S_{i}}]\rangle\rangle\leq\langle\langle[\mathbf{1}_{S_{i}}],[\mathbf{1}_{S_{i}}]\rangle\rangle=|S_{i}|,$

and so $\langle\langle[z_{i}],[\mathbf{1}_{S_{i}}]\rangle\rangle=|S_{i}|$. We
consider all $j\neq i$ in two cases.

Case I: $j\neq i$ and $|S_{j}|\leq|S_{i}|$. Considering (12), it suffices to
bound $\langle\langle[\mathbf{1}_{S_{i}}],[\mathbf{1}_{S_{j}}]\rangle\rangle$.
Letting $R$ denote the reversal operator defined by $Rf(x):=f(-x)$, then

$\langle\langle[\mathbf{1}_{S_{i}}],[\mathbf{1}_{S_{j}}]\rangle\rangle=\sup_{a\in\mathbb{R}^{d}}(R\mathbf{1}_{S_{i}}\star\mathbf{1}_{S_{j}})(a).$

Since $R\mathbf{1}_{S_{i}},\mathbf{1}_{S_{j}}\in L^{2}(\mathbb{R}^{d})$, it
holds that the convolution $R\mathbf{1}_{S_{i}}\star\mathbf{1}_{S_{j}}$ is
continuous, and since $S_{i}$ and $S_{j}$ are compact, the convolution has
compact support. Thus, the extreme value theorem gives that the convolution
achieves its supremum, meaning there exists $a\in\mathbb{R}^{d}$ such that

$\langle\langle[\mathbf{1}_{S_{i}}],[\mathbf{1}_{S_{j}}]\rangle\rangle=\langle\mathbf{1}_{S_{i}},T_{a}\mathbf{1}_{S_{j}}\rangle=|S_{i}\cap(S_{j}+a)|.$
(13)

Next, the assumptions $|S_{j}|\leq|S_{i}|$ and
$[\mathbf{1}_{S_{i}}]\neq[\mathbf{1}_{S_{j}}]$ together imply

$|S_{i}\cap(S_{j}+a)|<|S_{i}|.$ (14)

Indeed, equality in the bound $|S_{i}\cap(S_{j}+a)|\leq|S_{i}|$ is only
possible if $S_{i}\subseteq S_{j}+a$ (modulo null sets), but since
$|S_{j}|\leq|S_{i}|$ by assumption, this requires $S_{i}=S_{j}+a$ (modulo null
sets), which violates the assumption
$[\mathbf{1}_{S_{i}}]\neq[\mathbf{1}_{S_{j}}]$. Overall, we combine (12),
(13), and (14) to get

$\langle\langle[z_{i}],[\mathbf{1}_{S_{j}}]\rangle\rangle\leq\langle\langle[\mathbf{1}_{S_{i}}],[\mathbf{1}_{S_{j}}]\rangle\rangle=|S_{i}\cap(S_{j}+a)|<|S_{i}|=\langle\langle[z_{i}],[\mathbf{1}_{S_{i}}]\rangle\rangle.$

Case II: $j\neq i$ and $|S_{j}|\geq|S_{i}|$. If
$\langle\langle[z_{i}],[\mathbf{1}_{S_{j}}]\rangle\rangle\leq 0$, then

$\langle\langle[z_{i}],[\mathbf{1}_{S_{j}}]\rangle\rangle\leq
0<|S_{i}|=\langle\langle[z_{i}],[\mathbf{1}_{S_{i}}]\rangle\rangle,$

and so we are done. Now suppose
$\langle\langle[z_{i}],[\mathbf{1}_{S_{j}}]\rangle\rangle>0$. Considering

$\langle\langle[z_{i}],[\mathbf{1}_{S_{j}}]\rangle\rangle=\sup_{a\in\mathbb{R}^{d}}(Rz_{i}\star\mathbf{1}_{S_{j}})(a),$

then by continuity and compactness, the extreme value theorem produces
$a\in\mathbb{R}^{d}$ such that

$\langle\langle[z_{i}],[\mathbf{1}_{S_{j}}]\rangle\rangle=\langle
z_{i},T_{a}\mathbf{1}_{S_{j}}\rangle=|S_{i}\cap(S_{j}+a)|-|(B\setminus
S_{i})\cap(S_{j}+a)|.$ (15)

Since $\langle\langle[z_{i}],[\mathbf{1}_{S_{j}}]\rangle\rangle>0$, it follows
that $|S_{i}\cap(S_{j}+a)|>0$, i.e., $S_{i}\cap(S_{j}+a)$ is nonempty, which
in turn implies $S_{j}+a\subseteq B$. (This is why we defined $B$ to have
radius $3r$.) As before, $|S_{j}|\geq|S_{i}|$ and
$[\mathbf{1}_{S_{i}}]\neq[\mathbf{1}_{S_{j}}]$ together give
$|S_{i}\cap(S_{j}+a)|<|S_{j}|$. Thus,

$|(B\setminus
S_{i})\cap(S_{j}+a)|=|B\cap(S_{j}+a)|-|S_{i}\cap(S_{j}+a)|=|S_{j}+a|-|S_{i}\cap(S_{j}+a)|>0.$
(16)

We combine (15) and (16) to get

$\langle\langle[z_{i}],[\mathbf{1}_{S_{j}}]\rangle\rangle=|S_{i}\cap(S_{j}+a)|-|(B\setminus
S_{i})\cap(S_{j}+a)|<|S_{i}\cap(S_{j}+a)|\leq|S_{i}|=\langle\langle[z_{i}],[\mathbf{1}_{S_{i}}]\rangle\rangle,$

as claimed. ∎

For each $i$, assume $S_{i}$ is translated so that it is contained in the
smallest possible ball centered at the origin, and let $r_{i}$ denote the
radius of this ball. The proof of Theorem 24 gives that each template $z_{i}$
is supported in a closed ball of radius $R:=3\max_{i}r_{i}$. The fact that
these templates are localized bears some consequence for certain image
articulation manifolds [35]. In particular, for each
$\pi\colon\\{1,\ldots,k\\}\to\mathbb{N}\cup\\{0\\}$, let $M_{\pi}\subseteq
L^{2}(\mathbb{R}^{d})$ denote the manifold of images of the form

$\sum_{i=1}^{k}\sum_{j=1}^{\pi(i)}T_{a(i,j)}\mathbf{1}_{S_{i}},\qquad\text{where}\qquad\|a(i,j)-a(i^{\prime},j^{\prime})\|>4R\quad\forall(i,j)\neq(i^{\prime},j^{\prime}).$

Thanks to the $4R$ spacing, each translate of each template interacts with at
most one component of the image, and so for every $f\in M_{\pi}$, it holds
that

$\langle\langle[z_{i}],[f]\rangle\rangle=\max_{i^{\prime}:\pi(i^{\prime})>0}\langle\langle[z_{i}],[\mathbf{1}_{S_{i^{\prime}}}]\rangle\rangle.$

In particular, the same max filter bank can be used to determine the support
of $\pi$. As an example, if some multiset of characters are typed on a page in
a common font and with sufficient separation, then the max filter bank from
Theorem 24 that distinguishes the characters can be used to determine which
ones appear on the page.

###  5.2 Classifying mixtures of stationary processes

In this subsection, we focus on the case in which $V=\mathbb{R}^{n}$ and
$G\cong C_{n}$ is the group of circular translation operators. A natural
$G$-invariant probability distribution is a multivariate Gaussian with mean
zero and circulant covariance, and so we consider the task of classifying a
mixture of such distributions. One-dimensional textures can be modeled in this
way, especially if the covariance matrix has a small bandwidth so that distant
pixels are statistically independent. A standard approach for this problem is
to estimate the first- and second-order moments given a random draw. As we
will soon see, one can alternatively classify with high accuracy by
thresholding a single max filter. In what follows, we make use of Thompson’s
part metric on the set of positive definite matrices:

$d_{\infty}(A,B):=\|\log(A^{-1/2}BA^{-1/2})\|_{2\to 2}.$

We also let $A_{k}$ denote the leading $k\times k$ principal submatrix of $A$.

######  Theorem 25.

Fix $C>\log 2$, take $n,w\in\mathbb{N}$ such that
$k:=\lfloor\sqrt{n/2}\rfloor\geq w$, and consider any positive definite
$A,B\in\mathbb{R}^{n\times n}$ that are circulant with bandwidth $w$ and
satisfy

$d_{\infty}(A_{k},B_{k})\geq C.$ (17)

There exists $z\in\mathbb{R}^{n}$ supported on an interval of length $k$ and a
threshold $\theta\in\mathbb{R}$ such that for every mixture $\mathsf{M}$ of
$\mathsf{N}(0,A)$ and $\mathsf{N}(0,B)$, then given $x\sim\mathsf{M}$, the
comparison

$\langle\langle[z],[x]\rangle\rangle\gtrless\theta$

correctly classifies the latent mixture component of $x$ with probability
$1-o_{n\to\infty;C}(1)$.

###### Proof.

First, we observe that $d_{\infty}(A_{k},B_{k})\geq C$ is equivalent to

$\max\Big{\\{}\lambda_{\mathrm{max}}(A_{k}^{-1/2}B_{k}A_{k}^{-1/2}),\lambda_{\mathrm{max}}(B_{k}^{-1/2}A_{k}B_{k}^{-1/2})\Big{\\}}\geq
e^{C}.$

Without loss of generality, we may assume
$\lambda_{\mathrm{max}}(A_{k}^{-1/2}B_{k}A_{k}^{-1/2})\geq e^{C}$. Let
$v\in\mathbb{R}^{k}$ denote a corresponding unit eigenvector of
$A_{k}^{-1/2}B_{k}A_{k}^{-1/2}$, and define $z\in\mathbb{R}^{n}$ to be
supported in its first $k$ entries as the subvector $z_{k}:=A_{k}^{-1/2}v$.
Then

$\frac{z^{\top}Bz}{z^{\top}Az}=\frac{z_{k}^{\top}B_{k}z_{k}}{z_{k}^{\top}A_{k}z_{k}}=\frac{v^{\top}A_{k}^{-1/2}B_{k}A_{k}^{-1/2}v}{v^{\top}v}=\lambda_{\mathrm{max}}(A_{k}^{-1/2}B_{k}A_{k}^{-1/2})\geq
e^{C}.$ (18)

Consider $x_{1}\sim\mathsf{N}(0,A)$. Then
$\langle\langle[z],[x]\rangle\rangle=\max_{a\in C_{n}}\langle
T_{a}z,x_{1}\rangle$, where each $\langle T_{a}z,x_{1}\rangle$ has Gaussian
distribution with mean zero and variance $(T_{a}z)^{\top}A(T_{a}z)$, which in
turn equals $z^{\top}Az$ since $A$ is circulant. Denoting
$Z\sim\mathsf{N}(0,1)$, a union bound then gives

$\mathbb{P}\big{\\{}\langle\langle[z],[x_{1}]\rangle\rangle\geq t\big{\\}}\leq
n\cdot\mathbb{P}\big{\\{}(z^{\top}Az)^{1/2}\cdot Z\geq t\big{\\}}\leq
ne^{-t^{2}/(2z^{\top}Az)}$

for $t\geq 0$. This failure probability is $o_{n\to\infty;c_{1}}(1)$ by taking
$t:=\sqrt{c_{1}\cdot z^{\top}Az\cdot\log n}$ for any $c_{1}>2$. Next, consider
$x_{2}\sim\mathsf{N}(0,B)$, and take any subset $S\subseteq C_{n}$ consisting
of $k$ members of $C_{n}$ of pairwise distance at least $2k$. Then

$\langle\langle[z],[x_{2}]\rangle\rangle=\max_{a\in C_{n}}\langle
T_{a}z,x_{2}\rangle\geq\max_{a\in S}\langle T_{a}z,x_{2}\rangle.$

In this case, each $\langle T_{a}z,x_{2}\rangle$ has Gaussian distribution
with mean zero and variance $z^{\top}Bz$. Furthermore, since $k\geq w$, these
random variables have pairwise covariance zero, and since they are Gaussian,
they are therefore independent. For independent
$Z_{1},\ldots,Z_{k}\sim\mathsf{N}(0,1)$ and $t\geq 0$, we have

$\mathbb{P}\Big{\\{}\max_{i\in\\{1,\ldots,k\\}}Z_{i}\leq
t\Big{\\}}=\mathbb{P}\\{Z\leq
t\\}^{k}\leq(1-\tfrac{t}{1+t^{2}}\cdot\tfrac{1}{\sqrt{2\pi}}e^{-t^{2}/2})^{k},$

where the last step follows from a standard lower bound on the tail of a
standard Gaussian distribution. Take $t:=\sqrt{c_{2}\log k}$ to get

$\mathbb{P}\Big{\\{}\max_{i\in\\{1,\ldots,k\\}}Z_{i}\leq\sqrt{c_{2}\log
k}\Big{\\}}\leq\Big{(}1-(\tfrac{\sqrt{c_{2}\log k}}{1+c_{2}\log
k}\cdot\tfrac{1}{\sqrt{2\pi}}\cdot
k^{1-c_{2}/2})\cdot\tfrac{1}{k}\Big{)}^{k}\leq\exp\Big{(}-\tfrac{\sqrt{c_{2}\log
k}}{1+c_{2}\log k}\cdot\tfrac{1}{\sqrt{2\pi}}\cdot k^{1-c_{2}/2}\Big{)},$

which is $o_{k\to\infty;c_{2}}(1)$ when $c_{2}<2$. Overall, we simultaneously
have

$\langle\langle[z],[x_{1}]\rangle\rangle<\theta_{1}:=\sqrt{c_{1}\cdot
z^{\top}Az\cdot\log
n},\qquad\langle\langle[z],[x_{2}]\rangle\rangle>\theta_{2}:=\sqrt{\tfrac{1}{2}c_{2}\cdot
z^{\top}Bz\cdot\log n}$

with probability $1-o_{n\to\infty;c_{1},c_{2}}(1)$, provided $c_{1}>2>c_{2}$.
We select $c_{1}:=2(\tfrac{e^{C}}{2})^{1/2}>2$ and
$c_{2}:=2(\tfrac{e^{C}}{2})^{-1/2}<2$, and then (18) implies

$\frac{\theta_{2}^{2}}{\theta_{1}^{2}}=\frac{\frac{1}{2}c_{2}\cdot
z^{\top}Bz}{c_{1}\cdot
z^{\top}Az}=\frac{1}{e^{C}}\cdot\frac{z^{\top}Bz}{z^{\top}Az}\geq 1.$

Thus, $\theta_{1}\leq\theta_{2}$, and so the result follows by taking
$\theta:=(\theta_{1}+\theta_{2})/2$. ∎

Given a mixture of $k$ Gaussians with covariance matrices that pairwise
satisfy (17), then we can perform multiclass classification by a one-vs-one
reduction. Indeed, the binary classifier in Theorem 25 can be applied to all
$\binom{k}{2}$ pairs of Gaussians, in which case we correctly classify the
latent mixture component with high probability (provided $k$ is fixed and
$n\to\infty$).

Interestingly, max filters can also distinguish between stationary processes
with identical first- and second-order moments. For example,
$x\sim\mathsf{Unif}(\\{\pm 1\\}^{n})$ and $y\sim\mathsf{N}(0,I_{n})$ are both
stationary with mean zero and identity covariance. If $z\in\mathbb{R}^{n}$ is
a standard basis element, then with high probability, it holds that

$\langle\langle[z],[x]\rangle\rangle\leq 1\leq\sqrt{\log
n}\leq\langle\langle[z],[y]\rangle\rangle.$

This indicates that max filters incorporate higher-order moments.

###  5.3 Subgradients

In this subsection, we focus on the case in which $V=\mathbb{R}^{d}$ and $G$
is a closed subgroup of $\operatorname{O}(d)$. The previous subsections
carefully designed templates to classify certain data models. For real-world
classification problems, one is expected to train a classifier on a given
training set of labeled data. To do so, one selects a parameterized family of
classifiers and then locally minimizes some notion of training loss over this
family. This is feasible provided the classifier is a differentiable function
of the parameters. We envision a classifier in which the first layer is a max
filter bank, and so this subsection establishes how to differentiate a max
filter with respect to the template. This is made possible by convexity; see
Lemma 2(d).

Every convex function $f\colon\mathbb{R}^{d}\to\mathbb{R}$ has a
subdifferential $\partial f\colon\mathbb{R}^{d}\to 2^{\mathbb{R}^{d}}$ defined
by

$\partial f(x):=\\{u\in\mathbb{R}^{d}:f(x+h)\geq f(x)+\langle
h,u\rangle~{}~{}\forall h\in\mathbb{R}^{d}\\}.$

For a fixed $x\in\mathbb{R}^{d}$ and $u\in\partial f(x)$, it is helpful to
interpret the graph of $z\mapsto f(x)+\langle z-x,u\rangle$ as a supporting
hyperplane of the epigraph of $f$ at $x$. For example, the absolute value
function over $\mathbb{R}$ has the following subdifferential:

$\partial|\cdot|(x)=\left\\{\begin{array}[]{cl}\\{-1\\}&\text{if }x<0,\\\
{[}-1,1]&\text{if }x=0,\\\ \\{1\\}&\text{if }x>0.\end{array}\right.$

Indeed, this gives the slopes of the hyperplanes that support the epigraph of
$|\cdot|$ at $x\in\mathbb{R}$.

Following [74], we will determine the subdifferential of the max filtering map
by first finding its directional derivatives. In general, the directional
derivative of $f$ at $x$ in the direction of $v\neq 0$ is given by

$f^{\prime}(x;v):=\lim_{t\to 0^{+}}\frac{f(x+tv)-f(x)}{t}=\sup_{u\in\partial
f(x)}\langle u,v\rangle.$

The second equality is a standard result; see for example Theorem 23.4 in
[61]. It will be convenient to denote the set

$G(x,y):=\arg\max_{g\in G}\langle x,gy\rangle.$

Observe that $G(x,y)$ is a closed subset of $G$ since the map $g\mapsto\langle
x,gy\rangle$ is continuous.

######  Lemma 26.

Suppose $G$ is a closed subgroup of $\operatorname{O}(d)$ and select
$x,y,v\in\mathbb{R}^{d}$ with $v\neq 0$. Then

$\langle\langle[\cdot],[y]\rangle\rangle^{\prime}(x;v)=\max_{g\in
G(x,y)}\langle v,gy\rangle.$

###### Proof.

For each $t>0$ and $g\in G(x,y)$, we have

$\langle\langle[x+tv],[y]\rangle\rangle\geq\langle x+tv,gy\rangle=\langle
x,gy\rangle+t\langle v,gy\rangle=\langle\langle[x],[y]\rangle\rangle+t\langle
v,gy\rangle.$ (19)

For each $t>0$, select $g_{t}\in G(x+tv,y)$. Then

$\langle\langle[x+tv],[y]\rangle\rangle=\langle x+tv,g_{t}y\rangle=\langle
x,g_{t}y\rangle+t\langle
v,g_{t}y\rangle\leq\langle\langle[x],[y]\rangle\rangle+t\langle
v,g_{t}y\rangle.$ (20)

We rearrange and combine (19) and (20) to get

$\max_{g\in G(x,y)}\langle
v,gy\rangle\leq\frac{\langle\langle[x+tv],[y]\rangle\rangle-\langle\langle[x],[y]\rangle\rangle}{t}\leq\langle
v,g_{t}y\rangle$ (21)

for all $t>0$. Select a sequence $t_{n}\to 0^{+}$ such that $g_{t_{n}}$
converges to some $g^{\star}\in G$ (by passing to a subsequence if necessary).
Then continuity implies

$\langle
x+t_{n}v,g_{t_{n}}y\rangle=\langle\langle[x+t_{n}v],[y]\rangle\rangle\to\langle\langle[x],[y]\rangle\rangle.$

Furthermore,

$\displaystyle|\langle x+t_{n}v,g_{t_{n}}y\rangle-\langle
x,g^{\star}y\rangle|$ $\displaystyle\leq|\langle x,g_{t_{n}}y\rangle-\langle
x,g^{\star}y\rangle|+|\langle t_{n}v,g_{t_{n}}y\rangle|$
$\displaystyle\leq\|x\|\|g_{t_{n}}-g^{\star}\|_{2\to
2}\|y\|+t_{n}\|v\|\|y\|\to 0.$

It follows that $\langle
x,g^{\star}y\rangle=\langle\langle[x],[y]\rangle\rangle$, i.e., $g^{\star}\in
G(x,y)$. Taking limits of (21) then gives

$\max_{g\in G(x,y)}\langle v,gy\rangle\leq\lim_{t\to
0^{+}}\frac{\langle\langle[x+tv],[y]\rangle\rangle-\langle\langle[x],[y]\rangle\rangle}{t}\leq\lim_{t\to
0^{+}}\langle v,g_{t}y\rangle=\langle v,g^{\star}y\rangle\leq\max_{g\in
G(x,y)}\langle v,gy\rangle.\qed$

######  Theorem 27.

Suppose $G$ is a closed subgroup of $\operatorname{O}(d)$ and select
$x,y\in\mathbb{R}^{d}$. Then

$\partial\langle\langle[\cdot],[y]\rangle\rangle(x)=\operatorname{conv}\\{gy:g\in
G(x,y)\\}.$

###### Proof.

First, we claim that for every $g\in G(x,y)$, the vector $gy$ is in the
subdifferential $\partial\langle\langle[\cdot],[y]\rangle\rangle(x)$. Indeed,
for every $h\in\mathbb{R}^{d}$, we have

$\langle\langle[x+h],[y]\rangle\rangle\geq\langle x+h,gy\rangle=\langle
x,gy\rangle+\langle h,gy\rangle=\langle\langle[x],[y]\rangle\rangle+\langle
h,gy\rangle,$

where the last step uses the fact that $g\in G(x,y)$. Since
$\partial\langle\langle[\cdot],[y]\rangle\rangle(x)$ is convex, it follows
that the $\supseteq$ portion of the desired result holds. Next, suppose there
exists $w\in\partial\langle\langle[\cdot],[y]\rangle\rangle(x)$ such that
$w\not\in\operatorname{conv}\\{gy:g\in G(x,y)\\}$. Then there exists a
hyperplane that separates $w$ from $\\{gy:g\in G(x,y)\\}$, i.e., there is a
nonzero vector $v\in\mathbb{R}^{d}$ such that

$\langle v,gy\rangle<\langle v,w\rangle$

for every $g\in G(x,y)$. By Lemma 26, it follows that

$\langle\langle[\cdot],[y]\rangle\rangle^{\prime}(x;v)=\max_{g\in
G(x,y)}\langle v,gy\rangle<\langle
w,v\rangle\leq\sup_{u\in\partial\langle\langle[\cdot],[y]\rangle\rangle(x)}\langle
u,v\rangle=\langle\langle[\cdot],[y]\rangle\rangle^{\prime}(x;v),$

a contradiction. This establishes the $\subseteq$ portion of the desired
result. ∎

###  5.4 Random templates and limit laws

While the previous subsection was concerned with the differentiability needed
to optimize a feature map, it is well known that random feature maps suffice
for various tasks (e.g., Johnson–Lindenstrauss maps [46] and random kitchen
sinks [60]). In fact, our bilipschitz result (Theorem 18) uses random
templates. As such, one may be inclined to use random templates to produce
features for classification.

In this subsection, we focus on the case in which $V=\mathbb{R}^{d}$ and
$G\cong S_{d}$ is the group of $d\times d$ permutation matrices. Consider a
max filter bank consisting of independent standard gaussian templates
$z_{1},\ldots,z_{n}\in\mathbb{R}^{d}$. Then

$\langle\langle[z_{i}],[x]\rangle\rangle=\langle\operatorname{sort}(z_{i}),\operatorname{sort}(x)\rangle.$

When $d$ is large, we expect the vectors $\operatorname{sort}(z_{i})$ to
exhibit little variation, and so we are inclined to perform dimensionality
reduction. Figure 1 illustrates that the principal components of
$\\{\operatorname{sort}(z_{i})\\}_{i=1}^{n}$ exhibit a high degree of
regularity.

Figure 1:  (left) Draw $n=10,000$ independent standard gaussian random vectors
in $\mathbb{R}^{d}$ with $d=1,000$, sort the coordinates of each vector, and
then plot the top $6$ eigenvectors of the resulting sample covariance matrix.
(right) Discretize and plot the top $6$ eigenfunctions described in Theorem
28.

To explain this regularity, select $Q\colon(0,1)\to\mathbb{R}$ so that
$Q^{-1}$ is the cumulative distribution function of the standard normal
distribution. Take $z\in\mathsf{N}(0,I_{d})$ and put
$s:=\operatorname{sort}(z)$. Denoting $p_{i}:=\frac{i}{d+1}$ and
$q_{i}:=1-p_{i}$, then Section 4.6 in [30] gives

$\mathbb{E}[s_{i}]=Q(p_{i})+O(\tfrac{1}{d}),\qquad\operatorname{Cov}(s_{i},s_{j})=\tfrac{1}{d+2}p_{i}q_{j}Q^{\prime}(p_{i})Q^{\prime}(p_{j})+O(\tfrac{1}{d^{2}}),\qquad
i\leq j.$

The principal components of $\\{\operatorname{sort}(z_{i})\\}_{i=1}^{n}$
approximate the top eigenvectors of the covariance matrix, which in turn
approximate discretizations of eigenfunctions of the integral operator with
kernel $K\colon(0,1)^{2}\to\mathbb{R}$ defined by

$K(x,y):=\min\\{x,y\\}\cdot(1-\max\\{x,y\\})\cdot Q^{\prime}(x)\cdot
Q^{\prime}(y).$

The following result expresses these eigenfunctions in terms of the
probabilist’s Hermite polynomials, which are defined by the following
recurrence:

$p_{0}(x):=1,\qquad p_{n+1}(x):=xp_{n}(x)-p_{n}^{\prime}(x).$ (22)

######  Theorem 28.

The integral operator $L\colon L^{2}([0,1])\to L^{2}([0,1])$ defined by

$Lf(x):=\int_{0}^{1}K(x,y)f(y)dy$

has eigenvalue $\frac{1}{n+1}$ with corresponding eigenfunction $p_{n}\circ Q$
for each $n\in\mathbb{N}\cup\\{0\\}$.

As such, instead of max filtering with independent gaussian templates, one can
efficiently capture the same information by taking the inner product between
$\operatorname{sort}(x)$ and discretized versions of the eigenfunctions
$p_{n}\circ Q$. To reduce the dimensionality, one can simply use fewer
eigenfunctions. To prove Theorem 28, we will use the following lemma; here,
$\varphi\colon\mathbb{R}\to\mathbb{R}$ denotes the probability density
function of the standard normal distribution.

######  Lemma 29.

  * (a)

$p_{n}(x)=\frac{1}{n+1}p_{n+1}^{\prime}(x)$.

  * (b)

$\int_{0}^{x}p_{n+1}(Q(y))dy=-\varphi(Q(x))p_{n}(Q(x))$.

The proof of Lemma 29 follows quickly from the recurrence (22).

###### Proof of Theorem 28.

We compute $L(p_{n}\circ Q)$ by splitting the integral:

$\displaystyle L(p_{n}\circ Q)(x)$
$\displaystyle=(1-x)Q^{\prime}(x)\underbrace{\int_{0}^{x}yQ^{\prime}(y)p_{n}(Q(y))dy}_{I_{1}}+xQ^{\prime}(x)\underbrace{\int_{x}^{1}(1-y)Q^{\prime}(y)p_{n}(Q(y))dy}_{I_{2}}.$

For both $I_{1}$ and $I_{2}$, we integrate by parts with

$dv=Q^{\prime}(y)p_{n}(Q(y))dy=Q^{\prime}(y)\cdot\tfrac{1}{n+1}p_{n+1}^{\prime}(Q(y))\cdot
dy=\tfrac{1}{n+1}(p_{n+1}\circ Q)^{\prime}(y)dy,$

where the middle step follows from Lemma 29(a). Then Lemma 29(b) gives

$\displaystyle I_{1}$
$\displaystyle=yv\Big{|}_{0}^{x}-\int_{0}^{x}vdy=\tfrac{1}{n+1}\Big{(}xp_{n+1}(Q(x))+\varphi(Q(x))p_{n}(Q(x))\Big{)},$
$\displaystyle I_{2}$
$\displaystyle=(1-y)v\Big{|}_{x}^{1}+\int_{x}^{1}vdy=\tfrac{1}{n+1}\Big{(}-(1-x)p_{n+1}(Q(x))+\varphi(Q(x))p_{n}(Q(x))\Big{)}.$

These combine (and mostly cancel) to give

$L(p_{n}\circ
Q)(x)=\tfrac{1}{n+1}Q^{\prime}(x)p_{n}(Q(x))\varphi(Q(x))=\tfrac{1}{n+1}(p_{n}\circ
Q)(x).\qed$

##  6 Numerical examples

In this section, we use max filters as feature maps for various real-world
learning tasks.

######  Example 30 (Voting districts).

The one person, one vote principle insists that in each state, different
voting districts must have nearly the same number of constituents. This
principle is enforced with the help of a decennial redistricting process based
on U.S. Census data. Interestingly, each state assembly applies its own
process for redistricting, and partisan approaches can produce unwanted
gerrymandering. Historically, gerrymandering is detected by how contorted a
district’s shape looks; for example, the Washington Post article [44] uses a
particular geometric score to identify the top $10$ most gerrymandered voting
districts of the $113$th Congress.

As an alternative, we visualize the distribution of district shapes with the
help of max filtering. The shape files for all voting districts of the $116$th
Congress are available in [68]. We center each district at the origin and
scale it to have unit perimeter, and then we sample the boundary at $n=50$
equally spaced points. This results in a $2\times 50$ matrix representation of
each district. However, the same district shape may be represented by many
matrices corresponding to an action of $\operatorname{O}(2)$ on the left and
an action of $C_{n}$ that cyclically permutes columns. Hence, it is
appropriate to apply max filtering with $G\cong\operatorname{O}(2)\times
C_{n}$. It is convenient to identify $\mathbb{R}^{2\times 50}$ with
$\mathbb{C}^{50}$ so that the corresponding max filter is given by

$\displaystyle\langle\langle[z],[x]\rangle\rangle$
$\displaystyle=\max\Big{\\{}~{}\max_{a\in C_{n}}|z^{*}T_{a}x|,~{}\max_{a\in
C_{n}}|\overline{z}^{*}T_{a}x|~{}\Big{\\}}$
$\displaystyle=\max\Big{\\{}~{}\max_{a\in C_{n}}|(R\overline{z}\star
x)(a)|,~{}\max_{a\in C_{n}}|(Rz\star x)(a)|~{}\Big{\\}},$

which can be computed efficiently with the help of the fast Fourier transform.
We max filter with $100$ random templates to embed these districts in the
feature space $\mathbb{R}^{100}$, and then we visualize the result using PCA;
see Figure 2.

Figure 2:  Visualization of voting districts of the $116$th Congress obtained
by max filtering and principal component analysis. See Example 30 for details.

Interestingly, the principal components of this feature domain appear to be
interpretable. The first principal component (given by the horizontal axis)
seems to capture the extent to which the district is convex, while the second
principal component seems to capture the eccentricity of the district. Six of
the ten most gerrymandered districts from [44] are drawn in red, which appear
at relatively extreme points in this feature space. The four remaining
districts (NC-1, NC-4, NC-12, and PA-7) were redrawn by court order between
the $113$th and $116$th Congresses; the new versions of these districts are
drawn in blue. Unsurprisingly, the redrawn districts are not contorted, and
they appear at not-so-extreme points in the feature space.

######  Example 31 (ECG time series).

An electrocardiogram (ECG) uses electrodes placed on the skin to quantify the
heart’s electrical activity over time. With ECG data, a physician can
determine whether a patient has had a heart attack with about 70% accuracy
[51]. Recently, Makimoto et al. [51] trained a $6$-layer convolutional neural
network to perform this task with about 80% accuracy. The dataset they used is
available at [14], which consists of $(12+3)$-lead ECG data sampled at $1$ kHz
from $148$ patients who recently had a heart attack (i.e., myocardial
infarction or MI) and from $141$ patients who did not.

As an alternative to convolutional neural networks, we use max filters to
classify patients based on one second of ECG data (i.e., roughly one
heartbeat). In particular, we read the first $t=1000$ samples of all $15$
leads to obtain a $15\times t$ matrix $X$. We then lift this matrix to a
$15\times w\times(t-w+1)$ tensor with $w=30$, where each $15\times w$ slice
corresponds to a contiguous $15\times w$ submatrix of $X$. Then we normalize
each $1\times w\times 1$ subvector to have mean zero to get the tensor $Y$;
this discards any trend in the data. We account for time invariance by taking
$G$ to be the order-$(t-w+1)$ group of circular permutations of the $15\times
w$ slices of $Y$. Using this group, we max filter with $n=5$ different
templates and then classify with a support-vector machine. We constrain each
template to be supported on a single slice, and we train these templates
together with the support-vector machine classifier by gradient descent to
minimize hinge loss.

Figure 3:  We train max filter templates and a support-vector machine
classifier on electrocardiogram data to distinguish between patients who have
had a heart attack from those who have not. Above, we plot the most extreme
examples in the test set (those with heart attacks on the left, and those
without on the right). The blue windows illustrate the time segments of width
$w=30$ that align best with the templates. See Example 31 for details.

Following [51], we form the test set by randomly drawing $25$ examples from
the MI class and $25$ examples from the non-MI class, and then we form the
training set by randomly drawing $108$ non-test set examples from each class.
We perform $10$ trials of this experiment, and each achieve between $74$% and
$84$% accuracy on the test set. In particular, this is competitive with the
$6$-layer convolutional neural network in [51] despite accessing only a
fraction of the data. In Figure 3, we illustrate what the classifier considers
to be the most extreme examples in the test set. The time segments that align
best with the trained templates typically cover the heartbeat portion of the
ECG signal.

######  Example 32 (Textures).

In this example, we use max filtering to classify various textures,
specifically, those given in the Kylberg Texture Dataset v. 1.0, available in
[49]. This dataset consists of 28 classes of textures, such as blanket, grass,
and oatmeal. Each class consists of 160 distinct grayscale images. (We crop
these $576\times 576$ images down to $256\times 256$ for convenience.) A few
classifiers in the literature have been trained on this dataset [50, 3], but
in order to achieve a high level of performance, they leverage domain-specific
feature maps, augment the training set with other image sets, or simply
initialize with a pre-trained network. As an alternative, we apply max
filtering with random features, and we succeed with much smaller training
sets. For an honest comparison, we train different classifiers on training
sets of different sizes, and the results are displayed in Table 1. (Here, the
test set consists of $32$ points from each class.) We describe each classifier
in what follows.

Table 1:  Accuracy on test set versus size of training set (baseline = 0.035, see Example 32)  pts per class |  CNN |  LDA |  PCA-LDA |  our method
---|---|---|---|---
2 |  0.125 |  0.043 |  0.063 |  0.698
4 |  0.114 |  0.103 |  0.060 |  0.947
8 |  0.150 |  0.158 |  0.128 |  0.970
16 |  0.214 |  0.125 |  0.128 |  0.970
32 |  0.306 |  0.116 |  0.120 |  0.968
64 |  0.534 |  0.117 |  0.122 |  0.964
128 |  0.637 |  (OOM) |  0.128 |  0.962

First, we consider a convolutional neural network (CNN) with a standard
architecture. We pass the $256\times 256$ image through a convolutional layer
with a $3\times 3$ kernel before max pooling over $2\times 2$ patches. Then we
pass the result through another convolutional layer with a $3\times 3$ kernel,
and again max pool over $2\times 2$ patches. We end with a dense layer that
maps to $28$ dimensions, one for each class. We train by minimizing cross
entropy loss against the one-hot vector representation of each label. This
model has $1.7$ million parameters, which enables us to achieve perfect
accuracy on the training set.

As an alternative, we apply linear discriminant analysis (LDA). Training this
classifier involves a massive matrix operation, and this requires too much
memory in our setting when the training set has $128$ points per class. This
motivates the use of principal component analysis to reduce the dimensionality
of each $256\times 256$ image to $k=25$ principal components. The resulting
classifier (PCA-LDA) works for larger training sets, and even exhibits
improved accuracy when the training set is not too small.

For our method, we apply the same PCA-LDA method to a max filtering feature
domain. Surprisingly, we find that modding out by the entire group of pixel
permutations (i.e., just sorting the pixel values) performs reasonably well as
a feature map. Our method implements a multiscale version of this observation.
For each $\ell\in\\{0,1,\ldots,8\\}$, one may partition the $2^{8}\times
2^{8}$ image into $4^{8-\ell}$ square patches of width $2^{\ell}$. We perform
max filtering with the group $(S_{4^{\ell}})^{4^{8-\ell}}$ of all patch-
preserving pixel permutations for each $\ell\in\\{2,\ldots,8\\}$. For
simplicity, we are inclined to draw templates at random, but for computational
efficiency, we simulate random templates with the help of Theorem 28. In
particular, for each $\ell\in\\{2,\ldots,8\\}$ and $n\in\\{0,\ldots,5\\}$, we
sort the pixels in each $2^{\ell}\times 2^{\ell}$ patch and take their inner
products with a discretized version of $p_{n}\circ Q$ and sum over the
patches. This maps each $2^{8}\times 2^{8}$ image to a $42$-dimensional
feature vector. As before, we apply PCA with $k=25$ and then train an LDA
classifier on the result. In the end, this max filtering approach
significantly outperforms the other classifiers we considered, as seen in
Table 1.

##  7 Discussion

Max filtering offers a rich source of invariants that can be used for a
variety of machine learning tasks. In this section, we discuss several
opportunities for follow-on work.

Separating. In terms of sample complexity, Corollary 13 gives that a generic
max filter bank of size $2d$ separates all $G$-orbits in $\mathbb{R}^{d}$
provided $G\leq\operatorname{O}(d)$ is finite. If $G\leq\operatorname{O}(d)$
is not topologically closed, then no continuous invariant (such as a max
filter bank) can separate all $G$-orbits. If $G\leq\operatorname{O}(d)$ is
topologically closed, then by Theorem 3.4.5 in [56], $G$ is algebraic, and so
Theorem 12 applies. We suspect that max filtering separates orbits in such
cases, but progress on this front will likely factor through Problem 15(a).
For computational complexity, how well does the max filtering separation
hierarchy (1) separate isomorphism classes of weighted graphs? Judging by [2],
we suspect that max filtering with a template of order $k$ and treewidth $t$
can be computed with runtime $e^{O(k)}n^{t+1}\log n$, which is polynomial in
$n$ when $k$ is logarithmic and $t$ is bounded. Which classes of graphs are
separated by such max filters? Also, can max filtering be used to solve real-
world problems involving graph data, or is the exponent too large to be
practical?

Bilipschitz. The proof of Lemma 20 contains our approach to finding Lipschitz
bounds on max filter banks. Our upper Lipschitz bound follows immediately from
the fact that each individual max filter is Lipschitz, and it does not require
the group to be finite. This bound is optimal for $G=\operatorname{O}(d)$, but
it is known to be loose for small groups like $G=\\{\pm\operatorname{id}\\}$.
We suspect that our lower Lipschitz bound leaves a lot more room for
improvement. In particular, we use the pigeonhole principle to pass from a sum
to a maximum so that we can leverage projective uniformity. A different
approach might lead to a tighter bound that does not require $G$ to be finite.
More abstractly, we suspect that separating implies bilipschitz, though the
bounds might be arbitrarily bad; see Problem 19.

Randomness. Our separating and bilipschitz results (Corollary 13 and Theorem
18) are not given in terms of explicit templates. Meanwhile, our Mallat-type
stability result (Theorem 22) requires a localized template, and we can
interpret the templates used in our weighted graph separation hierarchy (1) as
being localized, too. We expect that there are other types of structured
templates that would reduce the computational complexity of max filtering
(much like the structured measurements studied in compressed sensing [63, 59,
47] and phase retrieval [10, 38, 13, 40]). It would be interesting to prove
separating or bilipschitz results in such settings. More generally, can one
construct explicit templates for which max filtering is separating or
bilipschitz for a given group? Going the other direction, the reader may have
noticed that the plots in Figure 1 deviate from each other at the edges of the
interval. We expect that such deviations decay as the dimension grows, but
this requires further analysis. How can one analyze the behavior of random max
filters for other groups?

Max filtering networks. Example 32 opens the door to a variety of innovations
with max filtering. Here, instead of fixing a single group to mod out by, we
applied max filtering with a family of different groups. We selected
permutations over patches due to the computational simplicity of using Theorem
28, and we arranged the patches in a hierarchical structure so as to capture
behavior at different scales. Is there any theory to explain the performance
of this architecture? Are there other useful architectures in this vicinity,
perhaps by combining with neural networks?

## Acknowledgments

This work was initiated at the SOFT 2021: Summer of Frame Theory virtual
workshop. The authors thank Boris Alexeev for helpful discussions and Soledad
Villar for bringing the article [37] to their attention.

## References

  * [1] B. Alexeev, A. S. Bandeira, M. Fickus, D. G. Mixon, Phase retrieval with polarization, SIAM J. Imaging Sci. 7 (2014) 35–66.
  * [2] N. Alon, R. Yuster, U. Zwick, Color-coding, J. ACM 42 (1995) 844–856.
  * [3] V. Andrearczyk, P. F. Whelan, Using filter banks in convolutional neural networks for texture classification, Pattern Recognit. Lett. 84 (2016) 63–69.
  * [4] R. Balan, P. Casazza, D. Edidin, On signal reconstruction without phase, Appl. Comput. Harmon. Anal. 20 (2006) 345–356.
  * [5] R. Balan, C. B. Dock, Lipschitz Analysis of Generalized Phase Retrievable Matrix Frames, arXiv:2109.14522 (2021).
  * [6] R. Balan, N. Haghani, M. Singh, Permutation Invariant Representations with Applications to Graph Deep Learning, arXiv:2203.07546 (2022).
  * [7] R. Balan, Y. Wang, Invertibility and robustness of phaseless reconstruction, Appl. Comput. Harmon. Anal. 38 (2015) 469–488.
  * [8] A. S. Bandeira, B. Blum-Smith, J. Kileel, A. Perry, J. Weed, A. S. Wein, Estimation under group actions: Recovering orbits from invariants, arXiv:1712.10163 (2017).
  * [9] A. S. Bandeira, J. Cahill, D. G. Mixon, A. A. Nelson, Saving phase: Injectivity and stability for phase retrieval, Appl. Comput. Harmon. Anal. 37 (2014) 106–125.
  * [10] A. S. Bandeira, Y. Chen, D. G. Mixon, Phase retrieval from power spectra of masked signals, Inform. Inference 3 (2014) 83–102.
  * [11] T. Bendory, D. Edidin, W. Leeb, N. Sharon, Dihedral multi-reference alignment, IEEE Trans. Inform. Theory 68 (2022) 3489–3499.
  * [12] J. Bochnak, M. Coste, M.-F. Roy, Real algebraic geometry, Springer, 2013.
  * [13] B. G. Bodmann, N. Hammen, Stable phase retrieval with low-redundancy frames, Adv. Comput. Math. 41 (2015) 317–331.
  * [14] R. Bousseljot, D. Kreiseler, A. Schnabel, Nutzung der EKG-Signaldatenbank CARDIODAT der PTB über das Internet, Biomed. Tech. 40 (1995) 317–318, https://www.physionet.org/content/ptbdb/1.0.0/.
  * [15] J. Bruna, S. Mallat, Classification with scattering operators, CVPR 2011, 1561–1566.
  * [16] J. Bruna, S. Mallat, Invariant scattering convolution networks, IEEE Trans. Pattern Anal. Mach. Intell. 35 (2013) 1872–1886.
  * [17] J. Cahill, P. Casazza, I. Daubechies, Phase retrieval in infinite-dimensional Hilbert spaces, Trans. Amer. Math. Soc., Ser. B 3 (2016) 63–76.
  * [18] J. Cahill, A. Contreras, A. Contreras-Hip, Classifying Signals Under a Finite Abelian Group Action: The Finite Dimensional Setting, arXiv:1911.05862 (2019).
  * [19] J. Cahill, A. Contreras, A. Contreras-Hip, Complete set of translation invariant measurements with Lipschitz bounds, Appl. Comput. Harmon. Anal. 49 (2020) 521–539.
  * [20] E. J. Candès, Y. C. Eldar, T. Strohmer, V. Voroninski, Phase retrieval via matrix completion, SIAM Rev. 57 (2015) 225–251.
  * [21] E. J. Candès, X. Li, M. Soltanolkotabi, Phase retrieval via Wirtinger flow: Theory and algorithms, IEEE Trans. Inform. Theory 61 (2015) 1985–2007.
  * [22] E. J. Candès, T. Strohmer, V. Voroninski, Phaselift: Exact and stable signal recovery from magnitude measurements via convex programming, Comm. Pure Appl. Math. 66 (2013) 1241–1274.
  * [23] H. J. Caulfield, R. Haimes, Generalized matched filtering, Appl. Opt. 19 (1980) 181–183.
  * [24] H. J. Caulfield, W. T. Maloney, Improved discrimination in optical character recognition, Appl. Opt. 8 (1969) 2354–2356.
  * [25] A. Cayley, On the theory of linear transformations, Cambridge Math. J. 4 (1845) 193–209.
  * [26] S. Chen, E. Dobriban, J. Lee, A group-theoretic framework for data augmentation, NeurIPS 2020, 21321–21333.
  * [27] Y. Chen, E. J. Candès, Solving random quadratic systems of equations is nearly as easy as solving linear systems, Comm. Pure Appl. Math. 70 (2017) 822–883.
  * [28] D. C. Cireşan, U. Meier, L. M. Gambardella, J. Schmidhuber, Deep, big, simple neural nets for handwritten digit recognition, Neural Comput. 22 (2010) 3207–3220.
  * [29] A. Conca, D. Edidin, M. Hering, C. Vinzant, An algebraic characterization of injectivity in phase retrieval, Appl. Comput. Harmon. Anal. 38 (2015) 346–356.
  * [30] H. A. David, H. N. Nagaraja, Order statistics, John Wiley & Sons, 2004.
  * [31] L. Demanet, P. Hand, Stable optimizationless recovery from phaseless linear measurements, J. Fourier Anal. Appl. 20 (2014) 199–221.
  * [32] H. Derksen, G. Kemper, Computational invariant theory, Springer, 2015.
  * [33] M. Domokos, Degree bound for separating invariants of abelian groups, Proc. Amer. Math. Soc. 145 (2017) 3695–3708.
  * [34] M. Domokos, Separating monomials for diagonalizable actions, arXiv:2202.07002 (2022).
  * [35] D. L. Donoho, C. Grimes, Image manifolds which are isometric to Euclidean space, J. Math. Imaging Vis. 23 (2005) 5–24.
  * [36] E. S. Dufresne, Separating Invariants, Ph.D. thesis, 2008.
  * [37] N. Dym, S. J. Gortler, Low Dimensional Invariant Embeddings for Universal Geometric Learning, arXiv:2205.02956 (2022).
  * [38] Y. C. Eldar, P. Sidorenko, D. G. Mixon, S. Barel, O. Cohen, Sparse phase retrieval from short-time Fourier measurements, IEEE Signal Process. Lett. 22 (2014) 638–642.
  * [39] F. Gao, G. Wolf, M. Hirn, Geometric scattering for graph data analysis, ICML 2019, 2122–2131.
  * [40] D. Gross, F. Krahmer, R. Kueng, Improved recovery guarantees for phase retrieval from coded diffraction patterns, Appl. Comput. Harmon. Anal. 42 (2017) 37–64.
  * [41] J. Gu, et al., Recent advances in convolutional neural networks, Pattern Recognit. 77 (2018) 354–377.
  * [42] D. Hilbert, Über die Theorie der algebraischen Formen, Math. Ann. 36 (1890) 473–534.
  * [43] T. N. Huang, S. Villar, A Short Tutorial on The Weisfeiler-Lehman Test And Its Variants, ICASSP 2021, 8533–8537.
  * [44] C. Ingraham, America’s most gerrymandered congressional districts, The Washington Post, May 15, 2014, https://www.washingtonpost.com/news/wonk/wp/2014/05/15/americas-most-gerrymandered-congressional-districts/.
  * [45] M. A. Iwen, S. Merhi, M. Perlmutter, Lower Lipschitz bounds for phase retrieval from locally supported measurements, Appl. Comput. Harmon. Anal. 47 (2019) 526–538.
  * [46] W. B. Johnson, J. Lindenstrauss, Extensions of Lipschitz mappings into a Hilbert space, Contemp. Math. 26 (1984) 189–206.
  * [47] F. Krahmer, S. Mendelson, H. Rauhut, Suprema of chaos processes and the restricted isometry property, Comm. Pure Appl. Math. 67 (2014) 1877–1904.
  * [48] A. Krizhevsky, I. Sutskever, G. Hinton, ImageNet Classification with Deep Convolutional Neural Networks, NeurIPS 2012, 1097–1105.
  * [49] G. Kylberg, The Kylberg Texture Dataset v. 1.0, http://www.cb.uu.se/~gustaf/texture/.
  * [50] G. Kylberg, I.-M. Sintorn, Evaluation of noise robustness for local binary pattern descriptors in texture classification, EURASIP J. Image Video Process. 2013 (2013) 1–20.
  * [51] H. Makimoto, et al., Performance of a convolutional neural network derived from an ECG database in recognizing myocardial infarction, Sci. Rep. 10 (2020) 1–9.
  * [52] S. Mallat, Group invariant scattering, Comm. Pure Appl. Math. 65 (2012) 1331–1398.
  * [53] P. Mitteroecker, P. Gunz, Advances in geometric morphometrics, Evol. Biol. 36 (2009) 235–247.
  * [54] C. Morris, M. Ritzert, M. Fey, W. L. Hamilton, J. E. Lenssen, G. Rattan, M. Grohe, Weisfeiler and Leman go neural: Higher-order graph neural networks, AAAI 2019, 4602–4609.
  * [55] D. Mumford, J. Fogarty, F. Kirwan, Geometric invariant theory, Springer, 1994.
  * [56] A. L. Onishchik, E. B. Vinberg, Lie groups and algebraic groups, Translated by D. A. Leites, Springer, 1990.
  * [57] M. Perlmutter, F. Gao, G. Wolf, M. Hirn, Understanding graph neural networks with asymmetric geometric scattering transforms, arXiv:1911.06253 (2019).
  * [58] A. Perry, J. Weed, A. S. Bandeira, P. Rigollet, A. Singer, The sample complexity of multireference alignment, SIAM J. Math. Data Science 1 (2019) 497–517.
  * [59] G. E. Pfander, H. Rauhut, J. A. Tropp, The restricted isometry property for time–frequency structured random matrices, Probab. Theory Related Fields 156 (2013) 707–737.
  * [60] A. Rahimi, B. Recht, Weighted sums of random kitchen sinks: Replacing minimization with randomization in learning, NeurIPS 2008, 858.
  * [61] R. T. Rockafellar, Convex Analysis, Princeton U. Press, 1970.
  * [62] Y. Rong, Y. Wang, Z. Xu, Almost everywhere injectivity conditions for the matrix recovery problem, Appl. Comput. Harmon. Anal. 50 (2021) 386–400.
  * [63] M. Rudelson, R. Vershynin, On sparse reconstruction from Fourier and Gaussian measurements, Comm. Pure Appl. Mathematics 61 (2008) 1025–1045.
  * [64] R. Sato, A survey on the expressive power of graph neural networks, arXiv:2003.04078 (2020).
  * [65] D. Silver, et al., Mastering the game of Go with deep neural networks and tree search, Nature 529 (2016) 484.
  * [66] P. Y. Simard, D. Steinkraus, J. C. Platt, Best practices for convolutional neural networks applied to visual document analysis, ICDAR 2003, 1–6.
  * [67] C. D. Sogge, Fourier integrals in classical analysis, Cambridge U. Press, 2017.
  * [68] N. Switzer, USA 2019 Congressional District Shape Files, https://www.kaggle.com/datasets/nswitzer/usa-2019-congressional-district-shape-files.
  * [69] S. Villar, D. W. Hogg, K. Storey-Fisher, W. Yao, B. Blum-Smith, Scalars are universal: Equivariant machine learning, structured like classical physics, NeurIPS 2021, 28848–28863.
  * [70] C. Vinzant, A small frame and a certificate of its injectivity, SampTA 2015, 197–200.
  * [71] I. Waldspurger, Exponential decay of scattering coefficients, SampTA 2017, 143–146.
  * [72] I. Waldspurger, A. d’Aspremont, S. Mallat, Phase recovery, maxcut and complex semidefinite programming, Math. Program. 149 (2015) 47–81.
  * [73] Y. Wang, Z. Xu, Generalized phase retrieval: Measurement number, matrix recovery and beyond, Appl. Comput. Harmon. Anal. 47 (2019) 423–446.
  * [74] G. A. Watson, Characterization of the subdifferential of some matrix norms, Linear Algebra Appl. 170 (1992) 33–45.
  * [75] B. Weisfeiler, A. Leman, The reduction of a graph to canonical form and the algebra which appears therein, Nauchno-Tekhn. Inform. 2 (1968) 12–16.
  * [76] J. Wood, Invariant pattern recognition: A review, Pattern Recognit. 29 (1996) 1–17.
  * [77] K. Xu, W. Hu, J. Leskovec, S. Jegelka, How Powerful are Graph Neural Networks?, ICLR 2019.
  * [78] D. Zhang, G. Lu, Review of shape representation and description techniques, Pattern Recognit. 37 (2004) 1–19.
  * [79] D. Zou, R. Balan, M. Singh, On Lipschitz bounds of general convolutional neural networks, IEEE Trans. Inform. Theory 66 (2019) 1738–1759.