# On the hardness of learning under symmetries

Bobak T. Kiani⋄ John A. Paulson School of Engineering and Applied Sciences,
Harvard; Department of Electrical Engineering and Computer Science, MIT;
e-mail<EMAIL_ADDRESS>Thien Le⋄ Department of Electrical Engineering
and Computer Science, MIT; e-mail<EMAIL_ADDRESS>Hannah Lawrence⋄ Department
of Electrical Engineering and Computer Science, MIT; e-mail<EMAIL_ADDRESS>Stefanie Jegelka Department of Electrical Engineering and Computer Science,
MIT; School of Computation, Information and Technology, TU Munich; e-mail:
<EMAIL_ADDRESS>Melanie Weber John A. Paulson School of Engineering and Applied
Sciences, Harvard; e-mail<EMAIL_ADDRESS>
###### Abstract

We study the problem of learning equivariant neural networks via gradient
descent. The incorporation of known symmetries (“equivariance”) into neural
nets has empirically improved the performance of learning pipelines, in
domains ranging from biology to computer vision. However, a rich yet separate
line of learning theoretic research has demonstrated that actually learning
shallow, fully-connected (i.e. non-symmetric) networks has exponential
complexity in the correlational statistical query (CSQ) model, a framework
encompassing gradient descent. In this work, we ask: are known problem
symmetries sufficient to alleviate the fundamental hardness of learning neural
nets with gradient descent? We answer this question in the negative. In
particular, we give lower bounds for shallow graph neural networks,
convolutional networks, invariant polynomials, and frame-averaged networks for
permutation subgroups, which all scale either superpolynomially or
exponentially in the relevant input dimension. Therefore, in spite of the
significant inductive bias imparted via symmetry, actually learning the
complete classes of functions represented by equivariant neural networks via
gradient descent remains hard.

$\diamond$$\diamond$footnotetext: co-first authors

##  1 Introduction

In recent years, the purview of machine learning has expanded to non-
traditional domains with geometric input types, from graphs to sets to point
clouds. Correspondingly, it is now common practice to tailor neural
architectures to the particular symmetries of the input – graph neural
networks are invariant to permutations of the input nodes, for example, while
networks operating on molecules and point clouds are invariant to permutation,
translation, and rotation. Empirically, encoding such structure has led to
computational benefits in applications (Wang et al., 2021; Batzner et al.,
2022; Bronstein et al., 2021). From the perspective of generalization,
previous research has quantified the benefits of imposing symmetry during
learning (Long and Sedghi, 2019; Bietti et al., 2021; Sannai et al., 2021; Mei
et al., 2021), often achieving rather tight bounds for simple models (Elesedy,
2021; Tahmasebi and Jegelka, 2023). At their core, these formal statements are
bounds on sample complexity, i.e. how much data is needed to learn a given
task. In contrast, the effect of symmetries on the computational complexity of
learning algorithms has not been previously studied. Broadly speaking,
generalization bounds are necessary but not sufficient to show efficiently
learnability, as there can be exponentially large gaps between sample
complexity and runtime lower bounds.

Indeed, an active line of research in learning theory studies the hardness of
learning fully-connected neural networks via “correlational statistical query”
(CSQ) algorithms (defined in Sec. 2.1), which notably encompass gradient
descent. In particular, for a Gaussian data distribution, Diakonikolas et al.
(2020) proved exponential lower bounds for learning shallow neural nets. Such
works provide valuable impossibility results, demonstrating that one cannot
hope to efficiently learn any function represented by a small neural net under
simple data distributions.

In this work, we ask: does invariance provide a restrictive enough subclass of
neural network functions to circumvent these impossibility results? Our main
finding, shown in various settings, is that an inductive bias towards symmetry
is not strong enough to alleviate the fundamental hardness of learning neural
nets via gradient descent. We view this result as a guide for future
theoretical work: additional function structure is generally needed, even in
the invariant setting, to achieve guarantees of efficient learnability. Our
proof techniques largely rely on extending techniques from the general setting
to symmetric function classes. In many of our lower bounds, imposition of
symmetry on a function class reduces the computational complexity of learning
by a factor proportional to the size of the group. Some of our lower bounds
use simple symmetrizations of existing hard families, while others require
bona fide new classes of hard invariant networks. For example, intuition from
message passing GNN architectures seems to indicate that complexity would grow
exponentially with the number of such message passing steps; however, our
results show that even a single message passing layer gathers sufficient
features to form exponentially hard functions (Liao et al., 2020).

#### Our Contributions

We consider various interrelated questions on the computational hardness of
learning, focusing primarily on how hard it is to learn data generated by
invariant architectures.

Question (primary). Given Gaussian i.i.d. inputs in $n$ variables, how hard is
it to learn the class of shallow single hidden layer invariant architectures
in the correlational statistical query model?

We answer this question in various settings. As a warm-up, we give general
(SQ) lower bounds for learning invariant Boolean functions in Sec. 3. Next, we
show that graph neural networks with a single message passing layer are
exponentially hard to learn, either in the number of nodes or feature
dimension (Theorem 3 and Theorem 4, respectively). Hardness in the feature
dimension results from careful symmetrization of the hard functions in
Diakonikolas et al. (2020). This technique did not extend to the node
dimension hardness result in Theorem 3, for which we instead constructed a
customized set of hard functions whose outputs are parities in the degree
counts of the graph. Since adjacency matrices of graphs are often Boolean
valued, this setting also offers general SQ lower bounds scaling exponentially
with the number of nodes $n$. For subgroups of permutations, we prove CSQ
hardness results for “symmetrized” (frame-averaged) neural network classes
encapsulating CNNs (Theorem 7 of Sec. 5). The proofs here are based on
constructions in Goel et al. (2020) and Diakonikolas et al. (2020), ensuring
that properly chosen frames retain the underlying structure necessary (e.g.,
sign invariance) to guarantee orthogonality of hard functions.

These lower bounds constitute our main results, but we provide two
complementary results to provide an even more complete picture. First, we show
that classes of sparse invariant polynomials can be learned efficiently, but
not by gradient descent, via a straightforward extension of the algorithm in
Andoni et al. (2014) for arbitrary polynomials (Sec. 6). This learning
algorithm is SQ, and we give various $n^{\Omega(\log n)}$ CSQ lower bounds for
learning over $n^{\Omega(\log n)}$ many orthogonal invariant polynomials of
degree at most $O(\log(n))$. Second, we take a classical computational
complexity perspective to establish the NP hardness of learning GNNs. The work
of Blum and Rivest (1988) established that even 3-node feedforward networks
are $\mathsf{NP}$ hard to train, and later expanded to average case complexity
for improper learning by various works (Daniely et al., 2014; Daniely, 2016;
Klivans and Kothari, 2014). Although it is perhaps not surprising that these
hardness results can be extended to invariant architectures, for the sake of
completeness (and since we were unaware of an extension in the literature), in
Sec. 4.3 we give an $\mathsf{NP}$ hardness result for proper learning of the
weights of a GNN architecture. Finally, Sec. 7 provides a few experiments
verifying that the hard classes of invariant functions we propose are indeed
difficult to learn.

###  1.1 Related work

A long line of research studies the hardness of learning neural networks. We
briefly review the most relevant works here, and refer the reader to App. A
for a more holistic review.

Our work is largely motivated to extend hardness results for feedforward
networks to equivariant and symmetric neural networks. Early results showed
that there exist rather artificial distributions of data for which feedforward
networks are hard to learn in worst-case settings (Judd, 1987; Blum and
Rivest, 1988; Livni et al., 2014). Interest later grew into studying hardness
for more natural settings via the statistical query framework (Kearns, 1998;
Shamir, 2018). Goel et al. (2020) showed superpolynomial CSQ lower bounds for
learning single hidden layer ReLU networks, subsequently strengthened to
exponential lower bounds in Diakonikolas et al. (2020). We use the hard class
of networks in Diakonikolas et al. (2020) as a basis for the frame averaged
networks in our work. Hardness results for learning feedforward neural
networks have also been shown by proving reductions to average case or
cryptographically hard problems (Chen et al., 2022a; Daniely and Vardi, 2020).

Though our work is focused on runtime lower bounds, some algorithms exist for
learning restricted classes of geometric networks efficiently such as CNNs
(Brutzkus and Globerson, 2017; Du et al., 2017, 2018a). Zhang et al. (2020);
Li et al. (2021) give algorithms learning one hidden layer GNNs via gradient
based algorithms, which are efficient when certain factors such as the
condition number of the weights are bounded. We should also note that for
practical separations between SQ and CSQ algorithms, Andoni et al. (2014) show
that learning sparse polynomials of degree $d$ in $n$ variables requires
$\Omega(n^{d})$ CSQ queries, but often only $\tilde{O}(nd)$ SQ queries. We
extend this to polynomials in the invariant polynomial ring in Sec. 6.

##  2 Background and Notation

We use $a$, ${\bm{a}}$, and ${\bm{A}}$ to denote scalars, vectors, and
matrices. ${\bm{I}}_{n}$ denotes the identity matrix of dimension $n$. We
denote groups by $G$ and graphs by ${\mathcal{G}}$. We denote the normal
distribution supported over $\mathbb{R}^{n}$ with mean
${\bm{v}}\in\mathbb{R}^{n}$ and covariance matrix
${\bm{M}}\in\mathbb{R}^{n\times n}$ as $\mathcal{N}({\bm{v}},{\bm{M}})$, often
abbreviated as $\mathcal{N}$ for $\mathcal{N}(\bm{0},{\bm{I}}_{n})$. We also
denote by $\mathbb{G}_{n}$ the set of graphs on $n$ vertices and, when clear
from context, $\mathcal{E}$ an arbitrary distribution over $\mathbb{G}_{n}$.
To ease notation, we often use $[n]$ to denote the set $\\{1,\dots,n\\}$.

Let $\mathcal{X}$ be an input space and $\mathcal{D}$ a distribution on
$\mathcal{X}$. Given functions $f,g:\mathcal{X}\to\mathbb{R}$, their inner
product is $\langle f,g\rangle_{\mathcal{D}}=\mathbb{E}_{\mathcal{D}}[fg]$
with corresponding norm $\|f\|_{\mathcal{D}}=\sqrt{\langle
f,f\rangle_{\mathcal{D}}}$. For functions $f,g:\mathcal{X}\to\\{-1,+1\\}$
whose outputs are Boolean, the classification error is
$\mathbb{P}_{{\bm{x}}\sim\mathcal{D}}[f({\bm{x}})\neq g({\bm{x}})]$.

###  2.1 SQ Learning framework

The well-studied statistical query (SQ) model offers a restricted query
complexity based model for proving hardness, which encapsulates most
algorithms in practice (Kearns, 1998; Reyzin, 2020). Given a joint
distribution $\mathcal{D}$ on input/output space
$\mathcal{X}\times\mathcal{Y}$, any SQ algorithm is composed of a set of
queries. Each query takes as input a function
$g:\mathcal{X}\times\mathcal{Y}\to[-1,1]$ and tolerance parameter $\tau>0$,
and returns a value $\operatorname{SQ}(g,\tau)$ in the range:

$\mathbb{E}_{({\bm{x}},y)\sim\mathcal{D}}\left[g({\bm{x}},y)\right]-\tau\leq\operatorname{SQ}(g,\tau)\leq\mathbb{E}_{({\bm{x}},y)\sim\mathcal{D}}\left[g({\bm{x}},y)\right]+\tau.$
(1)

A special class of queries are correlational statistical queries (CSQ), where
the query function $g:\mathcal{X}\to[-1,1]$ is only a function of ${\bm{x}}$,
and the oracle returns
$\mathbb{E}_{({\bm{x}},y)\sim\mathcal{D}}\left[g({\bm{x}})\cdot y\right]$ the
correlation of $g$ with $y$ up to error $\tau$.

Hardness is quantified as the number of queries required to learn a function
$y=c^{*}(x)$, drawn from a function class $\mathcal{C}$ up to a desired
maximal error. Notably, SQ hardness results imply hardness for any gradient
based algorithm, as gradients with respect to a loss can be captured by
statistical queries. Furthermore, only CSQ oracles are required for the mean-
squared error (MSE) loss.

######  Example 1 (GD from $\operatorname{CSQ}$).

Given a function $N_{\theta}:\mathcal{X}\to\mathbb{R}$ with parameters
$\theta$, gradients of the MSE loss with respect to $\theta$ can be estimated
with a correlational statistical query for each parameter:

$\mathbb{E}_{({\bm{x}},y)\sim\mathcal{D}}\left[\nabla_{\theta}\tfrac{1}{2}(y-N_{\theta}({\bm{x}}))^{2}\right]=\underbrace{\mathbb{E}\left[N_{\theta}({\bm{x}})\nabla_{\theta}N_{\theta}({\bm{x}})\right]}_{y-\text{
independent}}-\underbrace{\mathbb{E}\left[y\nabla_{\theta}N_{\theta}({\bm{x}})\right]}_{\text{CSQ
with}\nabla_{\theta}N_{\theta}}.$ (2)

Access to the distribution $\mathcal{D}$ is typically provided through
samples, and the tolerance $\tau$ in part captures the error in statistical
sampling of quantities such as the gradient. In fact, by drawing
$O(\log(1/\delta)/\tau^{2})$ samples, standard Hoeffding bounds guarantee
estimates of a given query are within tolerance $\tau$ with probability
$1-\delta$. In the SQ (CSQ) framework, the complexity of learning is
determined by the tolerance bound and number of queries (loosely corresponding
to “steps” in an algorithm) needed to learn a function class.

######  Definition 1 (SQ (CSQ) Learning).

Given a function class $\mathcal{C}$ and distribution $\mathcal{D}$ over input
space $\mathcal{X}$, an algorithm SQ (CSQ) learns $\mathcal{C}$ up to
classification error ($\ell_{2}$ squared loss) $\epsilon$ if, given only SQ
(CSQ) oracle access to $({\bm{x}},c^{*}({\bm{x}})),{\bm{x}}\sim\mathcal{D}$
for some unknown $c^{*}\in\mathcal{C}$, the algorithm outputs a function $f$
such that $\mathbb{P}_{{\bm{x}}\sim\mathcal{D}}\left[f({\bm{x}})\neq
c^{*}({\bm{x}})\right]\leq\epsilon$ (i.e.
$\|f-c^{*}\|_{\mathcal{D}}\leq\epsilon$).

SQ lower bounds for Boolean functions are typically shown by finding a set of
roughly orthogonal functions that lower bound the so-called statistical
dimension. For real-valued function classes, lower bounds on the statistical
dimension imply CSQ lower bounds. Applying more general SQ lower bounds is
challenging, unless one reduces the problem of learning a real-valued function
to that of learning a subset therein of (say) Boolean functions within the
class (Chen et al., 2022a). We cover this in more detail and also discuss
other learning theoretic frameworks in App. B.

##  3 Warm up: invariant Boolean functions

Before detailing our main results concerning real-valued functions, we explore
the foundational context of Boolean functions where the SQ formalism
originated (Kearns, 1998; Blum et al., 1994). This illustrative analysis will
serve as a warm-up for the more complex scenarios later. It is known that many
classes of Boolean functions have exponential query complexity (e.g., the
class of parity functions) arising from the orthogonality properties of
Boolean functions. Enforcing invariance under a given group $G$, the number of
orbits of the $2^{n}$ bitstrings gives an analogous hardness metric, though
care must be taken in dealing with the distribution of these orbits.

More formally, in the Boolean setting, let input distribution $\mathcal{D}$ be
uniform over $\mathcal{X}=\\{-1,+1\\}^{n}$. For a group $G$ with
representation $\rho$, let $\mathcal{O}_{\rho}=\\{\\{\rho(g)\cdot{\bm{x}}:g\in
G\\}:{\bm{x}}\in\\{-1,+1\\}^{n}\\}$ denote the orbits of the inputs under the
representation. For any given orbit $O\in\mathcal{O}_{\rho}$, let
${\bm{x}}_{O}$ denote an orbit representative (i.e., some fixed element in the
set $O$). We can represent any symmetric function
$f:\mathcal{O}_{n}\to\\{-1,+1\\}$ as a function of these orbits, where the
correlation of two functions $f,g$ is

$\langle
f,g\rangle_{\mathcal{D}}=\mathbb{E}\left[fg\right]=2^{-n}\sum\nolimits_{O_{k}\in\mathcal{O}_{\rho}}|O_{k}|f({\bm{x}}_{O_{k}})g({\bm{x}}_{O_{k}}).$
(3)

The probability $p_{\mathcal{O}_{\rho}}:\mathcal{O}_{\rho}\to[0,1]$ of a
random bitstring falling in orbit $O_{K}\in\mathcal{O}_{\rho}$ is equal to
$p_{\mathcal{O}_{\rho}}(O_{k})=|O_{k}|/2^{n}$. With this notation, we have a
general lower bound on the SQ hardness of learning invariant Boolean
functions.

######  Theorem 2 (Boolean SQ hardness).

For a given symmetry group $G$ with representation $\rho:G\to
GL(\\{-1,+1\\}^{n})$, let
$\|p_{\mathcal{O}_{\rho}}\|\coloneqq\big{(}{\sum_{O_{k}\in\mathcal{O}_{\rho}}\left(\frac{|O_{k}|}{2^{n}}\right)^{2}}\big{)}^{1/2}$
and let $\mathcal{H}_{\rho}$ be the class of symmetric Boolean functions,
defined as

$\mathcal{H}_{\rho}=\left\\{f:\\{-1,+1\\}^{n}\to\\{-1,+1\\}:\;\forall g\in
G,\forall{\bm{x}}\in\\{-1,+1\\}^{n}:f(\rho(g)\cdot{\bm{x}})=f({\bm{x}})\right\\}.$
(4)

Any SQ learner capable of learning $\mathcal{H}_{\rho}$ up to sufficiently
small classification error probability $\epsilon$ ($\epsilon<1/4$ suffices)
with queries of tolerance $\tau$ requires at least
$\tau^{2}\|p_{\mathcal{O}_{\rho}}\|^{-2}/2$ queries.

###### Proof sketch.

We show that $1-\|p_{\mathcal{O}_{\rho}}\|^{2}$ is the probability that, over
independently drawn inputs ${\bm{x}},{\bm{x}}^{\prime}$, the distribution
$(f({\bm{x}}),f({\bm{x}}^{\prime}))$ is uniformly random over $f$ drawn
uniformly from $\mathcal{H}_{\rho}$. From a proof technique in Chen et al.
(2022a), this connects the task to one of distinguishing distributions, for
which the SQ complexity is at least the stated amount. See App. C for complete
proof. ∎

For a more direct lower bound, we can use Hölder’s inequality such that
$\|p_{\mathcal{O}_{\rho}}\|^{2}\leq\frac{1}{2^{n}}\max_{O_{k}\in\mathcal{O}_{\rho}}|O_{k}|\leq|G|2^{-n}$.
This gives query complexity of at least $2^{n-1}|G|^{-1}\tau^{2}$. For
commonly studied groups, assuming $\tau=O(1)$, the SQ complexity generally
aligns with the number of orbits. Table 1 summarizes these hardness lower
bounds for commonly studied groups111Certain sub-classes of
$\mathcal{H}_{\rho}$ may suffice to attain hardness, much like how Parity
functions suffice to prove exponential lower bounds in the traditional Boolean
setting. We do not consider this strengthening here..

Group |  | $\max_{O_{k}\in\mathcal{O}_{\rho}}|O_{k}|$ | Query Complexity
---|---|---|---
Symmetric group on $n$ bits |  | $\frac{{\binom{n}{n/2}}}{2^{n}}=\frac{1}{\Theta(\sqrt{n})}$ | $\Theta(\sqrt{n})$
Symmetric group on $n\times n$ graphs |  | $\frac{n!}{2^{n^{2}}}=2^{-n^{2}+n\log n+O(n)}$ | $\Omega(2^{O(n^{2})})$
Cyclic group on $n$ bits |  | $n$ | $\Omega(2^{n}/n)$
Table 1: Query complexity of learning common invariant Boolean function
classes.

##  4 Lower bounds for GNNs

We now show statistical query lower bounds for graph neural networks, which
are invariant to node permutations, in three settings: (1) SQ lower bounds
scaling with the number of nodes $n$ for Erdős–Rényi distributed graphs with
trivial node features, (2) CSQ lower bounds scaling in the node feature
dimension for a fixed graph, and (3) a simple extension of $\mathsf{NP}$
hardness in a proper learning task (deciding if fitting a training set with a
given GNN is possible).

###  4.1 Hardness in number of nodes

Here, we consider a class of two hidden layer GNNs that are only a function of
the adjacency matrix, and show hardness in the number of nodes $n$. This
covers a commonly used procedure where one learns node-level equivariant
features in the first layer, aggregates these features, and passes them
through an MLP (Dwivedi et al., 2020). We consider a Boolean function class of
GNNs over the uniform distribution $\operatorname{Unif}(\\{0,1\\}^{n\times
n})$, which can be viewed as the Erdős–Rényi random graph model with
$p=0.5$.222Our exponential lower bounds arise from the wide range of possible
degrees of nodes in the graph. Though we do not generalize this result, this
fact means that it is likely extendable to restricted sparse Erdős–Rényi
models $G(n,p_{n})$ where $p_{n}=\omega(1/n)$ (i.e. average degree grows
arbitrarily with $n$). In this construction, node features
${\bm{x}}\in\mathbb{R}^{n}$ are trivial and always set to ${\bm{x}}=\bm{1}$.

####  $2$ hidden layer GNN family

We consider two layer GNNs $f=f^{(2)}\circ f^{(1)}$. These take the commonly
used form of message passing $f^{(1)}:\\{0,1\\}^{n\times
n}\to\mathbb{R}^{k_{1}}$ which aggregates $k_{1}$ permutation invariant
features for the graph, followed by a single hidden layer ReLU MLP
$f^{(2)}:\mathbb{R}^{k_{1}}\to\\{0,1\\}$.

$\begin{split}[f_{{\bm{a}},{\bm{b}}}^{(1)}({\bm{A}})]_{i}&=\bm{1}_{n}^{\top}\sigma\left(a_{i}+b_{i}{\bm{A}}{\bm{x}}\right)\quad\text{(output
of channel $i\in[k_{1}]$)}\\\
f_{{\bm{u}},{\bm{v}},{\bm{W}}}^{(2)}({\bm{h}})&=\sum\nolimits_{i=1}^{k_{2}}u_{i}\sigma(\langle{\bm{W}}_{:,i},{\bm{h}}\rangle+v_{i}),\end{split}$
(5)

with subscripts for trainable weights. For our hard class of functions,
$k_{1}=O(n)$ and $k_{2}=O(n)$, so there are at most $O(n^{2})$ parameters
(proportional to the number of edges).

#### Family of hard functions

We define the hard functions in terms of degree counts
${\bm{c}}_{{\bm{A}}}\in[n]^{n+1}$, where entry $[{\bm{c}}_{{\bm{A}}}]_{i}$
counts the number of nodes that have $i-1$ outgoing edges in ${\bm{A}}$:

$[{\bm{c}}_{{\bm{A}}}]_{i}=\sum\nolimits_{k=1}^{n}\mathbb{1}\left[[{\bm{A}}\bm{1}]_{k}=i-1\right].$
(6)

The hard functions in $\mathcal{H}_{ER,n}$ are enumerated over subsets
$S\subseteq[n+1]$ and a bit $b\in\\{0,1\\}$:

$\mathcal{H}_{ER,n}=\\{g_{S,b}:S\subseteq[n+1],b\in\\{0,1\\}\\},\quad
g_{S,b}({\bm{A}})=b+\sum\nolimits_{i\in S}[{\bm{c}}_{{\bm{A}}}]_{i}\mod 2.$
(7)

$g_{S,b}({\bm{A}})$ is a sort of parity function supported over integers in
${\bm{c}}_{{\bm{A}}}$ in $\mathcal{S}$. We show in Lemma 15 that the GNNs in
Eq. 5 can construct these functions. Additionally, the class
$\mathcal{H}_{ER,n}$ is Boolean, so hardness lower bounds apply in the general
SQ model.

######  Theorem 3 (SQ hardness of $\mathcal{H}_{ER,n}$).

Any SQ learner capable of learning $\mathcal{H}_{ER,n}$ up to classification
error probability $\epsilon$ sufficiently small ($\epsilon<1/4$ suffices) with
queries of tolerance $\tau$ requires at least
$\Omega\left(\tau^{2}\exp(n^{\Omega(1)})\right)$ queries.

###### Proof sketch.

We form networks that in the message passing layer calculate
${\bm{c}}_{{\bm{A}}}$, which is passed to the MLP calculating Eq. 7 as a sum
of ReLUs via hat-like functions that mimic parities. Based on concentration
properties of the Erdős–Rényi model, we show that for some $p<0.5$, at least
$\Theta(n^{p})$ entries of ${\bm{c}}_{\bm{A}}$ will have values scaling
roughly as $O(\sqrt{n})$ and the probability that any such entry is odd
converges to $1/2$. Also using concentration properties, we condition on a
high probability region (where probability of these entries being odd are
roughly independent from each other), and use hardness of learning parity
functions to get SQ lower bounds. ∎

###  4.2 Hardness in feature dimension

In this section, we derive correlational statistical query lower bounds (CSQ)
for commonly used graph neural network architectures (GNNs), and state the
corresponding learning-theoretic hardness result. The techniques in this
section extend those of Diakonikolas et al. (2020), resulting in an
exponential (in feature dimension) lower bound for learning one-hidden layer
GNNs.

####  $1$-hidden layer GNN

For a fixed directed, unweighted graph ${\mathcal{G}}$ with $n\in\mathbb{N}$
vertices, let ${\bm{A}}({\mathcal{G}})\in\mathbb{R}^{n\times n}$ be the
adjacency matrix or the Laplacian of ${\mathcal{G}}$. For input feature
dimension $d\in\mathbb{N}$ and width parameter $k\in\mathbb{N}$, we consider
the following set of functions:

$F_{{\mathcal{G}}}^{d,k}:=\left\\{f:\mathbb{R}^{n\times
d}\to\mathbb{R},\>f({\bm{X}})={\bm{1}}_{n}^{\top}\sigma({\bm{A}}({\mathcal{G}}){\bm{X}}{\bm{W}}){\bm{a}}\mid{\bm{W}}\in\mathbb{R}^{d\times
2k},{\bm{a}}\in\mathbb{R}^{2k}\right\\},$ (8)

for some nonlinearity $\sigma$. When the graph itself is part of the input, we
have the function class:

$F_{n}^{d,k}:=\left\\{f:\mathbb{R}^{n\times
d}\times\mathbb{G}_{n}\to\mathbb{R},\>f({\bm{X}},{\mathcal{G}})={\bm{1}}_{n}^{\top}\sigma({\bm{A}}({\mathcal{G}}){\bm{X}}{\bm{W}}){\bm{a}}\mid{\bm{W}}\in\mathbb{R}^{d\times
2k},{\bm{a}}\in\mathbb{R}^{2k}\right\\}.$ (9)

#### Family of hard functions

First, we define functions $f_{{\mathcal{G}}}:\mathbb{R}^{n\times
2}\to\mathbb{R}$ and $f_{n}:\mathbb{R}^{n\times
2}\times\bar{{\mathcal{G}}}_{n}\to\mathbb{R}$ as

$\displaystyle
f_{{\mathcal{G}}}({\bm{X}})=\bm{1}_{n}^{\top}\sigma\left({\bm{A}}({\mathcal{G}}){\bm{X}}{\bm{W}}^{*}\right){\bm{a}}^{*},\quad
f_{n}({\bm{X}},{\mathcal{G}})=\bm{1}_{n}^{\top}\sigma\left({\bm{A}}({\mathcal{G}}){\bm{X}}{\bm{W}}^{*}\right){\bm{a}}^{*},$
(10)

$\displaystyle\text{ where }\quad{\bm{W}}^{*}$
$\displaystyle:=\begin{bmatrix}\left(\cos(\pi j/k)\right)^{\top}_{j\in[2k]}\\\
\left(\sin(\pi
j/k)\right)^{\top}_{j\in[2k]}\end{bmatrix}\in\mathbb{R}^{2\times
2k},\qquad{\bm{a}}^{*}=\left((-1)^{j}\right)_{j\in[2k]}\in\mathbb{R}^{2k}.$
(11)

Given a set of matrices $\mathcal{B}\subset\mathbb{R}^{d\times 2}$, our family
of hard functions can now be written as:

$\displaystyle
C^{\mathcal{B}}_{{\mathcal{G}}}:=\left\\{g^{\bm{B}}_{{\mathcal{G}}}:{\bm{X}}\mapsto\frac{f_{{\mathcal{G}}}({\bm{X}}B)}{\|f_{{\mathcal{G}}}\|_{\mathcal{N}}}\mid
B\in\mathcal{B}\right\\},\quad
C^{\mathcal{B}}_{n}:=\left\\{g^{\bm{B}}_{n}:({\bm{X}},{\mathcal{G}})\mapsto\frac{f_{n}({\bm{X}}B,{\mathcal{G}})}{\|f_{n}\|_{\mathcal{N}\times\mathcal{E}}}\mid
B\in\mathcal{B}\right\\}.$ (12)

Below, we show exponential lower bounds 333To have a meaningful set of
functions, $g^{\bm{B}}_{\mathcal{G}}$ and $g^{\bm{B}}_{n}$ must not vanish,
which holds when $\sigma$ is not a low-degree polynomial ($\rm{ReLU}$
suffices, see Remark 12 of Diakonikolas et al. (2020)) and
${\mathcal{G}}\not\equiv\emptyset$ (graph with no edges). for learning
$F_{{\mathcal{G}}}^{d,k}$ and $F_{n}^{d,k}$:

######  Theorem 4 (Exponential CSQ lower bound for GNNs).

For any number of vertices $n$ independent of input dimension $d$ and width
parameter $k$, let $\epsilon>0$ be a sufficiently small error constant and
$\mathcal{E}$ a distribution over $\mathbb{G}_{n}$, independent of
$\mathcal{N}$, with $\Pr_{\mathcal{E}}(\\{\emptyset\\})<\Omega(1)<1$. Then
there exists a set $\mathcal{B}$ of size at least $2^{\Omega(d^{\Omega(1)})}$
such that any CSQ algorithm that queries from oracles of concept $f\in
C^{\mathcal{B}}_{{\mathcal{G}}}$ (resp. $C^{\mathcal{B}}_{n}$) and outputs a
hypothesis $h$ with $\|f-h\|_{\mathcal{N}}\leq\epsilon$ (resp.
$\|f-h\|_{\mathcal{N}\times\mathcal{E}}\leq\epsilon$) requires either
$2^{d^{\Omega(1)}}$ queries or at least one query with precision
$d^{-\Omega(k)}+2^{-d^{\Omega(1)}}$.

###### Proof sketch.

Full details of the proof can be found in App. E. From Lemma 17 of
Diakonikolas et al. (2020), there exists a set
$\mathcal{B}\subset\mathbb{R}^{d\times 2}$ of size at least
$2^{\Omega(d^{\Omega(1)})}$ such that

$\|{\bm{B}}^{\top}{\bm{B}}^{\prime}\|_{2}=O(d^{-\Omega(1)})<1\text{ if
}{\bm{B}}\neq{\bm{B}}^{\prime}$ (13)

and ${\bm{B}}^{\top}{\bm{B}}=\mathbf{1}_{2}$ which we use to index
$C_{{\mathcal{G}}}^{\mathcal{B}}$. We design a specialized class of invariant
orthogonal Hermite polynomials, such that the correlation in low degree
moments vanishes and is bounded as:

$\left|\langle
f_{{\mathcal{G}}}(\cdot{\bm{B}}),f_{{\mathcal{G}}}(\cdot{\bm{B}}^{\prime})\rangle_{\mathcal{N}}\right|\leq\|{\bm{B}}^{\top}{\bm{B}}^{\prime}\|^{k}_{2}\sum\nolimits_{m>k}\|f^{[m]}_{{\mathcal{G}}}\|^{2}_{\mathcal{N}}=\|{\bm{B}}^{\top}{\bm{B}}^{\prime}\|^{k}_{2}\|f_{{\mathcal{G}}}\|^{2}_{\mathcal{N}}.$
(14)

Finally, we use Eq. 13 to get the desired almost uncorrelatedness of elements
in $C_{{\mathcal{G}}}^{\mathcal{B}}$ and use a simple total probability
argument to extend this to $C_{n}^{\mathcal{B}}$. ∎

###  4.3 $\mathsf{NP}$ hardness of proper learning of GNNs

The classic work of Blum and Rivest (1988) showed that there exist datasets
for which determining whether or not a 3-node neural network can fit that
dataset is $\mathsf{NP}$ hard. It is straightforward to prove a similar result
for invariant networks. We map the task of training a single hidden layer GNN
to the $\mathsf{NP}$ hard problem of learning halfspaces with noise (Guruswami
and Raghavendra, 2009; Feldman et al., 2012). Here, one is given a set of $N$
labeled examples $\\{{\bm{x}}_{i},y_{i}\\}_{i=1}^{N}$ and must determine
whether there exists a halfspace parameterized by vector ${\bm{v}}$ and
constant $\theta$ that correctly classifies a specified fraction of the labels
such that
$\operatorname{sign}\left(\langle{\bm{v}},{\bm{x}}_{i}\rangle-\theta\right)=y_{i}$.
We map a GNN learning problem to this task, as informally described below and
formally detailed in App. D.

######  Proposition 5 ($\mathsf{NP}$ hardness of GNN training; informal).

There exists a sequence of datasets indexed by the number of nodes $n$, each
of the form $\\{{\bm{A}}_{i},{\bm{x}}_{i},y_{i}\\}_{i=1}^{N}$ where
$y_{i}\in\\{-1,+1\\}$, such that determining whether a GNN with parameters
${\bm{a}},{\bm{b}}\in\mathbb{R}^{k}$, $c\in\mathbb{R}$ and $k=O(n)$ of the
form

$f_{{\bm{a}},{\bm{b}},c}({\bm{x}},{\bm{A}})=c+\sum\nolimits_{i=1}^{k}\bm{1}^{\top}\operatorname{relu}\left({\bm{x}}+a_{i}{\bm{A}}{\bm{x}}\right)b_{i}$
(15)

can correctly classify a specific fraction of the data is $\mathsf{NP}$ hard.

##  5 CSQ Lower bound for CNNs and Frame Averaging

Several commonly used invariant architectures take the form of symmetrized
neural networks over permutation subgroups $G\leq S_{n}$, which can be
generalized as frame-averaging architectures (Puny et al., 2023). This
includes translation invariant CNNs, group convolutional networks, and
networks with preprocessing steps such as sorting. We study the hardness of
learning these invariant functions.

#### Frame averaging and its connection with CNNs.

Let $G$ be a subgroup of the permutation group $S_{n}$ that acts on
$\mathbb{R}^{n\times d}$ by permuting the rows. A frame is a function
$\mathcal{F}:\mathbb{R}^{n\times d}\to 2^{G}\backslash\emptyset$ satisfying
$\mathcal{F}(g\cdot{\bm{X}})=g\mathcal{F}({\bm{X}})$ set-wise. As noted in
Puny et al. (2023), given some frame $\mathcal{F}$, one can transform an
arbitrary function (or neural network) $h({\bm{X}})$ into an invariant
function by averaging its values as
$\frac{1}{|\mathcal{F}({\bm{X}})|}\sum_{g\in\mathcal{F}({\bm{X}})}h(g^{-1}{\bm{X}})$.
For instance, setting $\forall{\bm{X}}:\mathcal{F}({\bm{X}})=G$ recovers the
Reynolds (or group-averaging) operator. Given a fixed frame $\mathcal{F}$, we
will first show CSQ lower bounds for the following class of frame-averaged
one-hidden-layer fully connected nets:

$\displaystyle\mathcal{H}_{\mathcal{F}}:=\Bigg{\\{}f:\mathbb{R}^{n\times
d}\to\mathbb{R}\text{,
}f({\bm{X}})=\frac{1}{\sqrt{|\mathcal{F}({\bm{X}})|}}\sum_{g\in\mathcal{F}({\bm{X}})}{\bm{a}}^{\top}\sigma$
$\displaystyle({\bm{W}}^{\top}(g^{-1}{\bm{X}}))\mathbf{1}_{d}\mid{\bm{W}}\in\mathbb{R}^{n\times
k},{\bm{a}}\in\mathbb{R}^{k}\Bigg{\\}},$

for some nonlinearity $\sigma$.

######  Example 2 (Frame for CNN).

Set $d=1$ and let $G$ be the cyclic group acting on $\mathbb{R}^{n}$ via
cyclic shifts of its elements with frame $\mathcal{F}({\bm{X}})=G,\forall
X\in\mathbb{R}^{n}$. Then, $\mathcal{F}^{d}_{n}$ consists of CNNs with one
convolutional layer and $k$ hidden channels.

######  Remark 6 (Difficulty of frames).

Since the frame $\mathcal{F}({\bm{X}})$ may vary by datapoint ${\bm{X}}$, the
distribution $\text{Unif}(\mathcal{F}({\bm{X}})^{-1})\circ{\bm{X}}$ can be
significantly different from the original distribution over ${\bm{X}}$, even
if ${\bm{X}}\sim g{\bm{X}}$ for all $g$. For example, consider the frame
$\mathcal{F}({\bm{X}})$ containing all permutations that sort ${\bm{X}}$
lexicographically (where $|\mathcal{F}({\bm{X}})|>1$ if ${\bm{X}}$ contains a
repeated row). Even if ${\bm{X}}\sim\mathcal{N}$ with $\mathcal{N}$ invariant
to row-permutation, the resultant distribution is not row permutation
invariant.

We focus on certain cases where such effects are not too pronounced. For
instance, it is simple to check that if $\mathcal{F}({\bm{X}})$ is constant
over ${\bm{X}}$, then $\mathcal{F}({\bm{X}})=G\>\forall{\bm{X}}$ (Lemma 28).
In the following, we assume that $G$ is a polynomial-sized (in $n$) subgroup
of $S_{n}$.

#### Family of hard functions.

Recall the low-dimensional function from the proof of Theorem 4:

$f^{*}_{\exp}:\mathbb{R}^{2\times d}\to\mathbb{R}\text{ with
}f({\bm{X}})=({\bm{a}}^{*})^{\top}\sigma(({\bm{W}}^{*})^{\top}{\bm{X}})\mathbf{1}_{d}$
(16)

for some special parameter ${\bm{a}}^{*}\in\mathbb{R}^{2k}$ and
${\bm{W}}^{*}\in\mathbb{R}^{2\times 2k}$ in Eq. 11. We now define the family
of hard functions, indexed by a set of matrices
$\mathcal{B}\subset\mathbb{R}^{n\times 2}$ obtained from Lemma 29:

$\displaystyle
C^{\mathcal{B}}_{\mathcal{F}}=\Big{\\{}g_{{\bm{B}}}:\mathbb{R}^{n\times
d}\to\mathbb{R}\text{ with }g_{{\bm{B}}}({\bm{X}})=\frac{\sum_{g\in
G}f^{*}_{\exp}({\bm{B}}^{\top}g^{-1}{\bm{X}})}{\sqrt{|G|\cdot\|f^{*}_{\exp}\|_{\mathcal{N}}}}\mid{\bm{B}}\in\mathcal{B}$
$\displaystyle\subset\mathbb{R}^{n\times 2}\Big{\\}}.$ (17)

######  Theorem 7 (Exponential CSQ lower bound for polynomial-sized group
averaging).

For feature dimension $d$ independent of inputs $n$ and width parameter
$k=\Theta(n)$, let $\epsilon>0$ be a sufficiently small error constant . Then
there exists a set $\mathcal{B}^{\prime}$ of size at least
$2^{\Omega(d^{\Omega(1)})}/|G|^{2}$ such that: for any target $g\in
C^{\mathcal{B}^{\prime}}_{\mathcal{F}}$ with $\|g\|_{\mathcal{N}}=1$, any CSQ
algorithm outputting a hypothesis $h$ with $\|g-h\|_{\mathcal{N}}\leq\epsilon$
requires either $2^{n^{\Omega(1)}}/|G|^{2}$ queries or at least one query with
precision $|G|2^{-n^{\Omega(1)}}+\sqrt{|G|}n^{-\Omega(k)}$.

###### Proof sketch.

Using a union bound argument over the original set $\mathcal{B}$ from
Diakonikolas et al. (2020), and concentration of inner product between
permuted unit vectors, we obtain from Lemma 29 a set of orthogonal matrices
$\mathcal{B}^{\prime}$ such that both $\|{\bm{B}}({\bm{B}}^{\prime})^{T}\|$
and $\|g{\bm{B}}(g^{\prime}{\bm{B}}^{\prime})^{T}\|$ are small
$\forall{\bm{B}}\neq{\bm{B}}^{\prime}\in\mathcal{B}^{\prime}$, $\forall g\neq
g^{\prime}\in G$. This construction costs a factor of $|G|^{2}$ in
$|\mathcal{B}^{\prime}|$. The rest of the proof proceeds similarly; see Sec.
F.1. ∎

We also obtain superpolynomial lower bounds for more general frames using the
technique of Goel et al. (2020). Here, we drop the assumption that $G$ is
polynomially-sized. The hard functions are based on parity functions, similar
to Goel et al. (2020) (but with an extra input dimension for ${\bm{X}}$):

$f_{S}:\mathbb{R}^{n\times d}\to\mathbb{R}\text{ with
}f_{S}({\bm{X}})=\mathbf{1}_{d}^{\top}\Big{(}\sum_{{\bm{w}}\in\\{-1,1\\}^{m}}\chi({\bm{w}})\sigma(\langle{\bm{w}},{\bm{X}}_{S}\rangle/\sqrt{m})\Big{)},\qquad
S\subseteq_{m}[n].$ (18)

where $\subseteq_{m}$ indicates a size $m$ subset, $\chi$ is the parity
function $\chi({\bm{w}})=\prod_{i}w_{i}$ and
${\bm{X}}_{S}\in\mathbb{R}^{m\times d}$ denotes entries of ${\bm{X}}$ in $S$.
Here,
$\langle{\bm{w}},{\bm{X}}_{S}\rangle:=\sum_{i=1}^{m}w_{i}({\bm{X}}_{S})_{i,:}$
and $\sigma:\mathbb{R}^{d}\to\mathbb{R}^{d}$ is an activation function. The
hard invariant function class consists of frame-averages of the original
functions above:

$\mathcal{H}_{\mathcal{S},\mathcal{F}}:=\left\\{\tilde{f}_{S}:\mathbb{R}^{n\times
d}\to\mathbb{R}\text{ with
}f({\bm{X}})=\frac{1}{\sqrt{|\mathcal{F}({\bm{X}})|}}\sum_{g\in\mathcal{F}({\bm{X}})}f_{S}(g^{-1}{\bm{X}})\mid
S\subset_{m}[n]\right\\}.$ (19)

Below, we show a hardness result for this class inspired by that of Goel et
al. (2020), whenever the frame $\mathcal{F}$ is sign invariant:
$\mathcal{F}({\bm{X}})=\mathcal{F}({\bm{X}}\circ{\bm{z}})$ for almost all
${\bm{X}}$ and for all ${\bm{z}}\in\\{-1,1\\}^{n}$.

######  Theorem 8 (Superpolynomial CSQ lower bound for sign-invariant frame
averaging).

If $\mathcal{F}$ is a sign-invariant frame, then for any $c>0$, any CSQ
algorithm that queries from an oracle of concept
$f\in\mathcal{H}_{\mathcal{S},\mathcal{F}}$ and output a hypothesis $h$ with
$\|f-h\|_{\mathcal{N}}\leq\Theta(1)$ needs at least $\Omega(n^{\Omega(\log
n)}/M(n))$ queries or at least one query with precision $O(n^{-c/2})$, where
$M(n)$ is the size of the largest orbit of size $m\leq\log(n)$ subsets of
$[n]$ under $G$. If moreover $|\mathcal{F}({\bm{X}})|=1$ for almost all
${\bm{X}}$ (a “singleton” frame), then the same result holds with $M(n)$
replaced by $1$.

###### Proof sketch..

The proof hinges on the structure of $f_{S}$, and in particular that applying
$g$ to the input of $f_{S}$ yields $f_{gS}$. The sign-invariance of
$\mathcal{F}$ is necessary to ensure that the proof of Goel et al. (2020),
which involves introducing random sign-flips, goes through. See Lemma 33 for
the full proof. ∎

######  Example 3 (Singleton frames).

An example of such a frame is the lexicographical sorting of
${\bm{X}}\in\mathbb{R}^{n\times d}$ according to the absolute value of the
entries (to make the frame sign-invariant). With probability $1$, ${\bm{X}}$
does not have any repeated values, and thus there is a unique lexicographical
sort for ${\bm{X}}$ almost surely. In fact, for any group $G$, a singleton
frame exists if ${\bm{X}}$ has no self-symmetry with probability 1, i.e.
${\bm{X}}\neq g{\bm{X}}$ for all $g\in G$.

For the hard class function to be nonvanishing, we show in Corollary 39 that
at least a superpolynomial-sized subset of
$\mathcal{H}_{\mathcal{S},\mathcal{F}}$ has norm lower-bounded by
$\Omega(\textup{poly}(n)^{-1})$, assuming that $|\mathcal{F}({\bm{X}})|$ is
polynomial in $n$ almost surely (note that this does not mean $|G|$ has to be
polynomial in $n$).

##  6 Separation between SQ and CSQ for invariant function classes

A natural question is whether there exist real-valued function classes which
are hard to learn in the CSQ setting (and hence by gradient descent with mean
squared error loss), but efficient to learn via a more carefully constructed
_non-CSQ_ algorithm. In the general (non-invariant) setting, one such
separation is that of learning $k$-sparse polynomials, i.e., degree $d$
polynomials in $n$ variables expressible as sums of only $k$ unique monomials
over inputs drawn from a product distribution
$\mathcal{D}=\bigtimes_{i=1}^{n}\mu_{i}$. Here, an efficient SQ (but not CSQ)
algorithm called the GROWING-BASIS algorithm is known from Andoni et al.
(2014). CSQ algorithms require $\Omega(n^{d})$ queries corresponding to the
number of degree $d$ orthogonal polynomials, but the GROWING-BASIS algorithms
learns with $O(nkf(d))$ SQ queries, where $f(d)$ may be exponential in $d$ but
independent of $n$. We informally illustrate the invariant extension of this
here, leaving complete details to App. G.

To perform this extension, we apply GROWING-BASIS to more general spaces (so-
called ‘rings’) of invariant polynomials whose basis consists of independent
invariant functions or ‘ _generators_ ’
$\\{g_{i}\in\mathbb{R}[{\bm{x}}]^{G}\\}_{i\in[r]}$, for some $r\in\mathbb{N}$.
Independence here asserts that any $g_{i}$ cannot be written as a polynomial
in the others. The space of all polynomials in $g_{i}({\bm{x}})$ with degree
$d$ and sparsity $k$ (both in $g_{i}$’s expansion) is denoted by
$\mathbb{R}[g_{1}({\bm{x}}),\ldots,g_{r}({\bm{x}})]_{d,k}$. We have the
following separation result:

######  Lemma 9 (informal).

Let ${\bm{x}}\sim\mathcal{D}$ such that the induced distribution of
$(g_{i}({\bm{x}}))_{i\in[r]}$ is a product distribution. Then any CSQ
algorithm that learns
$f\in\mathbb{R}[g_{1}({\bm{x}}),\ldots,g_{r}({\bm{x}})]_{d,k}$ of degree at
most $d=O(\log n)$ to a small constant error requires $\Omega(r^{d})$ CSQ
queries of bounded tolerance. However, the GROWING-BASIS algorithm learns $f$
with $O(kr^{2}\log r)$ SQ queries of precision
$\Omega(1/(k\operatorname{poly}(r)))$.

A more rigorous exposition and proof is prepared in App. G. Finding
independent generators itself can be a challenging task; however, we provide
an example below, and more in App. G, for separations over groups where
$r=\Omega(n)$ generators are known (Derksen and Kemper, 2015).

######  Example 4.

The sign group $G$ flips the sign of input elements and is commonly used to
study eigenvectors (Lim et al., 2022). $G$ is indexed by vectors
${\bm{v}}\in\\{-1,+1\\}^{n}$ with representation
$\rho({\bm{v}})\cdot{\bm{x}}={\bm{x}}\odot{\bm{v}}$ where $\odot$ represents
pointwise multiplication. Consider generators $g_{1},\dots,g_{n}$ where
$g_{i}({\bm{x}})=|{\bm{x}}_{i}|$. For inputs
${\bm{x}}\sim\operatorname{Unif}([-1,+1]^{n})$, the distribution of
$[g_{1}({\bm{x}}),\dots,g_{n}({\bm{x}})]$ is a product distribution
$\operatorname{Unif}([0,1]^{n})$ with the “shifted” Legendre polynomials as
the orthogonal polynomials. For $d$ at most $O(\log n)$, any CSQ algorithm
requires $\Omega(n^{d})$ queries of bounded tolerance to learn these
polynomials, whereas the GROWING-BASIS algorithm learns in $O(kn^{2}\log n)$
SQ queries.

##  7 Experiments

(a) GNN learning $\mathcal{H}_{ER,n}$

(b) CNN learning $C^{\mathcal{B}}_{\mathcal{F}}$

Figure 1: Overparameterized GNN (a) and CNN (b) fail to learn functions from
the class $\mathcal{H}_{ER,n}$ and $C^{\mathcal{B}}_{\mathcal{F}}$
respectively by either failing to fit the training set or overfitting the
data. Plots are aggregated and averaged over five random realizations.

We train overparameterized GNNs and CNNs on the hard functions from Sec. 4 and
Sec. 5, respectively. Specifically, we attempt to learn relatively small
instances drawn from the function classes $\mathcal{H}_{ER,n}$ for graphs
(Theorem 3) and $C^{\mathcal{B}}_{\mathcal{F}}$ for translationally invariant
data (Theorem 7). For the function class $\mathcal{H}_{ER,n}$, we consider
$n=15$ node inputs and the target function $g_{S,b}$ where $S\subseteq[n]$ is
a random subset of size $7$ and $b$ is either $0$ or $1$ at random. For the
function class $C^{\mathcal{B}}_{\mathcal{F}}$, we choose a random orthogonal
matrix ${\bm{B}}$ and then symmetrize the function class as per Eq. 17.

Fig. 1(a) and 1(b) plot the performance of the GNN and CNN respectively. The
GNN is unable to fit even the training data consisting of $225$ $n=15$ node
graphs drawn uniformly from the Erdős–Rényi model (i.e., $p=0.5$). The CNN
trained on $n=50$ dimension inputs fits the training set of size $500$
appropriately, but fails to generalize in the mean squared error loss. The
GNN/CNN were overparameterized with $3$ layers of graph/cyclic convolution
followed by a two layer ReLU MLP on the aggregated invariant features. We
refer the reader to App. H for further details.

##  8 Discussion

Our study asked how hard it is to learn relatively simple function classes
constrained by symmetry. Setting aside mathematical formality, one may
conjecture that such function classes are good representations or
approximations of the data observed in nature. Our results indicate that such
an assumption is unlikely to be a realistic one — at least in worst-case
settings, where we show such functions would be exponentially hard to learn.
Provable guarantees of learning will have to incorporate further assumptions
in the model class or biases in training to better account for the practical
success of learning algorithms (Lawrence et al., 2021; Gunasekar et al., 2018;
Le and Jegelka, 2022).

We now state potential directions for future work and conjectures. It would be
interesting to extend our hardness results beyond the SQ setting or to apply
to continuous groups. For example, a natural question is whether symmetric
functions are cryptographically hard to learn, as shown in Chen et al. (2022a)
for two hidden layer ReLU networks. A line of work in probability and
statistics has identified barriers to many algorithms solving average-case
settings of classic optimization problems, such as max cut or largest
independent set. These include the overlap gap property (Gamarnik, 2021) and
bounds on low degree algorithms (Hopkins and Steurer, 2017; Brennan et al.,
2020). Whether these barriers extend to the types of problems studied in
geometric deep learning is an open question. We should note that there are
restricted settings where learning neural networks is possible (see App. A).
Recent work (Chen et al., 2022b; Chen and Narayanan, 2023) has shown that ReLU
networks with $O(1)$ hidden nodes are learnable in polynomial time. Extending
this result to invariant network classes with $O(1)$ channels would similarly
expand the set of learnable invariant networks.

## Acknowledgements

The authors thank Sitan Chen for insightful discussions and feedback. BTK and
MW were supported by the Harvard Data Science Initiative Competitive Research
Fund and NSF award 2112085. TL and SJ were supported by NSF awards 2134108 and
CCF-2112665 (TILOS AI Institute), and Office of Naval Research grant
N00014-20-1-2023 (MURI ML-SCOPE). HL is supported by the Fannie and John Hertz
Foundation and the NSF Graduate Fellowship under Grant No. 1745302.

## References

  * Abbe and Sandon (2023) Emmanuel Abbe and Colin Sandon.  Polynomial-time universality and limitations of deep learning.  _Communications on Pure and Applied Mathematics_ , 76(11):3493–3549, 2023.
  * Andoni et al. (2014) Alexandr Andoni, Rina Panigrahy, Gregory Valiant, and Li Zhang.  Learning sparse polynomial functions.  In _Proceedings of the twenty-fifth annual ACM-SIAM symposium on Discrete algorithms_ , pages 500–510. SIAM, 2014.
  * Anschuetz and Kiani (2022) Eric R Anschuetz and Bobak T Kiani.  Quantum variational algorithms are swamped with traps.  _Nature Communications_ , 13(1):7760, 2022.
  * Anschuetz et al. (2022) Eric R Anschuetz, Andreas Bauer, Bobak T Kiani, and Seth Lloyd.  Efficient classical algorithms for simulating symmetric quantum systems.  _arXiv preprint arXiv:2211.16998_ , 2022.
  * Bakshi et al. (2019) Ainesh Bakshi, Rajesh Jayaram, and David P Woodruff.  Learning two layer rectified neural networks in polynomial time.  In _Conference on Learning Theory_ , pages 195–268. PMLR, 2019.
  * Barvinok (2005) A Barvinok.  Measure concentration lecture notes.  _See http://www. math. lsa. umich. edu/barvinok/total710. pdf_ , 2005\.
  * Batzner et al. (2022) Simon Batzner, Albert Musaelian, Lixin Sun, Mario Geiger, Jonathan P Mailoa, Mordechai Kornbluth, Nicola Molinari, Tess E Smidt, and Boris Kozinsky.  E (3)-equivariant graph neural networks for data-efficient and accurate interatomic potentials.  _Nature communications_ , 13(1):2453, 2022.
  * Behboodi et al. (2022) Arash Behboodi, Gabriele Cesa, and Taco S Cohen.  A pac-bayesian generalization bound for equivariant networks.  _Advances in Neural Information Processing Systems_ , 35:5654–5668, 2022.
  * Bietti et al. (2021) Alberto Bietti, Luca Venturi, and Joan Bruna.  On the sample complexity of learning under geometric stability.  In A. Beygelzimer, Y. Dauphin, P. Liang, and J. Wortman Vaughan, editors, _Advances in Neural Information Processing Systems_ , 2021.  URL https://openreview.net/forum?id=vlf0zTKa5Lh.
  * Blum and Rivest (1988) Avrim Blum and Ronald Rivest.  Training a 3-node neural network is np-complete.  _Advances in neural information processing systems_ , 1, 1988.
  * Blum et al. (1994) Avrim Blum, Merrick Furst, Jeffrey Jackson, Michael Kearns, Yishay Mansour, and Steven Rudich.  Weakly learning dnf and characterizing statistical query learning using fourier analysis.  In _Proceedings of the twenty-sixth annual ACM symposium on Theory of computing_ , pages 253–262, 1994.
  * Brennan et al. (2020) Matthew Brennan, Guy Bresler, Samuel B Hopkins, Jerry Li, and Tselil Schramm.  Statistical query algorithms and low-degree tests are almost equivalent.  _arXiv preprint arXiv:2009.06107_ , 2020.
  * Bronstein et al. (2021) Michael M Bronstein, Joan Bruna, Taco Cohen, and Petar Veličković.  Geometric deep learning: Grids, groups, graphs, geodesics, and gauges.  _arXiv preprint arXiv:2104.13478_ , 2021.
  * Brutzkus and Globerson (2017) Alon Brutzkus and Amir Globerson.  Globally optimal gradient descent for a convnet with gaussian inputs.  In _International conference on machine learning_ , pages 605–614. PMLR, 2017.
  * Cao and Gu (2019) Yuan Cao and Quanquan Gu.  Tight sample complexity of learning one-hidden-layer convolutional neural networks.  _Advances in Neural Information Processing Systems_ , 32, 2019.
  * Chandrasekaran et al. (2012) Venkat Chandrasekaran, Benjamin Recht, Pablo A. Parrilo, and Alan S. Willsky.  The convex geometry of linear inverse problems.  _Foundations of Computational Mathematics_ , 12(6):805–849, 2012.  doi: 10.1007/s10208-012-9135-7.  URL https://doi.org/10.1007/s10208-012-9135-7.
  * Chatterjee (2014) Sourav Chatterjee.  _Superconcentration and Related Topics_.  Springer Cham, 2014.  http://example.com/history_of_time.pdf(visited 2016-01-01).
  * Chen and Narayanan (2023) Sitan Chen and Shyam Narayanan.  A faster and simpler algorithm for learning shallow networks.  _arXiv preprint arXiv:2307.12496_ , 2023.
  * Chen et al. (2022a) Sitan Chen, Aravind Gollakota, Adam Klivans, and Raghu Meka.  Hardness of noise-free learning for two-hidden-layer neural networks.  _Advances in Neural Information Processing Systems_ , 35:10709–10724, 2022a.
  * Chen et al. (2022b) Sitan Chen, Adam R Klivans, and Raghu Meka.  Learning deep relu networks is fixed-parameter tractable.  In _2021 IEEE 62nd Annual Symposium on Foundations of Computer Science (FOCS)_ , pages 696–707. IEEE, 2022b.
  * Chen et al. (2023) Sitan Chen, Zehao Dou, Surbhi Goel, Adam Klivans, and Raghu Meka.  Learning narrow one-hidden-layer relu networks.  In _The Thirty Sixth Annual Conference on Learning Theory_ , pages 5580–5614. PMLR, 2023.
  * Chen et al. (2020) Zhengdao Chen, Lei Chen, Soledad Villar, and Joan Bruna.  Can graph neural networks count substructures?  _Advances in neural information processing systems_ , 33:10383–10395, 2020.
  * Chen and Zhu (2023) Ziyu Chen and Wei Zhu.  On the implicit bias of linear equivariant steerable networks: Margin, generalization, and their equivalence to data augmentation.  _arXiv preprint arXiv:2303.04198_ , 2023.
  * Chizat and Bach (2020) Lenaic Chizat and Francis Bach.  Implicit bias of gradient descent for wide two-layer neural networks trained with the logistic loss.  In _Conference on Learning Theory_ , pages 1305–1338. PMLR, 2020\.
  * Cong et al. (2019) Iris Cong, Soonwon Choi, and Mikhail D Lukin.  Quantum convolutional neural networks.  _Nature Physics_ , 15(12):1273–1278, 2019.
  * Daniely (2016) Amit Daniely.  Complexity theoretic limitations on learning halfspaces.  In _Proceedings of the forty-eighth annual ACM symposium on Theory of Computing_ , pages 105–117, 2016.
  * Daniely and Vardi (2020) Amit Daniely and Gal Vardi.  Hardness of learning neural networks with natural weights.  _Advances in Neural Information Processing Systems_ , 33:930–940, 2020.
  * Daniely et al. (2014) Amit Daniely, Nati Linial, and Shai Shalev-Shwartz.  From average case complexity to improper learning complexity.  In _Proceedings of the forty-sixth annual ACM symposium on Theory of computing_ , pages 441–448, 2014.
  * Daniely et al. (2023) Amit Daniely, Nathan Srebro, and Gal Vardi.  Most neural networks are almost learnable.  _arXiv preprint arXiv:2305.16508_ , 2023.
  * Das et al. (2019) Abhimanyu Das, Sreenivas Gollapudi, Ravi Kumar, and Rina Panigrahy.  On the learnability of deep random networks.  _arXiv preprint arXiv:1904.03866_ , 2019.
  * De Palma et al. (2019) Giacomo De Palma, Bobak Kiani, and Seth Lloyd.  Random deep neural networks are biased towards simple functions.  _Advances in Neural Information Processing Systems_ , 32, 2019.
  * Derksen and Kemper (2015) Harm Derksen and Gregor Kemper.  _Computational invariant theory_.  Springer, 2015.
  * Diakonikolas et al. (2017) Ilias Diakonikolas, Daniel M. Kane, and Alistair Stewart.  Statistical query lower bounds for robust estimation of high-dimensional gaussians and gaussian mixtures.  In _2017 IEEE 58th Annual Symposium on Foundations of Computer Science (FOCS)_ , pages 73–84, 2017.  doi: 10.1109/FOCS.2017.16.
  * Diakonikolas et al. (2020) Ilias Diakonikolas, Daniel M Kane, Vasilis Kontonis, and Nikos Zarifis.  Algorithms and sq lower bounds for pac learning one-hidden-layer relu networks.  In _Conference on Learning Theory_ , pages 1514–1539. PMLR, 2020\.
  * Du et al. (2018a) Simon Du, Jason Lee, Yuandong Tian, Aarti Singh, and Barnabas Poczos.  Gradient descent learns one-hidden-layer cnn: Don’t be afraid of spurious local minima.  In _International Conference on Machine Learning_ , pages 1339–1348. PMLR, 2018a.
  * Du and Goel (2018) Simon S Du and Surbhi Goel.  Improved learning of one-hidden-layer convolutional neural networks with overlaps.  _arXiv preprint arXiv:1805.07798_ , 2018.
  * Du et al. (2017) Simon S Du, Jason D Lee, and Yuandong Tian.  When is a convolutional filter easy to learn?  _arXiv preprint arXiv:1709.06129_ , 2017.
  * Du et al. (2018b) Simon S Du, Yining Wang, Xiyu Zhai, Sivaraman Balakrishnan, Russ R Salakhutdinov, and Aarti Singh.  How many samples are needed to estimate a convolutional neural network?  _Advances in Neural Information Processing Systems_ , 31, 2018b.
  * Dwivedi et al. (2020) Vijay Prakash Dwivedi, Chaitanya K Joshi, Anh Tuan Luu, Thomas Laurent, Yoshua Bengio, and Xavier Bresson.  Benchmarking graph neural networks.  _arXiv preprint arXiv:2003.00982_ , 2020.
  * Elesedy (2021) Bryn Elesedy.  Provably strict generalisation benefit for invariance in kernel methods.  _Advances in Neural Information Processing Systems_ , 34:17273–17283, 2021.
  * Feldman et al. (2012) Vitaly Feldman, Venkatesan Guruswami, Prasad Raghavendra, and Yi Wu.  Agnostic learning of monomials by halfspaces is hard.  _SIAM Journal on Computing_ , 41(6):1558–1590, 2012.
  * Fey and Lenssen (2019) Matthias Fey and Jan Eric Lenssen.  Fast graph representation learning with pytorch geometric.  _arXiv preprint arXiv:1903.02428_ , 2019.
  * Gamarnik (2021) David Gamarnik.  The overlap gap property: A topological barrier to optimizing over random structures.  _Proceedings of the National Academy of Sciences_ , 118(41):e2108492118, 2021.
  * Geerts and Reutter (2022) Floris Geerts and Juan L Reutter.  Expressiveness and approximation properties of graph neural networks.  _arXiv preprint arXiv:2204.04661_ , 2022.
  * Goel et al. (2018) Surbhi Goel, Adam Klivans, and Raghu Meka.  Learning one convolutional layer with overlapping patches.  In _International Conference on Machine Learning_ , pages 1783–1791. PMLR, 2018.
  * Goel et al. (2020) Surbhi Goel, Aravind Gollakota, Zhihan Jin, Sushrut Karmalkar, and Adam Klivans.  Superpolynomial lower bounds for learning one-layer neural networks using gradient descent.  In _International Conference on Machine Learning_ , pages 3587–3596. PMLR, 2020.
  * Gunasekar et al. (2018) Suriya Gunasekar, Jason D Lee, Daniel Soudry, and Nati Srebro.  Implicit bias of gradient descent on linear convolutional networks.  _Advances in neural information processing systems_ , 31, 2018.
  * Guruswami and Raghavendra (2009) Venkatesan Guruswami and Prasad Raghavendra.  Hardness of learning halfspaces with noise.  _SIAM Journal on Computing_ , 39(2):742–765, 2009\.
  * Hopkins and Steurer (2017) Samuel B Hopkins and David Steurer.  Bayesian estimation from few samples: community detection and related problems.  _arXiv preprint arXiv:1710.00264_ , 2017.
  * Jegelka (2022) Stefanie Jegelka.  Theory of graph neural networks: Representation and learning.  _arXiv preprint arXiv:2204.07697_ , 2022.
  * Ji and Telgarsky (2020) Ziwei Ji and Matus Telgarsky.  Directional convergence and alignment in deep learning.  _Advances in Neural Information Processing Systems_ , 33:17176–17186, 2020.
  * Judd (1987) J Stephen Judd.  Learning in networks is hard.  In _Proc. of 1st International Conference on Neural Networks, San Diego, California, June 1987_. IEEE, 1987.
  * Kalimeris et al. (2019) Dimitris Kalimeris, Gal Kaplun, Preetum Nakkiran, Benjamin Edelman, Tristan Yang, Boaz Barak, and Haofeng Zhang.  Sgd on neural networks learns functions of increasing complexity.  _Advances in neural information processing systems_ , 32, 2019.
  * Kearns (1998) Michael Kearns.  Efficient noise-tolerant learning from statistical queries.  _Journal of the ACM (JACM)_ , 45(6):983–1006, 1998.
  * Kibble (1945) WF Kibble.  An extension of a theorem of mehler’s on hermite polynomials.  In _Mathematical Proceedings of the Cambridge Philosophical Society_ , volume 41, pages 12–15. Cambridge University Press, 1945.
  * Klivans and Kothari (2014) Adam Klivans and Pravesh Kothari.  Embedding hard learning problems into gaussian space.  In _Approximation, Randomization, and Combinatorial Optimization. Algorithms and Techniques (APPROX/RANDOM 2014)_. Schloss Dagstuhl-Leibniz-Zentrum fuer Informatik, 2014.
  * Kolountzakis et al. (2005) Mihail N Kolountzakis, Evangelos Markakis, and Aranyak Mehta.  Learning symmetric k-juntas in time n^ o (k).  _arXiv preprint math/0504246_ , 2005.
  * Lawrence et al. (2021) Hannah Lawrence, Kristian Georgiev, Andrew Dienes, and Bobak T Kiani.  Implicit bias of linear equivariant networks.  _arXiv preprint arXiv:2110.06084_ , 2021.
  * Le and Jegelka (2022) Thien Le and Stefanie Jegelka.  Training invariances and the low-rank phenomenon: beyond linear networks.  In _International Conference on Learning Representations_ , 2022.  URL https://openreview.net/forum?id=XEW8CQgArno.
  * Li et al. (2021) Qunwei Li, Shaofeng Zou, and Wenliang Zhong.  Learning graph neural networks with approximate gradient descent.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 35, pages 8438–8446, 2021.
  * Liao et al. (2020) Renjie Liao, Raquel Urtasun, and Richard Zemel.  A pac-bayesian approach to generalization bounds for graph neural networks.  _arXiv preprint arXiv:2012.07690_ , 2020.
  * Lim et al. (2022) Derek Lim, Joshua Robinson, Lingxiao Zhao, Tess Smidt, Suvrit Sra, Haggai Maron, and Stefanie Jegelka.  Sign and basis invariant networks for spectral graph representation learning.  _arXiv preprint arXiv:2202.13013_ , 2022.
  * Lipton et al. (2005) Richard J Lipton, Evangelos Markakis, Aranyak Mehta, and Nisheeth K Vishnoi.  On the fourier spectrum of symmetric boolean functions with applications to learning symmetric juntas.  In _20th Annual IEEE Conference on Computational Complexity (CCC’05)_ , pages 112–119. IEEE, 2005.
  * Livni et al. (2014) Roi Livni, Shai Shalev-Shwartz, and Ohad Shamir.  On the computational efficiency of training neural networks.  _Advances in neural information processing systems_ , 27, 2014.
  * Long and Sedghi (2019) Philip M Long and Hanie Sedghi.  Generalization bounds for deep convolutional neural networks.  _arXiv preprint arXiv:1905.12600_ , 2019.
  * Loukas (2019) Andreas Loukas.  What graph neural networks cannot learn: depth vs width.  _arXiv preprint arXiv:1907.03199_ , 2019.
  * Lv (2021) Shaogao Lv.  Generalization bounds for graph convolutional neural networks via rademacher complexity.  _arXiv preprint arXiv:2102.10234_ , 2021.
  * Macdonald (1979) I G Macdonald.  _Symmetric functions and Hall polynomials / by I_.  Oxford University Press, Oxford : New York, 1979.
  * Mei et al. (2021) Song Mei, Theodor Misiakiewicz, and Andrea Montanari.  Learning with invariances in random features and kernel models.  In _Conference on Learning Theory_ , pages 3351–3418. PMLR, 2021\.
  * Meyer et al. (2023) Johannes Jakob Meyer, Marian Mularski, Elies Gil-Fuster, Antonio Anna Mele, Francesco Arzani, Alissa Wilms, and Jens Eisert.  Exploiting symmetry in variational quantum machine learning.  _PRX Quantum_ , 4(1):010328, 2023.
  * Mohri et al. (2018) Mehryar Mohri, Afshin Rostamizadeh, and Ameet Talwalkar.  _Foundations of machine learning_.  MIT press, 2018.
  * Morris et al. (2019) Christopher Morris, Martin Ritzert, Matthias Fey, William L Hamilton, Jan Eric Lenssen, Gaurav Rattan, and Martin Grohe.  Weisfeiler and leman go neural: Higher-order graph neural networks.  In _Proceedings of the AAAI conference on artificial intelligence_ , volume 33, pages 4602–4609, 2019.
  * Mossel et al. (2003) Elchanan Mossel, Ryan O’Donnell, and Rocco P Servedio.  Learning juntas.  In _Proceedings of the thirty-fifth annual ACM symposium on Theory of computing_ , pages 206–212, 2003.
  * Nguyen et al. (2022) Quynh T Nguyen, Louis Schatzki, Paolo Braccia, Michael Ragone, Patrick J Coles, Frederic Sauvage, Martin Larocca, and M Cerezo.  Theory for equivariant quantum neural networks.  _arXiv preprint arXiv:2210.08566_ , 2022.
  * Oymak and Soltanolkotabi (2018) Samet Oymak and Mahdi Soltanolkotabi.  End-to-end learning of a convolutional neural network via deep tensor decomposition.  _arXiv preprint arXiv:1805.06523_ , 2018.
  * Parberry (1994) Ian Parberry.  _Circuit complexity and neural networks_.  MIT press, 1994.
  * Paszke et al. (2019) Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  Pytorch: An imperative style, high-performance deep learning library.  In _Advances in Neural Information Processing Systems 32_ , pages 8024–8035. Curran Associates, Inc., 2019.
  * Pesah et al. (2021) Arthur Pesah, Marco Cerezo, Samson Wang, Tyler Volkoff, Andrew T Sornborger, and Patrick J Coles.  Absence of barren plateaus in quantum convolutional neural networks.  _Physical Review X_ , 11(4):041011, 2021.
  * Petrache and Trivedi (2023) Mircea Petrache and Shubhendu Trivedi.  Approximation-generalization trade-offs under (approximate) group equivariance.  _arXiv preprint arXiv:2305.17592_ , 2023.
  * Pisier (1986) Gilles Pisier.  Probabilistic methods in the geometry of banach spaces.  In Giorgio Letta and Maurizio Pratelli, editors, _Probability and Analysis_ , pages 167–241, Berlin, Heidelberg, 1986. Springer Berlin Heidelberg.  ISBN 978-3-540-40955-7.
  * Puny et al. (2023) Omri Puny, Derek Lim, Bobak T Kiani, Haggai Maron, and Yaron Lipman.  Equivariant polynomials for graph neural networks.  _arXiv preprint arXiv:2302.11556_ , 2023.
  * Ragone et al. (2022) Michael Ragone, Paolo Braccia, Quynh T Nguyen, Louis Schatzki, Patrick J Coles, Frederic Sauvage, Martin Larocca, and M Cerezo.  Representation theory for geometric quantum machine learning.  _arXiv preprint arXiv:2210.07980_ , 2022.
  * Reyzin (2020) Lev Reyzin.  Statistical queries and statistical algorithms: Foundations and applications.  _arXiv preprint arXiv:2004.00557_ , 2020.
  * Sannai et al. (2021) Akiyoshi Sannai, Masaaki Imaizumi, and Makoto Kawano.  Improved generalization bounds of group invariant/equivariant deep networks via quotient feature spaces.  In _Uncertainty in Artificial Intelligence_ , pages 771–780. PMLR, 2021.
  * Schatzki et al. (2022) Louis Schatzki, Martin Larocca, Frederic Sauvage, and Marco Cerezo.  Theoretical guarantees for permutation-equivariant quantum neural networks.  _arXiv preprint arXiv:2210.09974_ , 2022.
  * Sedghi et al. (2016) Hanie Sedghi, Majid Janzamin, and Anima Anandkumar.  Provable tensor methods for learning mixtures of generalized linear models.  In _Artificial Intelligence and Statistics_ , pages 1223–1231. PMLR, 2016.
  * Shah et al. (2020) Harshay Shah, Kaustav Tamuly, Aditi Raghunathan, Prateek Jain, and Praneeth Netrapalli.  The pitfalls of simplicity bias in neural networks.  _Advances in Neural Information Processing Systems_ , 33:9573–9585, 2020.
  * Shamir (2018) Ohad Shamir.  Distribution-specific hardness of learning neural networks.  _The Journal of Machine Learning Research_ , 19(1):1135–1163, 2018.
  * Skolik et al. (2023) Andrea Skolik, Michele Cattelan, Sheir Yarkoni, Thomas Bäck, and Vedran Dunjko.  Equivariant quantum circuits for learning on weighted graphs.  _npj Quantum Information_ , 9(1):47, 2023.
  * Sokolic et al. (2017) Jure Sokolic, Raja Giryes, Guillermo Sapiro, and Miguel Rodrigues.  Generalization error of invariant classifiers.  In _Artificial Intelligence and Statistics_ , pages 1094–1103. PMLR, 2017.
  * Song et al. (2017) Le Song, Santosh Vempala, John Wilmes, and Bo Xie.  On the complexity of learning neural networks.  _Advances in neural information processing systems_ , 30, 2017.
  * Song et al. (2021) Min Jae Song, Ilias Zadik, and Joan Bruna.  On the cryptographic hardness of learning single periodic neurons.  _Advances in neural information processing systems_ , 34:29602–29615, 2021.
  * Spencer (2014) Joel Spencer.  _Asymptopia_ , volume 71.  American Mathematical Soc., 2014.
  * Szörényi (2009) Balázs Szörényi.  Characterizing statistical query learning: simplified notions and proofs.  In _International Conference on Algorithmic Learning Theory_ , pages 186–200. Springer, 2009.
  * Tahmasebi and Jegelka (2023) Behrooz Tahmasebi and Stefanie Jegelka.  The exact sample complexity gain from invariances for kernel regression on manifolds.  _arXiv preprint arXiv:2303.14269_ , 2023.
  * Talagrand (1994) Michel Talagrand.  On Russo’s Approximate Zero-One Law.  _The Annals of Probability_ , 22(3):1576 – 1587, 1994.  doi: 10.1214/aop/1176988612.  URL https://doi.org/10.1214/aop/1176988612.
  * Valiant (1984) Leslie G Valiant.  A theory of the learnable.  _Communications of the ACM_ , 27(11):1134–1142, 1984.
  * Valle-Perez et al. (2018) Guillermo Valle-Perez, Chico Q Camargo, and Ard A Louis.  Deep learning generalizes because the parameter-function map is biased towards simple functions.  _arXiv preprint arXiv:1805.08522_ , 2018.
  * Vardi (2023) Gal Vardi.  On the implicit bias in deep-learning algorithms.  _Communications of the ACM_ , 66(6):86–93, 2023\.
  * Vempala and Wilmes (2019) Santosh Vempala and John Wilmes.  Gradient descent for one-hidden-layer neural networks: Polynomial convergence and sq lower bounds.  In _Conference on Learning Theory_ , pages 3115–3117. PMLR, 2019\.
  * Verma and Zhang (2019) Saurabh Verma and Zhi-Li Zhang.  Stability and generalization of graph convolutional neural networks.  In _Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_, pages 1539–1548, 2019.
  * Wang et al. (2021) Dian Wang, Robin Walters, and Robert Platt.  So(2)-equivariant reinforcement learning.  In _International Conference on Learning Representations_ , 2021.
  * Wu et al. (2020) Zonghan Wu, Shirui Pan, Fengwen Chen, Guodong Long, Chengqi Zhang, and S Yu Philip.  A comprehensive survey on graph neural networks.  _IEEE transactions on neural networks and learning systems_ , 32(1):4–24, 2020.
  * Xu et al. (2018) Keyulu Xu, Weihua Hu, Jure Leskovec, and Stefanie Jegelka.  How powerful are graph neural networks?  _arXiv preprint arXiv:1810.00826_ , 2018.
  * Yun et al. (2020) Chulhee Yun, Shankar Krishnan, and Hossein Mobahi.  A unifying view on implicit bias in training linear neural networks.  _arXiv preprint arXiv:2010.02501_ , 2020.
  * Zhang et al. (2023) Bohang Zhang, Shengjie Luo, Liwei Wang, and Di He.  Rethinking the expressive power of gnns via graph biconnectivity.  _arXiv preprint arXiv:2301.09505_ , 2023.
  * Zhang et al. (2020) Shuai Zhang, Meng Wang, Sijia Liu, Pin-Yu Chen, and Jinjun Xiong.  Fast learning of graph neural networks with guaranteed generalizability: one-hidden-layer case.  In _International Conference on Machine Learning_ , pages 11268–11277. PMLR, 2020.
  * Zhang et al. (2017) Yuchen Zhang, Jason Lee, Martin Wainwright, and Michael I Jordan.  On the learnability of fully-connected neural networks.  In _Artificial Intelligence and Statistics_ , pages 83–91. PMLR, 2017\.
  * Zhong et al. (2017) Kai Zhong, Zhao Song, and Inderjit S Dhillon.  Learning non-overlapping convolutional neural networks with multiple kernels.  _arXiv preprint arXiv:1711.03440_ , 2017.
  * Zhou and Feng (2018) Pan Zhou and Jiashi Feng.  Understanding generalization and optimization performance of deep cnns.  In _International Conference on Machine Learning_ , pages 5960–5969. PMLR, 2018.

## Table of Contents

[sections] [sections]l1

##  Appendix A Extended related works

#### Restricted learning algorithms for CNNs and GNNs

Various works show that convolutional and graph neural networks can be
efficiently learned under various assumptions on the form of the network. For
given input distributions, Brutzkus and Globerson [2017], Du et al. [2017]
give an efficient learning algorithm for learning a single convolution filter
and Zhong et al. [2017] give a polynomial time algorithm for learning a sum of
a polynomial number of convolution kernels applied to an input. These can be
viewed as a single hidden layer network where the last layer weights are all
equal to one. Du and Goel [2018] (improving on Goel et al. [2018]) give a
polynomial time algorithm for learning single layer convolutional networks
(not necessarily invariant) with filters that have stride at least half the
filter width. These CNNs have a single channel but are not invariant to the
cyclic group due to the presence of weights in the second layer that vary by
each patch of the filter. Similarly, Du et al. [2018a] show that gradient
descent can also learn networks of a similar form with non-overlapping patches
in the random setting (weights are random) with high probability. Oymak and
Soltanolkotabi [2018] learn convolutional networks of depth $D$ layers with a
single kernel per layer and make a similar assumption that stride is at least
as large as width. Zhang et al. [2020], Li et al. [2021] provide learning
algorithms for a class of one hidden layer GNNs via gradient descent based
algorithms whose sample complexity and runtime depend on factors like the
condition number or norm of the weights of the target function. These
algorithms are only guaranteed to be polynomial in runtime for the restricted
instances when such factors are appropriately bounded.

#### Learning feedforward neural networks

The study of the hardness of learning feedforward nueral networks has a rich
history dating back many decades. Judd [1987], Blum and Rivest [1988] show
that in the proper learning setting, it is $\mathsf{NP}$ complete to find a
set of weights of a given network that fits a given training set. These
results were later expanded in [Zhang et al., 2017] to more realistic
settings. Many works study the hardness of improperly learning the class of
ReLU feedforward networks. In the SQ setting, Goel et al. [2018], Diakonikolas
et al. [2020], Song et al. [2017] show that at least a superpolynomial number
of queries are needed to learn even single hidden layer networks in the
correlational statistical query model. Chen et al. [2022a] show that networks
with two hidden layers are even hard in the general SQ setting by using the
networks to round inputs to boolean inputs and applying hardness results from
Boolean learning theory and cryptography. We should remark that there are also
other papers that reduce the task of learning feedforward neural networks to
average case or cryptographically hard problems [Song et al., 2021, Daniely
and Vardi, 2020]. To sidestep these hardness results and provide proofs of
learnability in the PAC setting, a number of works make assumptions on the
networks or inputs/outputs to give efficient PAC learning algorithms for
single hidden layer neural networks [Bakshi et al., 2019, Goel et al., 2018,
Shamir, 2018, Sedghi et al., 2016, Vempala and Wilmes, 2019]. Such assumptions
include bounds on the condition number [Bakshi et al., 2019], approximation
guarantees of the network by polynomials [Vempala and Wilmes, 2019],
positivity of the second layer weights Diakonikolas et al. [2020], and others.
In perhaps the most general setting, Chen and Narayanan [2023], Chen et al.
[2023] give a polynomial time learning algorithm for learning single hidden
layer ReLU networks with $O(1)$ hidden nodes (i.e. runtime is exponential in
the number of nodes but polynomial in other parameters such as the input size
and error).

#### Random neural network learnability

Though the class of neural networks may be challenging to learn, various works
study the average case setting where one is interested in learning random
neural networks or random features from neural networks. In fact, prior work
has shown that the complexity of such random neural networks, measured via
statistical measures, continuity, or robustness, is rather low though this is
by no means a guarantee of learnability [Kalimeris et al., 2019, De Palma et
al., 2019, Shah et al., 2020, Valle-Perez et al., 2018]. Various hardness
results exist for this setting of learning random neural networks. Das et al.
[2019] gives SQ hardness results for learning neural networks with sign
activation that shows query complexity increasing exponentially with the depth
of the network. Daniely and Vardi [2020] show that neural networks with random
weights are hard to learn if the input distribution is allowed to depend on
the weights based on reductions to random K-SAT. Nevertheless, recent work by
Daniely et al. [2023] shows that random constant depth neural networks with
ReLU activation can be learned in error $\epsilon$ in time scaling as
$d^{\operatorname{polylog}(1/\epsilon)}$ where $d$ is the size of the network.

#### Sample complexity and generalization bounds

Separate from the question of computational complexity of learning is the
question of how many samples are needed to learn a function class. Specific to
(exactly or approximately) invariant neural network architectures, there exist
papers studying the sample complexity of learning convolutional neural
networks [Du et al., 2018b, Zhou and Feng, 2018, Cao and Gu, 2019] and graph
neural networks [Zhang et al., 2020].

Various papers consider how generalization bounds improve when enforcing
equivariance. These include generalization bounds based on covering numbers or
the complexity of the invariant function space Sokolic et al. [2017], Petrache
and Trivedi [2023], based on invariant kernel method algorithms Bietti et al.
[2021], Tahmasebi and Jegelka [2023], Elesedy [2021], and based on equivariant
versions of norm-based PAC-Bayesian generalization bounds [Behboodi et al.,
2022]. For the group of translations, there also exist papers studying
generalization bounds for architectures with convolutional layers though these
architectures are not strictly invariant [Long and Sedghi, 2019, Zhou and
Feng, 2018]. Similarly, there are various generalization bounds as well for
the class of graph neural networks [Lv, 2021, Verma and Zhang, 2019].

#### Other related topics

Though not directly a statement of computational hardness, various papers give
no-go theorems for learning function classes by studying expressivity
limitations of graph neural networks [Jegelka, 2022, Wu et al., 2020, Loukas,
2019]. These expressivity results provide a set of functions that a graph
neural network cannot express and thus by definition also cannot learn to
arbitrary accuracy. The most well known set of results are related to
limitations of graph neural networks in distinguishing non-isomorphic graphs
via the Weisfeiler-Lehman hierarchy [Xu et al., 2018, Morris et al., 2019,
Geerts and Reutter, 2022]. Other work show expressivity limitations of GNNs by
studying their power in expressing invariant polynomials [Puny et al., 2023],
identifying graph biconnectivity [Zhang et al., 2023], counting substructures
[Chen et al., 2020], and through various other means.

Various works study the implicit bias of neural networks primarily to help
address why neural networks can even learn in overparameterized settings
[Vardi, 2023, Chizat and Bach, 2020, Ji and Telgarsky, 2020]. Here, the aim of
these works is to show that gradient descent will converge to specific
regularized functions that fit a given training set. Implicit bias results
have been proven for linear but multi-layer classes of convolutional neural
networks [Gunasekar et al., 2018, Yun et al., 2020, Le and Jegelka, 2022] and
equivariant neural networks [Lawrence et al., 2021, Chen and Zhu, 2023].
Generally, these results show that gradient descent on the parameter space of
such neural networks is implicitly regularized under some norm or semi-norm
depending on properties of the group. Extending these results to more
realistic settings is likely challenging.

In Boolean learning theory, there are various results on the problem of
learning a symmetric $k$-junta [Mossel et al., 2003, Lipton et al., 2005,
Kolountzakis et al., 2005]. Here, we are promised that the Boolean function
only depends on $k$ variables, and within that $k$ variable subset, the
function is symmetric and thus only depends on the number of $1$s in that
subset. For this problem, Kolountzakis et al. [2005] achieve a runtime
$n^{O(k\log k)}$ that is polynomial in $n$ for fixed $k$. Note, that the
functions we study in the Boolean setting are assumed to be symmetric over the
whole input space and not just on a size $k$ subset of the inputs bits which
drastically simplifies the problem. In fact, the hardness of the symmetric
$k$-junta learning problem largely arises from challenges in finding which $k$
bits are in the support of the Boolean function.

Finally, in quantum machine learning, there has been a recent interest in
studying quantum variational algorithms or quantum neural networks that are
constrained under symmetries [Meyer et al., 2023, Ragone et al., 2022, Skolik
et al., 2023, Cong et al., 2019]. Some recent theoretical work has quantified
the hardness and sample complexity associated with these models that are
generally linear though acting on high dimensional Hilbert spaces [Pesah et
al., 2021, Schatzki et al., 2022, Anschuetz and Kiani, 2022, Nguyen et al.,
2022, Anschuetz et al., 2022]. Since these quantum models typically operate in
a different regime than classical models and are typically focused on quantum
tasks, they are out of the scope of our current work.

##  Appendix B Extended background into learning frameworks

There are various models used to rigorously quantify the hardness of learning.
We recommend the review in Reyzin [2020] for an overview of the SQ formalism
and its connections to other learning models. Here, we briefly mention the PAC
model and overview further the background into the SQ framework from the main
text.

Perhaps the most widely used framework for quantifying hardness of learning is
the provably approximately correct (PAC) model [Valiant, 1984].

######  Definition 10 (PAC Learning [Mohri et al., 2018]).

A concept class $\mathcal{C}$ is PAC-learnable if there exists an algorithm
such that for any $\epsilon>0$ and $\delta>0$, for all distributions
$\mathcal{D}$ on $\mathcal{X}$, and for any target concept
$c^{*}\in\mathcal{C}$, the algorithm takes at most
$m=\operatorname{poly}(1/\epsilon,1/\delta,n)$ samples drawn i.i.d. from
$({\bm{x}},c({\bm{x}}))$ with ${\bm{x}}\sim\mathcal{D}$, and returns a
function $f$ satisfying $\|f-c^{*}\|_{\mathcal{D}}\leq\epsilon$ with
probability at most $1-\delta$. If the algorithm runs in time at most
$\operatorname{poly}(1/\epsilon,1/\delta,n)$, then $\mathcal{C}$ is
_efficiently PAC learnable_. The algorithm is a _proper learner_ if it returns
$f\in\mathcal{C}$, and otherwise denoted an _improper learner_.

Proving computational hardness in the PAC model typically requires reducing
learning tasks to cryptographically or complexity theoretically hard problems.
Unconditional hardness results for learning with neural networks would imply
$\mathsf{P}\neq\mathsf{NP}$: Abbe and Sandon [2023] show that the task of
training neural networks via gradient descent is $\mathsf{P}$-Complete, i.e.
any poly-time algorithm in $\mathsf{P}$ can be reduced to the task of training
a given poly-sized neural network with stochastic gradient descent.

Given the above, the statistical query formalism sidesteps these challenges by
constraining algorithms to be composed of a set of noisy queries. Perhaps the
most important restriction in the SQ framework is that algorithms do not have
access to individual data points but only noisy queries averaged across the
input/output distribution. The classic example of an efficient non-SQ
algorithm is Gaussian elimination for learning parities [Blum et al., 1994].
Nevertheless, the SQ formalism naturally captures most algorithms used in
practice including gradient descent as shown for example in Example 1 of the
main text for gradients over the mean squared error loss.

For boolean functions, correlational statistical queries can capture any
general statistical query since any query can be split into its $+1$ or $-1$
output cases [Reyzin, 2020]. For example, assume one performs an SQ query with
query function $q:\mathcal{X}\times\\{-1,+1\\}\to[-1,+1]$, then

$\mathbb{E}_{({\bm{x}},y)\sim\mathcal{D}}\left[g({\bm{x}},y)\right]=\mathbb{E}\left[g({\bm{x}},+1)\frac{1+y}{2}+g({\bm{x}},-1)\frac{1-y}{2}\right],$
(20)

which decomposes into two correlational statistical queries and two functions
independent of the target.

We remarked in the main text that generalizing real-valued hardness results
proven in the CSQ framework to the SQ framework is in general challenging. In
fact, Vempala and Wilmes [2019] more rigorously underscored this challenge
(see their Proposition 4.1) by providing a rather unnatural, yet efficient, SQ
algorithm that learns any finite real-valued function class satisfying a non-
degeneracy assumption that the set of points $f(x)=g(x)$ is measure zero for
any pair of functions $f,g$ in the class. Namely, they show the following no-
go theorem.

######  Theorem 11 (Limitations on $\operatorname{SQ}$ hardness for real-
valued functions [Vempala and Wilmes, 2019]).

Given a concept class $\mathcal{C}$ of functions $f:\mathcal{X}\to\mathbb{R}$
of size $|\mathcal{C}|$ where for all $f,g\in\mathcal{C}$ the set of inputs
where $f(x)=g(x)$ is measure zero, the family $\mathcal{C}$ can be learned
with $O(\log|\mathcal{C}|)$ queries to $\operatorname{SQ}$ with a constant
error tolerance $\tau$.

###### Proof.

We aim to choose a query $g:\mathcal{C}\to[-1,1]$ which maps the target
function $f\in\mathcal{C}$ to a unique value in the range $[-1,1]$. For our
target function $f\in\mathcal{C}$, let us call this value $v$. To specify the
form of this query, it suffices to send each pair $(x,f(x))$ to its
corresponding unique value in $[-1,+1]$, i.e. output $v$ for every input
$(x,f(x))$. Since outputs of disjoint functions in $\mathcal{C}$ are equal
only on a support of measure zero, this mapping can be performed with
probability one. Thus, we can query the value $v$ and receive an output, say
$v^{\prime}$ up to error $\tau$. Repeating this procedure iteratively to only
include functions in the tolerance range $[v^{\prime}-\tau,v^{\prime}+\tau]$,
we can select the target function after $O(\log|\mathcal{C}|)$ queries. ∎

Loosely, the above theorem shows that real-valued classes of functions are
only hard in SQ settings when there exist sufficiently many hard functions
with finitely many outputs. Chen et al. [2022a] achieved such a result for
two-hidden layer networks by constructing networks that round inputs to the
nearest boolean bitstring.

##  Appendix C Boolean statistical query setting

As discussed in Sec. 3, the Boolean function setting provides a nice
illustration of the types of results obtainable via SQ lower bounds. For
further background into the boolean SQ setting, we recommend the review in
Reyzin [2020]. Here, we will prove Theorem 2 restated below for convenience.

See 2

For a simple proof of the above statement, we will follow the proof technique
in Appendix B of Chen et al. [2022a]. The above statement will follow as a
consequence of the pairwise independence of functions in $\mathcal{H}_{\rho}$
defined below.

######  Definition 12 (Boolean Case of Definition C.1 of Chen et al. [2022a]).

Let $\mathcal{C}$ be a hypothesis class mapping $\\{-1,+1\\}^{n}$ to
$\\{-1,+1\\}$, and let $\mathcal{D}$ be a distribution on $\mathcal{X}$.
$\mathcal{C}$ is a $(1-\eta)$-pairwise independent function family if with
probability $1-\eta$ over the choice of ${\bm{x}},{\bm{x}}^{\prime}$ drawn
independently from $\mathcal{D}$, the distribution of
$(f({\bm{x}}),f({\bm{x}}^{\prime}))$ for $f$ drawn uniformly at random from
$\mathcal{C}$ is the product distribution
$\operatorname{Unif}(\\{-1,+1\\}^{2})$.

Pairwise independent function families offer a convenient means to bound the
number of queries needed to distinguish the outputs of a given unknown
function $f^{*}$ versus uniformly random outputs from the function class.

######  Theorem 13 (Theorem C.4 of Chen et al. [2022a]).

Let function class $\mathcal{C}$ mapping $\\{-1,+1\\}^{n}$ to $\\{-1,+1\\}$ be
a $(1-\eta)$-pairwise independent function family w.r.t. a distribution
$\mathcal{D}$ on $\\{-1,+1\\}^{n}$. For any $f\in\mathcal{C}$, let
$\mathcal{D}_{f}$ denote the distribution of $({\bm{x}},f({\bm{x}}))$ where
${\bm{x}}\sim\mathcal{D}$. Let
$\mathcal{D}_{\operatorname{Unif}(\mathcal{C})}$ denote the distribution of
$({\bm{x}},y)$ where ${\bm{x}}\sim\mathcal{D}$ and $y=f({\bm{x}})$ for
$f\sim\operatorname{Unif}(\mathcal{C})$. Any SQ learner able to distinguish
the labeled distribution $\mathcal{D}_{f^{*}}$ for an unknown
$f^{*}\in\mathcal{C}$ from the randomly labeled distribution
$\mathcal{D}_{\operatorname{Unif}(\mathcal{C})}$ using bounded queries of
tolerance $\tau$ requires at least $\frac{\tau^{2}}{2\eta}$ such queries.

Finally, we are ready to prove Theorem 2.

###### Proof of Theorem 2.

Since the symmetric function is constant on orbits, we can represent any
symmetric function as $f:\mathcal{O}_{n}\to\\{-1,+1\\}$ and the correlation
between two functions $f,g$ can be written as

$\langle
f,g\rangle_{SQ}=\mathbb{E}\left[fg\right]=2^{-n}\sum_{O_{k}\in\mathcal{O}_{\rho}}|O_{k}|f({\bm{x}}_{O_{k}})g({\bm{x}}_{O_{k}}).$
(21)

The probability $p_{\mathcal{O}_{\rho}}:\mathcal{O}_{\rho}\to[0,1]$ of a
uniformly random bitstring falling in orbit $O_{K}\in\mathcal{O}_{\rho}$ is
equal to $p_{\mathcal{O}_{\rho}}(O_{k})=|O_{k}|/2^{n}$. Note that the
distribution of $(f({\bm{x}}_{O}),f({\bm{x}}_{O^{\prime}}))$ is equal to the
uniform distribution whenever ${\bm{x}}_{O}\neq{\bm{x}}_{O^{\prime}}$ since we
have a uniform distribution over all possible symmetric functions. This occurs
with probability

$\eta=\sum_{O_{k}\in\mathcal{O}_{\rho}}\left(\frac{|O_{k}|}{2^{n}}\right)^{2}=\|p_{\mathcal{O}_{\rho}}\|^{2}.$
(22)

Therefore, $\mathcal{H}_{\rho}$ is a $(1-\eta)$-pairwise independent function
family with parameter $\eta$ as above.

To apply Theorem 13, note that for a given function
$f^{*}\in\mathcal{H}_{\rho}$, to distinguish $\mathcal{D}_{f^{*}}$ from
$\mathcal{D}_{\operatorname{Unif}(\mathcal{H}_{\rho})}$, it suffices to return
a function $f$ that has classification error at most $1/4$ with respect to
$f^{*}$. Since distinguishing the two distributions requires
$\frac{\tau^{2}}{2\eta}$ queries, we arrive at the final result. ∎

##  Appendix D Complexity theoretic hardness extension

In this section, we would like to show a reduction from the task of training
an invariant neural network on a given dataset to a previously $\mathsf{NP}$
hard problem. This extends the classic result of Blum and Rivest [1988]
showing that there exists a sequence of datasets such that determining whether
or not a 3-node neural network can fit that dataset is hard. Though we
consider extending this to a simple graph neural network here for simplicity,
we remark that we chose this only for convenience and similar extensions can
be found for other classes of networks.

The particular $\mathsf{NP}$ hard problem we consider is the decision version
of the learning halfspaces with noise problem [Guruswami and Raghavendra,
2009, Feldman et al., 2012].

######  Problem 1 (Learning halfspaces with noise [Guruswami and Raghavendra,
2009]).

Given a set of $N$ labeled examples $\\{{\bm{x}}_{i},y_{i}\\}_{i=1}^{N}$ where
${\bm{x}}_{i}\in\\{0,1\\}^{n}$ and $y_{i}\in\\{-1,+1\\}$ and constants
$\delta,\epsilon$ where $0<\delta\leq\epsilon<1$, we distinguish between the
following two cases:

  * •

Case 1: There exists a vector ${\bm{v}}\in\mathbb{R}^{n}$, constant
$\theta\in\mathbb{R}$, and set $S\subseteq[N]$ of size at least
$|S|\geq(1-\delta)N$ such that
$y_{i}\left(\langle{\bm{v}},{\bm{x}}_{i}\rangle-\theta\right)\geq 0$ for all
$i\in S$.

  * •

Case 2: For all ${\bm{v}}\in\mathbb{R}^{n}$, $\theta\in\mathbb{R}$, and sets
$S\subseteq[N]$ such that $|S|\geq(1-\epsilon)N$, there exists some $i\in S$
such that $y_{i}\left(\langle{\bm{v}},{\bm{x}}_{i}\rangle-\theta\right)<0$.

Consider the following class of message passing GNNs with a single hidden
layer and $k$ channels which acts on an adjacency matrix
${\bm{A}}\in\\{0,1\\}^{n\times n}$ and node features
${\bm{x}}\in\mathbb{R}^{n}$:

$f_{{\bm{a}},{\bm{b}},c}({\bm{x}},{\bm{A}})=c+\sum_{i=1}^{k}\bm{1}^{\top}\operatorname{relu}\left({\bm{x}}+a_{i}{\bm{A}}{\bm{x}}\right)b_{i}.$
(23)

We will map the below problem of fitting the weights ${\bm{a}},{\bm{b}},c$ of
the GNNs above to a given dataset consisting of a set of graphs and their
corresponding outputs.

######  Problem 2.

Given a labeled set of graphs, node features, and outputs
$\\{{\bm{A}}_{i},{\bm{x}}_{i},y_{i}\\}_{i=1}^{N}$ where
${\bm{A}}_{i}\in\\{0,1\\}^{n\times n}$, ${\bm{x}}_{i}\in\mathbb{R}^{n}$ and
$y_{i}\in\\{-1,+1\\}$ and constants $\delta,\epsilon$ where
$0<\delta\leq\epsilon<1$, distinguish between the following two cases:

  * •

Case 1: There exists weights
$a_{1},\dots,a_{k},b_{1},\dots,b_{k},c\in\mathbb{R}$ and set $S\subseteq[N]$
of size at least $|S|\geq(1-\delta)N$ such that
$y_{i}f_{{\bm{a}},{\bm{b}},c}({\bm{x}}_{i},{\bm{A}}_{i})\geq 0$ for all $i\in
S$.

  * •

Case 2: For all weights $a_{1},\dots,a_{k},b_{1},\dots,b_{k},c\in\mathbb{R}$,
and sets $S\subseteq[N]$ such that $|S|\geq(1-\epsilon)N$, there exists some
$i\in S$ such that
$y_{i}f_{{\bm{a}},{\bm{b}},c}({\bm{x}}_{i},{\bm{A}}_{i})<0$.

######  Proposition 14 ($\mathsf{NP}$ hardness of GNN training).

Training the GNN in Eq. 23 to solve 2 is $\mathsf{NP}$ hard.

###### Proof.

We show that solving 2 is equivalent to solving 1 which is previously shown to
be $\mathsf{NP}$ hard [Guruswami and Raghavendra, 2009, Feldman et al., 2012].
Given a dataset $\\{{\bm{x}}_{i},y_{i}\\}_{i=1}^{N}$ for 1 over $n$-bit
Boolean strings, we map this to a dataset over graphs
$\\{{\bm{A}}_{i},{\bm{x}}_{i}^{\prime},y_{i}^{\prime}\\}_{i=1}^{N}$ of
$\sum_{k=1}^{n}k=n(n+1)/2$ nodes where ${\bm{x}}_{i}^{\prime}=-\bm{1}$ for all
inputs and $y_{i}=y_{i}^{\prime}$ are the same for both datasets. Adjacency
matrix ${\bm{A}}_{i}$ is constructed by including a disconnected $k$-clique
(including self loops) in the graph whenever $[{\bm{x}}_{i}]_{k}=1$. Since
there are $n(n+1)/2$ nodes, there are sufficient nodes to construct any such
graph. E.g., for ${\bm{x}}=(1,0,1)^{\top}$, a corresponding adjacency matrix
is

${\bm{A}}=\begin{bmatrix}1&0&0&0&0&0\\\ 0&0&0&0&0&0\\\ 0&0&0&0&0&0\\\
0&0&0&1&1&1\\\ 0&0&0&1&1&1\\\ 0&0&0&1&1&1\\\ \end{bmatrix}.$ (24)

We also assume that the GNNs have $k=n$ channels. In this setting, we have
that for any $i\in[N]$

$\begin{split}f_{{\bm{a}},{\bm{b}},c}({\bm{x}}_{i}^{\prime},{\bm{A}}_{i})&=c+\sum_{j=1}^{n}\bm{1}^{\top}\operatorname{relu}\left({\bm{x}}_{i}^{\prime}+a_{j}{\bm{A}}_{i}{\bm{x}}_{i}^{\prime}\right)b_{j}\\\
&=c+\sum_{j=1}^{n}b_{j}\sum_{\ell=1}^{n}\ell\mathbb{1}[\exists\;\ell-\text{clique}\in{\bm{A}}_{i}]\operatorname{relu}(-1-\ell
a_{j}),\end{split}$ (25)

where $\mathbb{1}[\exists\;\ell-\text{clique}\in{\bm{A}}_{i}]$ is the
indicator function for whether or not a $\ell$-clique is in the graph
${\bm{A}}_{i}$. The above is clearly a linear function in the vector
${\bm{r}}_{i}\in\\{0,1\\}^{n}$ where
$[{\bm{r}}_{i}]_{\ell}=\mathbb{1}[\exists\;\ell-\text{clique}\in{\bm{A}}_{i}]$.
Furthermore, the graph is chosen so that ${\bm{r}}_{i}={\bm{x}}_{i}$. Thus,
any choice of weights ${\bm{a}},{\bm{b}},c$ corresponds to a halfspace for
${\bm{x}}_{i}$.

It only remains to be shown that any such halfspace can be written in the form
of Eq. 25. Note that

$\begin{split}f_{{\bm{a}},{\bm{b}},c}({\bm{x}}_{i},{\bm{A}}_{i})&=c^{\prime}+\begin{bmatrix}b_{1}&b_{2}&\cdots&b_{n}\end{bmatrix}{\bm{M}}{\bm{r}}_{i},\end{split}$
(26)

where $c^{\prime}$ is an updated constant and

${\bm{M}}=\begin{bmatrix}\operatorname{relu}(-1-a_{1})&2\operatorname{relu}(-1-2a_{1})&\cdots&n\operatorname{relu}(-1-na_{1})\\\
\operatorname{relu}(-1-a_{2})&2\operatorname{relu}(-1-2a_{2})&\cdots&n\operatorname{relu}(-1-na_{2})\\\
\vdots&\vdots&\ddots&\vdots\\\
\operatorname{relu}(-1-a_{n})&2\operatorname{relu}(-1-2a_{n})&\cdots&n\operatorname{relu}(-1-na_{n})\\\
\end{bmatrix}.$ (27)

Setting $a_{k}=(1/2-k)^{-1}$ makes the above matrix have non-zero entries on
the diagonal and every entry above the diagonal. Such a matrix is invertible
by Gaussian elimination so every possible halfspace can be constructed by
proper choice of the weights $b_{1},\dots,b_{n},c^{\prime}$. ∎

##  Appendix E GNN hardness proofs

###  E.1 Lower bound for Erdős–Rényi distributed graphs

Here, we will consider a model over graphs of $n$ nodes where adjacency
matrices ${\bm{A}}\sim\operatorname{Unif}(\\{0,1\\}^{n\times n})$ are drawn
uniformly from all possible directed graphs. This distribution can be viewed
as an Erdős–Rényi $G(n,p)$ over directed graphs including self edges with
$p=0.5$. Node features are fixed in this model to ${\bm{x}}\in\mathbb{R}^{n}$
where we set ${\bm{x}}=\bm{1}$ in all instances (i.e., node features are
trivial). Note, that we consider Boolean functions here w.l.o.g. as
$\\{0,1\\}$ valued as opposed to $\\{-1,+1\\}$ valued as in previous sections
since using $\\{0,1\\}$ valued functions in this section makes the story
easier to follow.

Our aim is to prove Theorem 3 copied below.

See 3

The function class $\mathcal{H}_{ER,n}$ that achieves the desired lower bound
are Boolean valued functions that depend on the counts of the degrees or
number of neighbors of the nodes. To describe this class, we first introduce
some notation. For a given graph ${\bm{A}}\in\\{0,1\\}^{n\times n}$, let
${\bm{d}}_{{\bm{A}}}\in[n]^{n}$ be a vector with entry
$[{\bm{d}}_{{\bm{A}}}]_{i}$ equal to the number of outgoing edges (or degree)
for node $i$. From here, let ${\bm{c}}_{{\bm{A}}}\in[n]^{n+1}$ denote the
count of the degrees where entry $[{\bm{c}}_{{\bm{A}}}]_{i}$ counts the number
of nodes that have $i-1$ outgoing edges in ${\bm{A}}$ or more formally

$[{\bm{c}}_{{\bm{A}}}]_{i}=\sum_{k=1}^{n}\mathbb{1}\left[[{\bm{d}}_{{\bm{A}}}]_{k}=i-1\right].$
(28)

Note that ${\bm{c}}_{{\bm{A}}}$ is a length $n+1$ vector since the number of
outgoing edges can be anything from $\\{0,1,\dots,n\\}$. Also note that this
vector is an invariant vector as permuting the nodes leaves the counts
invariant. Given a subset $S\subseteq[n+1]$ and a bit $b\in\\{0,1\\}$, the
functions we aim to learn take the form

$g_{S,b}({\bm{A}})=b+\sum_{i\in S}[{\bm{c}}_{{\bm{A}}}]_{i}\mod 2,$ (29)

which are a sort of parity function supported over the elements of the degree
counts ${\bm{c}}_{{\bm{A}}}$ in $\mathcal{S}$. Therefore, our hypothesis class
$\mathcal{H}$ consists of $2^{n+2}$ functions:

$\mathcal{H}_{ER,n}=\\{g_{S,b}:S\subseteq[n+1],b\in\\{0,1\\}\\}.$ (30)

First, we will show below that the class $\mathcal{H}_{ER,n}$ is exponentially
hard to learn in the SQ setting. Later in Lemma 15, we show that the functions
in $\mathcal{H}$ can be constructed by the GNN architecture thus completing
our proof.

###### Proof of Theorem 3.

To simplify the function $g_{S,b}$ even further, let
$\widehat{{\bm{c}}}_{{\bm{A}}}\in\\{0,1\\}^{n+1}$ be a boolean valued version
of ${\bm{c}}_{{\bm{A}}}$ where entry
$[\widehat{{\bm{c}}}_{{\bm{A}}}]_{i}=[{\bm{c}}_{{\bm{A}}}]_{i}\mod 2$. Then
the function $g_{S,b}$ is equivalently a parity function in
$\widehat{{\bm{c}}}_{{\bm{A}}}$:

$g_{S,b}({\bm{A}})=b+\sum_{i\in S}[{\bm{c}}_{{\bm{A}}}]_{i}\mod
2=(-1)^{b}\prod_{i\in S}(-1)^{[\widehat{{\bm{c}}}_{{\bm{A}}}]_{i}}.$ (31)

We will prove an SQ lower bound using Theorem 13 based on the
$(1-\eta)$-pairwise independence property of the function class
$\mathcal{H}_{ER,n}$ in Definition 12. Namely, note that for two graphs
${\bm{A}},{\bm{A}}^{\prime}\in\\{0,1\\}^{n\times n}$, if
$\widehat{{\bm{c}}}_{{\bm{A}}}\neq\widehat{{\bm{c}}}_{{\bm{A}}^{\prime}}$,
then the distribution of $(g_{S,b}({\bm{A}}),g_{S,b}({\bm{A}}^{\prime}))$ for
$g_{S,b}$ drawn uniformly from $\mathcal{H}_{ER,n}$ is equal to the product
distribution $\operatorname{Unif}(\\{0,1\\}^{2})$. This follows from the fact
that $g_{S,b}({\bm{A}})=(-1)^{b}\prod_{i\in
S}(-1)^{[\widehat{{\bm{c}}}_{{\bm{A}}}]_{i}}$ is a parity function of the
entries in $\widehat{{\bm{c}}}_{{\bm{A}}}$. Since $\mathcal{H}_{ER,n}$ spans
uniformly over all possible parity functions in $\\{0,1\\}^{n+1}$, the
pairwise independence properties follows from this whenever
$\widehat{{\bm{c}}}_{{\bm{A}}}\neq\widehat{{\bm{c}}}_{{\bm{A}}^{\prime}}$.

Thus, the parameter $1-\eta$ in the pairwise independence is equal to the
probability that
$\widehat{{\bm{c}}}_{{\bm{A}}}\neq\widehat{{\bm{c}}}_{{\bm{A}}^{\prime}}$. We
will establish that this probability is at least $1-O(\exp(-n^{\Omega(1)}))$
at the end of this proof thus guaranteeing
$\eta^{-1}=\Omega\left(\exp(n^{\Omega(1)})\right)$. Once this is established,
we can apply Theorem 13 and note that for a given function
$f^{*}\in\mathcal{H}_{ER,n}$, to distinguish $\mathcal{D}_{f^{*}}$ from
$\mathcal{D}_{\operatorname{Unif}(\mathcal{H}_{\rho})}$ as in Theorem 13, it
suffices to return a function $f$ that has classification error at most $1/4$
with respect to $f^{*}$. Since distinguishing the two distributions requires
$\frac{\tau^{2}}{2\eta}$ queries, we arrive at the final result by noting that
$\eta^{-1}=\Omega\left(\exp(n^{\Omega(1)})\right)$.

To finish this proof, we must establish that

$\mathbb{P}_{{\bm{A}}\sim\operatorname{Unif}(\\{0,1\\}^{n\times
n})}[\widehat{{\bm{c}}}_{{\bm{A}}}\neq\widehat{{\bm{c}}}_{{\bm{A}}^{\prime}}]\geq
1-O(\exp(-n^{\Omega(1)})).$ (32)

We will show this result by proving that for some constant $C>0$ and $p<1/2$,
for large enough $n$, there are at least $Cn^{p}$ entries of
$\widehat{{\bm{c}}}_{{\bm{A}}^{\prime}}$ such that the probability of those
entry being even or odd is at most $1/2+o(1)$. To continue, let us denote this
region

$\Delta_{p}=\\{0,1,\dots,n\\}\cap[n/2-Cn^{p},n/2+Cn^{p}],\quad p<1/2,$ (33)

as the set of integers within the range $[n/2-Cn^{p},n/2+Cn^{p}]$. This range
is chosen to coincide with the largest entries of ${\bm{c}}_{{\bm{A}}}$.

The degree of any node is independently distributed as binomial with
distribution $\operatorname{bin}(n,0.5)$. To asymptotically expand $i$ around
the peak of the binomial coefficient, let us index $i=n/2+r$ for given $r$
such that $i\in\Delta_{p}$. Therefore, for any such $n/2+r=i\in\Delta_{p}$, we
have that any given node $j\in[n]$ has degree ${\bm{d}}_{k}=i$ with
probability

$\mathbb{P}[{\bm{d}}_{k}=i]=2^{-n}{\binom{n}{n/2+r}}=\frac{2}{\sqrt{\pi
n}}(1+O(1/n)).$ (34)

In the above, we use the asymptotic approximation [Spencer, 2014] for the
binomial coefficient around its peak $\binom{n}{n/2}$ of

${\binom{n}{\frac{n+r}{2}}}=(1+O(1/n))\sqrt{\frac{2}{\pi
n}}2^{n}\exp\left(-\frac{r^{2}}{2n}+O(r^{3}/n^{2})\right).$ (35)

Thus, the marginal distribution of any entry $i\in\Delta_{p}$ of
$[{\bm{c}}_{{\bm{A}}}]_{i}\sim\operatorname{bin}\left(n,\frac{2}{\sqrt{\pi
n}}(1+O(1/n))\right)$ is also binomial distributed with probability converging
to $\frac{2}{\sqrt{\pi n}}$. For entry $i\in\Delta_{p}$, this implies that
$\mathbb{E}[{\bm{c}}_{{\bm{A}}}]_{i}=O(\sqrt{n})$. From the Chernoff
concentration bound for the binomial distribution, we also have that for any
$i\in\Delta_{p}$ and for any $\delta>0$:

$\mathbb{P}\left[\left|[{\bm{c}}_{{\bm{A}}}]_{i}-\mathbb{E}[{\bm{c}}_{{\bm{A}}}]_{i}\right|\geq
n^{1/4+\delta}\right]\leq
2\exp\left(-\frac{n^{1/2+2\delta}(\mathbb{E}[{\bm{c}}_{{\bm{A}}}]_{i})^{-1}}{2}\right)=2\exp\left(-\Omega(n^{2\delta})\right).$
(36)

By a union bound over all $i\in\Delta_{p}$, we therefore have that with
exponentially high probability,
$\left|[{\bm{c}}_{{\bm{A}}}]_{i}-\mathbb{E}[{\bm{c}}_{{\bm{A}}}]_{i}\right|\leq
n^{1/4+\delta}$ for all $i\in\Delta_{p}$.

Summarizing the previous results, we have that there are $2Cn^{p}$ entries of
$[{\bm{c}}_{{\bm{A}}}]$ where for any entry $i\in\Delta_{p}$, it holds that
for some arbitrarily small $\delta>0$

$\frac{2}{\sqrt{\pi}}\sqrt{n}(1+O(1/n))-n^{1/4+\delta}\leq[{\bm{c}}_{{\bm{A}}}]_{i}\leq\frac{2}{\sqrt{\pi}}\sqrt{n}(1+O(1/n))+n^{1/4+\delta}.$
(37)

Since the probability that $[{\bm{c}}_{{\bm{A}}}]_{i}$ is equal to any number
$m$ is at most $O(n^{-1/4})$, the probability that $[{\bm{c}}_{{\bm{A}}}]_{i}$
is odd is at most $1/2+O(n^{-1/4})$. We now aim to show that conditioned on
other entries of ${\bm{c}}_{{\bm{A}}}$, the probability that a given entry of
${\bm{c}}_{{\bm{A}}}$ is odd or even is roughly independent of this
conditioning.

Given a set of known values $S_{known}\subset\Delta_{p}$ equal to
$[{\bm{c}}_{{\bm{A}}}]_{k}=m_{k}$ for $k\in S_{known}$, the marginal
distribution of any entry of $[{\bm{c}}_{{\bm{A}}}]_{j}$ for $j\notin
S_{known}$ conditioned on $[{\bm{c}}_{{\bm{A}}}]_{k}$ for $k\in S_{known}$ is
also binomial. We now estimate the parameters of this binomial distribution.
Given the convergence guarantee in Eq. 37, $\sum_{k\in
S_{known}}[{\bm{c}}_{{\bm{A}}}]_{k}=O(|S_{known}|n^{1/2})=O(n^{1/2+p})$ which
is $o(n)$ since $p<1/2$. Therefore for any set $S_{known}\subseteq\Delta_{p}$,
there are at least $n-o(n)$ remaining nodes whose degree or number of outgoing
edges are left to distribute. These remaining nodes cannot have degrees whose
values fall in $S_{known}$, but note that since $|S_{known}|\leq O(n^{p})$ and
$p<1/2$, there are still $\Omega(\sqrt{n})$ buckets (values of the possible
degrees) remaining each with probability proportional to $\Theta(1/\sqrt{n})$
by Eq. 35. Therefore, from the convergence guarantee of Eq. 37, we have

$([{\bm{c}}_{{\bm{A}}}]_{j}\;|\;[{\bm{c}}_{{\bm{A}}}]_{k}=m_{k},k\in
S_{known})\sim\operatorname{bin}\left(n-o(n),\Theta(1/\sqrt{n})\right).$ (38)

Thus, the conditional probability that $[{\bm{c}}_{{\bm{A}}}]_{j}$ is equal to
any number $m$ is at most $O(n^{-1/4})$, and this implies that the probability
that $[{\bm{c}}_{{\bm{A}}}]_{j}$ is odd is at most $1/2+O(n^{-1/4})$ even when
conditioned on entries in $S_{known}$.

Finally, let us calculate the desired probability that
$\widehat{{\bm{c}}}_{{\bm{A}}}=\widehat{{\bm{c}}}_{{\bm{A}}^{\prime}}$. For
ease of notation, let us index the elements of $\Delta_{p}$ as
$r_{1},\dots,r_{N}$. We have for large enough $n$ that

$\begin{split}\mathbb{P}[\widehat{{\bm{c}}}_{{\bm{A}}}=\widehat{{\bm{c}}}_{{\bm{A}}^{\prime}}]&\leq\mathbb{P}\left[[\widehat{{\bm{c}}}_{{\bm{A}}}]_{r_{1}}=[\widehat{{\bm{c}}}_{{\bm{A}}^{\prime}}]_{r_{1}},\dots,[\widehat{{\bm{c}}}_{{\bm{A}}}]_{r_{N}}=[\widehat{{\bm{c}}}_{{\bm{A}}^{\prime}}]_{r_{N}}\right]\\\
&=\mathbb{P}\left[[\widehat{{\bm{c}}}_{{\bm{A}}}]_{r_{1}}=[\widehat{{\bm{c}}}_{{\bm{A}}^{\prime}}]_{r_{1}}\right]\;\mathbb{P}\bigl{[}[\widehat{{\bm{c}}}_{{\bm{A}}}]_{r_{2}}=[\widehat{{\bm{c}}}_{{\bm{A}}^{\prime}}]_{r_{2}}\;|\;[\widehat{{\bm{c}}}_{{\bm{A}}}]_{r_{1}}=[\widehat{{\bm{c}}}_{{\bm{A}}^{\prime}}]_{r_{1}}\bigr{]}\times\\\
&\quad\times\cdots\mathbb{P}\bigl{[}[\widehat{{\bm{c}}}_{{\bm{A}}}]_{r_{N}}=[\widehat{{\bm{c}}}_{{\bm{A}}^{\prime}}]_{r_{N}}\;|\;[\widehat{{\bm{c}}}_{{\bm{A}}}]_{r_{1}}=[\widehat{{\bm{c}}}_{{\bm{A}}^{\prime}}]_{r_{1}},\dots,[\widehat{{\bm{c}}}_{{\bm{A}}}]_{r_{N-1}}=[\widehat{{\bm{c}}}_{{\bm{A}}^{\prime}}]_{r_{N-1}}\bigr{]}\\\
&\leq\left(\frac{1}{2}+O(n^{-1/4})\right)^{|\Delta_{p}|}\\\
&=O(\exp(-\Omega(|\Delta_{p}|)))\\\ &=O(\exp(-\Omega(n^{p}))).\end{split}$
(39)

Since $\mathbb{P}_{{\bm{A}}\sim\operatorname{Unif}(\\{0,1\\}^{n\times
n})}[\widehat{{\bm{c}}}_{{\bm{A}}}\neq\widehat{{\bm{c}}}_{{\bm{A}}^{\prime}}]=1-\mathbb{P}_{{\bm{A}}\sim\operatorname{Unif}(\\{0,1\\}^{n\times
n})}[\widehat{{\bm{c}}}_{{\bm{A}}}=\widehat{{\bm{c}}}_{{\bm{A}}^{\prime}}]$,
this completes the proof. ∎

Finally, we need to show that the functions $g_{S,b}$ can be constructed as
ReLU networks. Before proceeding to give this construction, we should note
that it is known that any Boolean function that can be computed in time
$O(T(n))$ can also be expressed by a neural network of size $O(T(n^{2}))$
[Parberry, 1994, Shamir, 2018]. Our construction will be specific to the GNN
class we consider and show that the boolean functions in $g_{S,b}$ are
similarly efficiently constructible.

######  Lemma 15 ($g_{S,b}$ as GNN).

For a GNN $f(\cdot)$ of the form of Eq. 5 with hidden layer widths equal to
$k_{1}=O(n)$ and $k_{2}=O(n)$, there exist weights
${\bm{a}},{\bm{b}}\in\mathbb{R}^{k_{1}},{\bm{u}},{\bm{v}}\in\mathbb{R}^{k_{2}},{\bm{W}}\in\mathbb{R}^{k_{2}\times
k_{1}}$ such that for any $S\subseteq[n+1]$, $b\in\\{0,1\\}$,
$g_{S,b}({\bm{A}})=f({\bm{A}};{\bm{a}},{\bm{b}},{\bm{u}},{\bm{v}},{\bm{W}})$
for all inputs ${\bm{A}}\in\\{0,1\\}^{n\times n}$.

###### Proof.

Let us fix a given $S\subseteq[n+1]$, $b\in\\{0,1\\}$. The function we want to
represent is

$g_{S,b}({\bm{A}})=b+\sum_{i\in S}[{\bm{c}}_{{\bm{A}}}]_{i}\mod 2,$ (40)

where as a reminder, $[{\bm{c}}_{{\bm{A}}}]_{i}$ counts the number of nodes
that have $i-1$ outgoing edges (see Eq. 28).

For the construction, we will set $k_{1}=k_{2}=n+1$. As a reminder, the two
layers of the GNN take the form:

$\begin{split}[f_{{\bm{a}},{\bm{b}}}^{(1)}({\bm{A}})]_{i}&=\bm{1}_{n}^{\top}\sigma\left(a_{i}+b_{i}{\bm{A}}{\bm{x}}\right)\quad\text{(output
of channel $i\in[k_{1}]$)}\\\
f_{{\bm{u}},{\bm{v}},{\bm{W}}}^{(2)}({\bm{h}})&=\sum_{i=1}^{k_{2}}u_{i}\sigma(\langle{\bm{W}}_{:,i},{\bm{h}}\rangle+v_{i}).\end{split}$
(41)

Throughout this construction, the node features ${\bm{x}}=\bm{1}$ are
constant.

For the first layer, we will choose $a_{i}=2-i$ and $b_{i}=1$ for all $i$.
This first layer can only be a function of
${\bm{d}}_{{\bm{A}}}={\bm{A}}\bm{1}$ where entry $i$ is a count of the number
of outgoing edges. Since entries in ${\bm{d}}_{{\bm{A}}}$ take values in
$\\{0,1,\dots,n\\}$, we can equivalently write this first layer output as a
function of the counts ${\bm{c}}_{{\bm{A}}}$. In fact, this first layer is
linear as a function of ${\bm{c}}_{{\bm{A}}}$ and equal to

${\bm{h}}={\bm{M}}{\bm{c}}_{{\bm{A}}}=\begin{bmatrix}1&2&\cdots&n&n+1\\\
0&1&\cdots&n-1&n\\\ \vdots&\vdots&\ddots&\vdots&\vdots\\\ 0&0&\cdots&1&2\\\
0&0&\cdots&0&1\\\ \end{bmatrix}{\bm{c}}_{{\bm{A}}},$ (42)

where the matrix ${\bm{M}}$ is clearly full rank and invertible since it
strictly positive along the diagonal and in entries above the diagonal.

Figure 2: Sample form of function $h(x)$ used in constructing
$g_{S,b}({\bm{A}})$ as a GNN. For the construction, we will have $x=\sum_{i\in
S}[{\bm{c}}_{{\bm{A}}}]_{i}$.

Our aim in the second layer is to take the sum $b+\sum_{i\in
S}[{\bm{c}}_{{\bm{A}}}]_{i}$ and take its parity (i.e. $\mod 2$). To achieve
this via sums of ReLUs $\sigma$ we can consider the function $h(\cdot)$ of the
form below:

$h(x)=\sigma(x+b)-2\sigma(x+b-1)+2\sigma(x+b-2)-\cdots+(-1)^{n}\sigma(x+b-n).$
(43)

A sample version of this plot is in Fig. 2. Setting $x=b+\sum_{i\in
S}[{\bm{c}}_{{\bm{A}}}]_{i}$ achieves the end goal with $O(n)$ hidden nodes.

To provide actual values for the second layer parameters, let
$\bm{1}_{S}\in{0,1}^{n+1}$ denote the vector where entry $i$ is equal to $1$
if $i\in S$ and equal to $0$ otherwise. Then, in the second layer, we have

  * •

${\bm{W}}_{:,i}=({\bm{M}}^{-1})^{\top}\bm{1}_{S}$ for all columns $i\in[n+1]$

  * •

$a_{1}=1$ and $a_{i}=2(-1)^{i-1}$ for all $i\in\\{2,\dots,n+1\\}$

  * •

$v_{i}=b-i+1$ for all $i\in[n+1]$

Using the values above, we have

$\begin{split}f_{{\bm{u}},{\bm{v}},{\bm{W}}}^{(2)}({\bm{h}})&=\sigma(\langle({\bm{M}}^{-1})^{\top}\bm{1}_{S},{\bm{M}}{\bm{c}}_{{\bm{A}}}\rangle+b)+\sum_{i=2}^{k_{2}}(-1)^{i-1}2\sigma(\langle({\bm{M}}^{-1})^{\top}\bm{1}_{S},{\bm{M}}{\bm{c}}_{{\bm{A}}}\rangle+b-i+1)\\\
&=\sigma\left(b+\sum_{i\in
S}[{\bm{c}}_{{\bm{A}}}]_{i}\right)+\sum_{i=2}^{k_{2}}(-1)^{i-1}2\sigma\left(\sum_{i\in
S}[{\bm{c}}_{{\bm{A}}}]_{i}+b-i+1\right)\\\
&=h\left([{\bm{c}}_{{\bm{A}}}]_{i}\right).\end{split}$ (44)

As $\sum_{i\in S}[{\bm{c}}_{{\bm{A}}}]_{i}$ only takes values in
$\\{0,1,\dots,n\\}$, the above is supported on all of its inputs and this
completes the proof. ∎

###  E.2 Invariant Hermite polynomial

#### Multivariate Hermite polynomials

To handle the graph setting with scaling feature dimension, we first extend
the technical machinery related to Hermite polynomials to the invariant
settings we consider.

Recall from Diakonikolas et al. [2020] that the $d$-variate normalized Hermite
polynomial is defined for every multi-index
$J=(J_{1},\ldots,J_{d})\in\mathbb{N}^{d}$ as:

$H_{J}:\mathbb{R}^{d}\to\mathbb{R}:(x_{1},x_{2},\ldots,x_{d})\mapsto\prod_{j=1}^{d}H_{J_{j}}(x_{j}),\qquad
H_{i}:\mathbb{R}\to\mathbb{R}:x\mapsto\frac{H_{e_{i}}(x)}{\sqrt{i!}},i\in\mathbb{N},$
(45)

where $H_{e_{i}}$ is the usual probabilist Hermite polynomial.

#### Graph-invariant Hermite polynomial

Fix a graph ${\mathcal{G}}$ on $n$ vertices and let
${\bm{A}}:={\bm{A}}({\mathcal{G}})\in\mathbb{R}^{n\times n}$ be a graph shift
operator (say the Laplacian). Define, for any multi-index
$J\in\mathbb{N}^{d}$:

$H^{{\bm{A}}}_{J}:\mathbb{R}^{n\times
d}\to\mathbb{R}:{\bm{X}}\mapsto\frac{1}{\sqrt{n}}\sum_{v=1}^{n}H_{J}(\left({\bm{A}}{\bm{X}}\right)_{v}).$
(46)

We have the following orthogonality property:

######  Lemma 16 (Orthogonality of invariant Hermite polynomial).

For a fixed ${\bm{A}}\in\mathbb{R}^{n\times n}$ such that
$\sum_{u=1}^{n}{\bm{A}}_{vu}^{2}=1$ for all $v\in[n]$, $\langle
H^{{\bm{A}}}_{J},H^{{\bm{A}}}_{J^{\prime}}\rangle_{L^{2}(\mathcal{N})}=\delta_{J=J^{\prime}}\cdot
c_{J}$ for all $J,J^{\prime}\in\mathbb{N}^{d}\backslash\\{0\\}$ for some
$c_{J}\in\Theta(1)$ that depends only on ${\bm{A}}$.

###### Proof.

Let $\Sigma=\Sigma(\rho)=\begin{bmatrix}1&\rho\\\ \rho&1\end{bmatrix}$ for
some $\rho\in(-1,1)$. We will also use Mehler’s formula [Kibble, 1945] which
states:

$\frac{p_{\mathcal{N}(0,\Sigma)}(x,y)}{p_{\mathcal{N}}(x)p_{\mathcal{N}}(y)}=\sum_{m=0}^{\infty}\rho^{m}H_{m}(x)H_{m}(y),$
(47)

where $p_{\mathcal{N}(0,\Sigma)}$ is the pdf of the bivariate normal
distribution with mean $0$ and covariance $\Sigma$ and $p_{\mathcal{N}}$ is
the pdf of the univariate standard normal distribution.

We first show orthogonality of Hermite polynomials evaluated on dependent
random variables. For all $i,j\in\mathbb{N}$:

$\displaystyle\mathbb{E}_{(x,y)\sim\mathcal{N}_{2}(0,\Sigma)}[H_{i}(x)H_{j}(y)]$
$\displaystyle=\int_{\mathbb{R}}\int_{\mathbb{R}}H_{i}(x)H_{j}(y)p_{\mathcal{N}(0,\Sigma)}(x,y)\textup{d}x\textup{d}y$
(48)
$\displaystyle=\int_{\mathbb{R}}\int_{\mathbb{R}}H_{i}(x)H_{j}(y)p_{\mathcal{N}}(x)p_{\mathcal{N}}(y)\cdot\frac{p_{\mathcal{N}(0,\Sigma)}(x,y)}{p_{\mathcal{N}}(x)p_{\mathcal{N}}(y)}\textup{d}x\textup{d}y$
(49)
$\displaystyle=\int_{\mathbb{R}}\int_{\mathbb{R}}H_{i}(x)H_{j}(y)\sum_{m=0}^{\infty}\rho^{m}H_{m}(x)H_{m}(y)\textup{d}P(x)\textup{d}P(y)$
(50)
$\displaystyle=\sum_{m=0}^{\infty}\rho^{m}\left(\int_{\mathbb{R}}H_{i}(x)H_{m}(x)\textup{d}P(x)\right)\left(\int_{\mathbb{R}}H_{j}(y)H_{m}(y)\textup{d}P(y)\right)$
(51) $\displaystyle=\sum_{m=0}^{\infty}\rho^{m}\delta_{i=m}\delta_{j=m}$ (52)
$\displaystyle=\rho^{i}\delta_{i=j}.$ (53)

We are now ready to present the proof of Lemma 16. Pick
$J,J^{\prime}\in\mathbb{N}^{d}\backslash\\{0\\}$. We have:

$\displaystyle\langle
H_{J}^{{\bm{A}}},H_{J^{\prime}}^{{\bm{A}}}\rangle_{L^{2}(\mathcal{N})}$
$\displaystyle=\frac{1}{n}\sum_{v,v^{\prime}\in[n]}\mathbb{E}_{{\bm{X}}\sim\mathcal{N}^{n\times
d}}[H_{J}(({\bm{A}}{\bm{X}})_{v})H_{J^{\prime}}(({\bm{A}}{\bm{X}})_{v^{\prime}})]$
(54)
$\displaystyle=\frac{1}{n}\sum_{v,v^{\prime}\in[n]}\mathbb{E}_{{\bm{Z}}\sim
P}[H_{J}({\bm{z}}_{v})H_{J^{\prime}}({\bm{z}}_{v^{\prime}})],$ (55)

where ${\bm{z}}_{v}$ is the $v$-th row of ${\bm{Z}}\in\mathbb{R}^{n\times d}$
and $P$ is a multivariate normal distribution over $\mathbb{R}^{n\times d}$
with mean $0$ and covariance:

$\displaystyle\mathrm{Cov}({\bm{Z}}_{vj},{\bm{Z}}_{v^{\prime}j^{\prime}})=\begin{cases}0,&\text{
if }j\neq j^{\prime}\\\
\sum_{u=1}^{n}{\bm{A}}_{vu}{\bm{A}}_{v^{\prime}u}=:\rho_{vv^{\prime}},&\text{
otherwise.}\end{cases}$ (56)

Furthermore, for any $v$ and any $j\neq j^{\prime}$, ${\bm{Z}}_{vj}$ is
independent of ${\bm{Z}}_{v^{\prime}j^{\prime}}$. This is because columns of
${\bm{X}}$ are independent from one another. Note that when $v=v^{\prime}$ and
$j=j^{\prime}$ the variance is $\sum_{u=1}^{n}{\bm{A}}_{vu}^{2}=1$ by our
assumption on ${\bm{A}}$.

By the dependence structure of $P$, one can write:

$\displaystyle\mathbb{E}_{{\bm{Z}}\sim
P}[H_{J}({\bm{z}}_{v})H_{J^{\prime}}({\bm{z}}_{v^{\prime}})]$
$\displaystyle=\prod_{j=1}^{d}\mathbb{E}_{{\bm{Z}}\sim
P}[H_{J_{j}}({\bm{Z}}_{vj})H_{J^{\prime}_{j}}({\bm{Z}}_{v^{\prime}j})]$ (57)
$\displaystyle=\prod_{j=1}^{d}\mathbb{E}_{({\bm{Z}}_{vj},{\bm{Z}}_{v^{\prime}j})\sim\mathcal{N}(0,\Sigma(\rho_{vv^{\prime}}))}[H_{J_{j}}({\bm{Z}}_{vj})H_{J^{\prime}_{j}}({\bm{Z}}_{v^{\prime}j})]$
(58)
$\displaystyle=\prod_{j=1}^{d}\rho_{vv^{\prime}}^{J_{i}}\delta_{J_{i}=J^{\prime}_{i}}$
(59) $\displaystyle=\rho_{vv^{\prime}}^{|J|}\delta_{J=J^{\prime}}.$ (60)

Therefore,

$\displaystyle\langle
H_{J}^{{\bm{A}}},H_{J^{\prime}}^{{\bm{A}}}\rangle_{L^{2}(\mathcal{N})}$
$\displaystyle=\delta_{J=J^{\prime}}\cdot\frac{1}{n}\sum_{v,v^{\prime}\in[n]}\rho_{vv^{\prime}}^{|J|}=\delta_{J=J^{\prime}}\cdot
c_{J},$ (61)

where $c_{J}=\frac{1}{n}\sum_{v,v^{\prime}\in[n]}\rho_{vv^{\prime}}^{|J|}$. ∎

######  Remark 17.

The quantity $\rho_{vv^{\prime}}$ deserves an explanation. Essentially,
$\rho_{vv^{\prime}}$ measures how similar two nodes $v$ and $v^{\prime}$ are,
in term of composition of its incoming edges. Algebraically, it is the
$(vv^{\prime})$-th entry of the gram matrix ${\bm{A}}{\bm{A}}^{\top}$. We can
get the correct normalization of ${\bm{A}}$ in the premise of Lemma 16 by
starting with a graph ${\mathcal{G}}$, adding a self-loop to each node of
${\mathcal{G}}$ and set
${\bm{A}}=({\bm{D}}^{{\dagger}})^{\frac{1}{2}}{\bm{M}}$ where ${\bm{M}}$ is
the adjacency matrix of ${\mathcal{G}}$ (with self-loops) and ${\bm{D}}$ is
the in-degree matrix, we can ensure that $\rho_{vv^{\prime}}\in[0,1]$ for each
$v\neq v^{\prime}$:

$\rho_{vv^{\prime}}=\sum_{u=1}^{n}{\bm{A}}_{vu}{\bm{A}}_{v^{\prime}u}=\frac{|N_{v}\cap
N_{v^{\prime}}|}{\sqrt{\deg(v)\deg(v^{\prime})}}\leq 1,$ (62)

where $N_{v}$ is the set of vertices with an incoming edge to $v\in[n]$.

This immediately gives a bound $\max_{J}c_{J}\leq n^{2}$. Finally, it is also
worth noting that $c_{J}$ does not vanish under the above normalization, since
adding self-loops before the normalization to ${\mathcal{G}}$ ensures that
$\rho_{vv}=1$ for all $v\in[n]$ even if $v$ was an isolated vertex in the
original graph. Therefore, without loss of generality, we can assume that
$c_{J}\in\Theta(1)$.

It is possible to use a different normalization, such as the Laplacian, as
long as one can ensure that $\rho_{vv^{\prime}},c_{J}=\Theta(1)$. The
requirement that rows of ${\bm{A}}$ sums to $1$ is also not strictly necessary
but does greatly simplify our exposition.

Note that we are not claiming a complete orthonormal basis for the space
$L^{2}(\mathbb{R}^{n\times d},\mathcal{N}^{n\times d})$. However, all of the
functions we care about will lie in
$\mathcal{H}^{(d)}:=\textup{span}\left(H_{J}^{{\bm{A}}}\right)_{J\in\mathbb{N}^{d}}$.
For a fixed $n$ and a function $f\in\mathcal{H}^{(d)}$, we define
$\hat{f}_{J}=\langle f,H_{J}^{{\bm{A}}}\rangle_{L^{2}(\mathcal{N})}$ and the
degree $m$ part of $f$:
$f^{[m]}({\bm{X}})=\sum_{|J|=m}\hat{f}_{J}H_{J}^{{\bm{A}}}({\bm{X}})$.

Facts about gradients of Hermite polynomial also carry over to invariant
Hermite polynomials:

######  Lemma 18.

For a fixed graph ${\mathcal{G}}$ on $n$ vertices and corresponding
${\bm{A}}$, we have:

  1. 1.

$\frac{\partial}{\partial{\bm{X}}_{vj}}H_{M}^{{\bm{A}}}(X)=\frac{1}{\sqrt{n}}\sum_{v^{\prime}=1}^{n}\sqrt{M_{j}}H_{M-E_{j}}(({\bm{A}}{\bm{X}})_{v^{\prime}}){\bm{A}}_{v^{\prime}v}$,
for all $v\in[n],j\in[d],M\in\mathbb{N}^{d}$, where $E_{j}$ is the vector with
$1$ at the $j^{\prime}$-th position and $0$ everywhere else.

  2. 2.

$\mathbb{E}_{{\bm{X}}\sim\mathcal{N}^{n\times d}}[\langle\nabla
H_{M}^{{\bm{A}}}({\bm{X}}),\nabla
H_{L}^{{\bm{A}}}({\bm{X}})\rangle]=\delta_{M=L}\cdot|M|c_{M}$.

  3. 3.

$\mathbb{E}_{{\bm{X}}\sim\mathcal{N}^{n\times
d}}[\langle\nabla^{l}H_{M}^{{\bm{A}}}({\bm{X}}),\nabla^{l}H_{L}^{{\bm{A}}}({\bm{X}})\rangle]=\delta_{M=L}\cdot|M|(|M|-1)\ldots(|M|-l+1)c_{M}$.

In the above,
$c_{M}=n^{-1}\sum_{v^{\prime},v^{\prime\prime}=1}^{n}\rho_{v^{\prime}v^{\prime\prime}}^{|M|}$.

###### Proof.

  1. 1.

We have:

$\displaystyle\frac{\partial}{\partial{\bm{X}}_{vj}}H_{M}^{{\bm{A}}}({\bm{X}})$
$\displaystyle=\frac{1}{\sqrt{n}}\sum_{v^{\prime}=1}^{n}(\nabla
H_{M}(({\bm{A}}{\bm{X}})_{v^{\prime}}))_{j}{\bm{A}}_{v^{\prime}v}$ (63)
$\displaystyle=\frac{1}{\sqrt{n}}\sum_{v^{\prime}=1}^{n}\sqrt{M_{j}}H_{M-E_{j}}(({\bm{A}}{\bm{X}})_{v^{\prime}}){\bm{A}}_{v^{\prime}v},$
(64)

where $E_{j^{\prime}}$ is the vector with $1$ at the $j^{\prime}$-th position
and $0$ everywhere else. The second equality is a known fact for the
multivariate Hermite polynomials (See Appendix C of Diakonikolas et al.
[2020]).

  2. 2.

Sum the above over all pairs $v,v^{\prime}$ and take expectation to get:

$\displaystyle\mathbb{E}_{{\bm{X}}\sim\mathcal{N}^{n\times d}}[\langle\nabla
H_{M}^{{\bm{A}}}({\bm{X}}),\nabla H_{L}^{{\bm{A}}}({\bm{X}})\rangle]$ (65)
$\displaystyle=$
$\displaystyle\sum_{v=1}^{n}\sum_{j=1}^{d}\frac{1}{n}\mathbb{E}_{{\bm{X}}\sim\mathcal{N}^{n\times
d}}\left[\sum_{v^{\prime},v^{\prime\prime}=1}^{n}\left(\nabla
H_{M}(({\bm{A}}{\bm{X}})_{v^{\prime}})\right)_{j}{\bm{A}}_{v^{\prime}v}\left(\nabla
H_{L}(({\bm{A}}{\bm{X}})_{v^{\prime\prime}})\right)_{j}{\bm{A}}_{v^{\prime\prime}v}\right].$
(66)

By the orthogonality property of Hermite multivariate polynomial under
dependent multivariate Gaussian (equation 60), we can compute:

$\displaystyle\mathbb{E}_{{\bm{X}}\sim\mathcal{N}^{n\times
d}}[H_{M-E_{j}}(({\bm{A}}{\bm{X}})_{v^{\prime}})H_{L-E_{j}}(({\bm{A}}{\bm{X}})_{v^{\prime\prime}})]=\rho_{v^{\prime}v^{\prime\prime}}^{|M|-1}\delta_{M=L},$
(67)

and,

$\displaystyle\mathbb{E}_{{\bm{X}}\sim\mathcal{N}^{n\times d}}[\langle\nabla
H_{M}^{{\bm{A}}}({\bm{X}}),\nabla H_{L}^{{\bm{A}}}({\bm{X}})\rangle]$
$\displaystyle=\frac{1}{n}\sum_{v,v^{\prime},v^{\prime\prime}=1}^{n}{\bm{A}}_{v^{\prime}v}{\bm{A}}_{v^{\prime\prime}v}\rho_{v^{\prime}v^{\prime\prime}}^{|M|-1}\delta_{M=L}\sum_{j=1}^{d}M_{j}$
(68)
$\displaystyle=\delta_{M=L}\cdot\frac{|M|}{n}\sum_{v^{\prime},v^{\prime\prime}=1}^{n}\rho_{v^{\prime}v^{\prime\prime}}^{|M|-1}\sum_{v=1}^{n}{\bm{A}}_{v^{\prime}v}{\bm{A}}_{v^{\prime\prime}v}$
(69)
$\displaystyle=\delta_{M=L}\cdot\frac{|M|}{n}\sum_{v^{\prime},v^{\prime\prime}=1}^{n}\rho_{v^{\prime}v^{\prime\prime}}^{|M|}.$
(70)

  3. 3.

Recursively apply the first part to get:

$\displaystyle\frac{\partial^{l}}{\partial{\bm{X}}_{v_{1}j_{1}},\ldots,{\bm{X}}_{v_{l}j_{l}}}H_{M}^{G}({\bm{X}})$
(71) $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{n}}\sum_{v^{\prime}=1}^{n}\left(\prod_{p=1}^{l}{\bm{A}}_{v^{\prime}v_{p}}\sqrt{\left(M-\sum_{p^{\prime}=1}^{p}E_{j_{p^{\prime}}}\right)_{j_{p}}}\right)H_{M-\sum_{p=1}^{l}e_{j_{p}}}(({\bm{A}}{\bm{X}})_{v^{\prime}}).$
(72)

Summing indices and passing through the expectation to get:

$\displaystyle\mathbb{E}_{{\bm{X}}\sim\mathcal{N}^{n\times
d}}[\langle\nabla^{l}H_{M}^{{\bm{A}}}({\bm{X}}),\nabla^{l}H_{L}^{{\bm{A}}}({\bm{X}})\rangle]$
(73) $\displaystyle=$
$\displaystyle\frac{\delta_{M=L}}{n}\left(\sum_{j_{1},\ldots,j_{l}}\left(M-\sum_{p^{\prime}=1}^{p}E_{j_{p^{\prime}}}\right)_{j_{p}}\right)\sum_{v^{\prime},v^{\prime\prime}}\left(\sum_{v=1}^{n}{\bm{A}}_{v^{\prime}v}{\bm{A}}_{v^{\prime\prime}v}\right)^{l}\rho_{v^{\prime}v^{\prime\prime}}^{|M|-l}$
(74) $\displaystyle=$
$\displaystyle\delta_{M=L}|M|(|M|-1)\ldots(|M|-l)\cdot\frac{1}{n}\sum_{v^{\prime},v^{\prime\prime}=1}^{n}\rho_{v^{\prime}v^{\prime\prime}}^{|M|}.$
(75)

∎

######  Corollary 19.

Fix $k,d,e\in\mathbb{N}$. Let $p,q$ be linear combinations of
$H_{M}^{{\bm{A}}}:\mathbb{R}^{d}\to\mathbb{R}$ with $M\in\mathbb{N}^{d}$
varying among all multi-indices of size $k=|M|$. Then,

$\left\langle\nabla^{k}p({\bm{X}}),\nabla^{k}q({\bm{X}})\right\rangle=k!\cdot\mathbb{E}_{{\bm{X}}\sim\mathcal{N}^{n\times
d}}[p({\bm{X}})q({\bm{X}})].$ (76)

###### Proof.

Let $p({\bm{X}})=\sum_{M:|M|=k}b_{M}H_{M}^{{\bm{A}}}({\bm{X}})$ and
$q(X)=\sum_{K:|K|=k}f_{K}H_{K}^{{\bm{A}}}({\bm{X}})$. Then:

$\mathbb{E}_{{\bm{X}}\sim\mathcal{N}^{n\times
d}}[p({\bm{X}})q({\bm{X}})]=\sum_{M:|M|=k}f_{M}b_{M}c_{M},$ (77)

where as before,
$c_{M}=n^{-1}\sum_{v^{\prime},v^{\prime\prime}=1}^{n}\rho_{v^{\prime}v^{\prime\prime}}^{|M|}$.

From the previous lemma,

$\displaystyle\langle\nabla^{k}p({\bm{X}}),\nabla^{k}q({\bm{X}})\rangle$
$\displaystyle=\mathbb{E}_{{\bm{X}}\sim\mathcal{N}^{n\times
d}}\left[\left\langle\sum_{M:|M|=k}b_{M}\nabla^{k}H_{M}^{{\bm{A}}}({\bm{X}}),\sum_{K:|K|=k}f_{M}\nabla^{k}H_{K}^{{\bm{A}}}({\bm{X}})\right\rangle\right]$
(78) $\displaystyle=\sum_{M,K:|M|=|K|=k}b_{M}f_{M}k!\delta_{M=K}c_{M}$ (79)
$\displaystyle=k!\sum_{M:|M|=k}b_{M}f_{M}c_{M}$ (80)
$\displaystyle=k!\cdot\mathbb{E}_{{\bm{X}}\sim\mathcal{N}^{n\times
d}}[p({\bm{X}})q({\bm{X}})].$ (81)

∎

###  E.3 Proof of Theorem 4

The proof of Theorem 4 is derived from the construction in Diakonikolas et al.
[2020]. Hence, much of the setup will be taken from their construction.
However, care must be taken to use the proper invariant Hermite polynomial
basis we derived in the previous subsection. For completeness, we extend their
formalism here as needed and note where proofs are directly derived from their
work.

We first verify that the orthogonal system of invariant Hermite polynomial
subsumes $1$-hidden-layer GNNs. We formally show this below to continue the
proof later on in the basis of the invariant Hermite polynomials. As a
reminder, we defined
$\mathcal{H}^{(d)}:=\textup{span}\left(H_{J}^{{\bm{A}}}\right)_{J\in\mathbb{N}^{d}}$.

######  Lemma 20 (Inner product decomposition).

For every $k\in\mathbb{N}$ and ${\mathcal{G}}\in\mathbb{G}_{n}$,
$F^{d,k}_{{\mathcal{G}}}\subset\mathcal{H}^{(d)}$.

###### Proof.

Pick $f\in F^{d,k}_{{\mathcal{G}}}$, with a corresponding ${\bm{W}}$ and
${\bm{a}}$. Let
$g:\mathbb{R}^{d}\to\mathbb{R}:{\bm{x}}\mapsto\sum_{c=1}^{2k}a_{c}\sigma(\langle{\bm{x}},{\bm{w}}_{c}\rangle)$
where ${\bm{w}}_{c}$ is the $c$-th column of ${\bm{W}}$. Then we have
$f({\bm{X}})=\sum_{v=1}^{n}g(({\bm{A}}{\bm{X}})_{v})$ where the subscript
represents the row of the corresponding matrix. Since $\rho\in
L^{2}(\mathcal{N})$, $g$ is in $L^{2}(\mathcal{N})$ and we can write uniquely:

$g({\bm{x}})=\sum_{J\in\mathbb{N}^{d}}\hat{g}_{J}H_{J}({\bm{x}}).$ (82)

Therefore,

$f({\bm{X}})=\sum_{v=1}^{n}g(({\bm{A}}{\bm{X}})_{v})=\sum_{v=1}^{n}\sum_{J\in\mathbb{N}^{d}}\hat{g}_{J}H_{J}(({\bm{A}}{\bm{X}})_{v})=\sum_{J\in\mathbb{N}^{d}}(\hat{g}_{J}\sqrt{n})H_{J}^{{\bm{A}}}({\bm{X}})\in\mathcal{H}^{(d)}.$
(83)

∎

Now we verify that the special function $f_{{\mathcal{G}}}$ in our
construction of the family of hard functions has vanishing low degree moments:

######  Lemma 21 (Low degree moment vanishes).

For any $k\in\mathbb{N}$, for any ${\mathcal{G}}$ and appropriate ${\bm{A}}$,
$f_{{\mathcal{G}}}\in\mathcal{H}^{(2)}$. Furthermore,
$(\hat{f_{{\mathcal{G}}}})_{J}=0$ for every $|J|<k$.

###### Proof.

The first part follows from specializing Lemma 20 for $d=2$. The second part
follows from the fact that the corresponding $g$ function in the proof of
Lemma 20 for $f$ is the function $f_{\varphi,\sigma}$ in Diakonikolas et al.
[2020], whose low degree moments vanish. ∎

It is also easy to verify that $C_{{\mathcal{G}}}^{\mathcal{B}}\subset
F^{d,k}_{{\mathcal{G}}}$:

######  Lemma 22.

For any $k,d,n\in\mathbb{N}$, for some activation $\sigma$ that is not a low-
degree polynomial, $C_{{\mathcal{G}}}^{\mathcal{B}}\subset
F^{d,k}_{{\mathcal{G}}}$.

###### Proof.

Pick ${\bm{B}}\in\mathcal{B}$ arbitrarily and let
$h({\bm{X}})=f_{{\mathcal{G}}}({\bm{X}}{\bm{B}})$. We have:

$\displaystyle
h({\bm{X}})=f_{{\mathcal{G}}}({\bm{X}}{\bm{B}})=\bm{1}_{n}^{\top}\sigma\left({\bm{A}}({\bm{X}}{\bm{B}}){\bm{W}}^{*}\right){\bm{a}}^{*}=\bm{1}_{n}^{\top}\sigma\left({\bm{A}}{\bm{X}}({\bm{B}}{\bm{W}}^{*})\right){\bm{a}}^{*}.$
(84)

Since ${\bm{B}}{\bm{W}}^{*}$ is in $\mathbb{R}^{d\times 2k}$ and
${\bm{a}}^{*}$ is in $\mathbb{R}^{2k}$, we have $h\in
F^{d,k}_{{\mathcal{G}}}$. Functions in $\mathcal{F}_{k}^{\mathcal{B}}$ also do
not vanish, as long as ${\bm{A}}$ does not have vanishing rank and the
activation function $\sigma$ is not a low-degree polynomial, by the same
argument as in Diakonikolas et al. [2020]. ∎

From Diakonikolas et al. [2020], we know the following facts about matrices in
$\mathbb{R}^{d\times 2}$:

######  Lemma 23 (Lemma 4.7 in Diakonikolas et al. [2020]).

Fix $c\in(0,\frac{1}{2})$. Then there is a set $S\subset\mathbb{R}^{d\times
2}$ of size at least $2^{\Omega(d^{c})}$ such that for each
${\bm{A}},{\bm{B}}\in S$, $\|{\bm{A}}^{\top}{\bm{B}}\|_{2}\leq
O(d^{c-\frac{1}{2}})$.

The following Lemma generalizes Lemma 4.6 in Diakonikolas et al. [2020]:

######  Lemma 24.

Fix a function $\mathcal{H}^{(2)}\ni p:\mathbb{R}^{n\times 2}\to\mathbb{R}$
and ${\bm{U}},{\bm{V}}\in\mathbb{R}^{d\times 2}$ such that
${\bm{U}}^{\top}{\bm{U}}={\bm{V}}^{\top}{\bm{V}}=\mathbf{1}_{2}$. Then:

$\left|\langle
p(\cdot{\bm{U}})p(\cdot{\bm{V}})\rangle_{L^{2}(\mathcal{N})}\right|\leq\sum_{m=0}^{\infty}\left(2\sqrt{n}\|{\bm{U}}^{\top}{\bm{V}}\|_{2}\right)^{m}\cdot\|p^{[m]}\|^{2}_{L^{2}(\mathcal{N})},$
(85)

where $p^{[m]}$ is the degree $m$ part of $p$.

###### Proof.

Let $f,g:\mathbb{R}^{n\times d}\to\mathbb{R}$ be defined as
$f({\bm{X}})=p({\bm{X}}{\bm{U}})$ and $g({\bm{X}}):=p({\bm{X}}{\bm{V}})$. Let
$f^{[m]},g^{[m]}$ and $p^{[m]}$ be the degree $m$ part of the corresponding
function.

From Corollary 19,

$\displaystyle\sum_{m=0}^{\infty}\frac{1}{m!}\langle\nabla^{m}f^{[m]}({\bm{X}}),\nabla^{m}g^{[m]}({\bm{X}})\rangle$
$\displaystyle=\sum_{m=0}^{\infty}\mathbb{E}_{{\bm{X}}\sim\mathcal{N}^{n\times