On the Hölder Stability of Multiset and Graph Neural Networks

## Abstract

Abstract Extensive research efforts have been put into characterizing and constructing maximally separating multiset and graph neural networks.
However, recent empirical evidence suggests the notion of separation itself doesn’t capture several interesting phenomena. On the one hand, the quality of this separation may be very weak, to the extent that the embeddings of "separable" objects might even be considered identical when using fixed finite precision. On the other hand, architectures which aren’t capable of separation in theory, somehow achieve separation when taking the network to be wide enough. In this work, we address both of these issues, by proposing a novel pair-wise separation quality analysis framework which is based on an adaptation of Lipschitz and Hölder stability to parametric functions. The proposed framework, which we name Hölder in expectation , allows for separation quality analysis, without restricting the analysis to embeddings that can separate all the input space simultaneously. We prove that common sum-based models are lower-Hölder in expectation, with an exponent
that decays rapidly with the network’s depth . Our analysis leads to adversarial examples of graphs which can be separated by three 1-WL iterations, but cannot be separated in practice by standard maximally powerful Message Passing Neural Networks (MPNNs). To remedy this, we propose two novel MPNNs with improved separation quality, one of which is lower Lipschitz in expectation. We show these MPNNs can easily classify our adversarial examples, and compare favorably with standard MPNNs on standard graph learning tasks.

## 1 Introduction

Motivated by a multitude of applications, including molecular systems (Gilmer et al., ), social networks (Borisyuk et al., ), recommendation systems (Gao et al., ) and more, permutation invariant deep learning for both multisets and graphs have gained increasing interest in recent years. This in turn has inspired several theoretical works analyzing common permutation invariant models, and their expressive power and limitations.

For multiset data, it is known that simple summation-based multiset functions are injective, and as a result, can approximate all continuous functions on multisets (Zaheer et al., ). These results have been discussed and strengthened in many different recent publications (Wagstaff et al., ; Amir et al., ; Tabaghi & Wang, ; Wang et al., ).

For graph neural networks (GNNs) the situation is more delicate, as all known graph neural networks with polynomial complexity have limited expressive power. Our focus in this paper will be on the Message Passing Neural Network (MPNN) (Gilmer et al., ) paradigm, which includes a variety of popular GNNs (Xu et al., ; Gilmer et al., ; Kipf & Welling, ).
In their seminal works, Xu et al. (); Morris et al. () analyze the separation capabilities of MPNNs, showing that an MPNN $f$ can separate two graphs $G$ and $H$ (i.e., $f(G)\neq f(H)$) only if the Weisfeiler-Lehman (WL) isomorphism test (Weisfeiler & Lehman, ) can also do so.
Accordingly, maximally expressive MPNNs are those which are able to separate all graph pairs which are WL-separable. They further show that MPNNs which employ injective multiset functions are maximally expressive.

The theoretical ability of a permutation invariant network to separate a pair of objects (multisets/graphs) is a necessary condition for all learning tasks which require such separation, e.g., binary classification tasks where the two objects have opposite labels. However, while current theory ensures separation, it tells us little about the separation quality, so that embeddings of "separable" objects may be extremely similar, to the extent that in some cases, graphs which can be theoretically separated are completely identical on a computer with fixed finite precision. This pheonomena was observed for graphs with analytic activations and very small width, see (Bravo et al., ), figure 2 top.

In a non-parametric setting, separation quality of a function $f$ can be studied via bi-Hölder stability guarantees:
for metric spaces $(X,d_{X})$ and $(Y,d_{Y})$, a function $f:X\to Y$ is $\beta$ upper-Hölder and $\alpha$ lower-Hölder if there exist some positive constants $c,C$ such that

$$ $cd_{X}(x,x^{\prime})^{\alpha}\leq d_{Y}(f(x),f(x^{\prime}))\leq Cd_{X}(x,x^{ \prime})^{\beta},\forall x,x^{\prime}\in X.$ $$

The upper and lower Hölder conditions guarantee that embedding distances in $Y$ will not be much larger, or much smaller, than distances in the original space $X$ (in our case a space of multisets or graphs).
A function $f$ will have higher separation quality the closer the Hölder exponents are to one. When $\beta=1$ we say $f$ is (upper) Lipschitz, and when $\alpha=1$ we say that $f$ is lower Lipschitz.

Bi-Lipschitz stability analysis is a central topic in the study of frames Balan (); Alexeev et al. (), phase retrieval Bandeira et al. (); Cheng et al. (), and several group-invariant learning scenarios Cahill et al. (; ). For multisets, Bi-Lipschitz stability results are known for multiset-functions based on max-filters Cahill et al. () or sorting Balan et al. (). However, this is not the case for more common sum-based multiset functions: Amir et al. () showed that *Relu-sum* multiset functions, which are based on summation of point-wise applied ReLU networks (see (1) below) are never even injective. Multiset functions which use *smooth activations* instead, are injective but not lower-Lipschitz. Accordingly, a remaining challenge, which we will address in this paper, is characterizing the lower-Hölder stability of sum-based multiset functions.

For MPNN architectures for graphs, there are even less stability results. The upper Lipschitz-ness of the MPNN GIN (Xu et al., ) and similar architectures was established in (Chuang & Jegelka, ), however lower Lipschitz or Hölder guarantees have not been addressed to date. Indeed, a recent survey (Morris et al., ) lists bi-Lipschitz guarantees for MPNNs as one of the future challenges for theoretical GNN research. We will address this challenge in this paper as well.

As a first step for our stability analysis, we need to establish how to extend notions of Hölder stability to parametric functions. One simple approach is requiring the parametric function to be *uniformly* Hölder, with the same exponent and constant regardless of the parameter choice. Indeed, we will show that all parameteric functions we consider in this paper are uniformly *upper*-Lipschitz.

In contrast, we believe the notion of uniform *lower*-Hölder to be too stringent. For example, as mentioned above, multiset networks based on ReLU activations are never injective, and so can’t be uniformly lower-Hölder. Nonetheless, numerical estimates of the stability of such networks show that wide ReLU networks have comparable or even better stability than injective smooth-activation-based multiset functions (see Amir et al. () Figure 2). Accordingly, we resort to a probabilistic framework of *Hölder stability in expectation*, as presented in section 2.

Figure: Figure 1: Separation quality on a single adversarial multiset-pair constructed as described in appendix G.3.
Refer to caption: extracted/5914493/figs/relu_vs_smooth_per_separation.png

##### Main results: multisets

With respect to our new relaxed notion of expected lower-Hölder, we show that ReLU-sum multiset networks have an expected lower-Hölder exponent of $\alpha=3/2$. Surprisingly, while smooth activations lead to injectivity, we find that their expected Hölder exponent is much worse: at best it is equal to the maximal number of elements in the multiset data, $\alpha\geq n$. This scenario is summarized in Figure 1: For a given pair of multisets, smooth activations will separate even with a very small network width, but the separation quality can be very poor. In contrast, shallow Relu-sum networks may not attain separation. However, as width increases separation will be obtained (with high probability), and the expected quality of separation will surpass the quality attained with smooth activations.

The relatively moderate exponent $\alpha=3/2$ of ReLU networks is guaranteed only when the range of the bias of the network contains the range of the multiset features, an assumption which may be problematic to fulfill in practice. To address this, we suggest an *Adaptive ReLU* summation network, where the bias is adapted to the values in the multiset. This network is guaranteed the same $3/2$ exponent while overcoming the range issue.

Finally, we consider multiset functions based on linear functions and column-wise sorting. These were shown to be bi-Lipschitz, when wide enough, in (Balan et al., ). We show that they are lower Lipschitz *in expectation* even with a width of 1.

**Table 1: Summary of lower-Hölder exponent bounds for multiset (top) and graph (bottom) models.**
| Model | bounds |
| --- | --- |
| relu sum | $\alpha=3/2$ |
| adaptive relu | $\alpha=3/2$ |
| smooth sum | $n\leq\alpha$ |
| sort based | $\alpha=1$ |
| ReluMPNN | $1+\frac{K+1}{2}\leq\alpha$ |
| SmoothMPNN | $2^{K+1}\leq\alpha$ |
| SortMPNN | $\alpha=1$ |

##### Main Results: MPNNs

MPNNs are constructed upon multiple application of multiset functions to node neighborhoods. We consider four MPNNs based on the four multiset functions we analyzed: SortMPNN, ReluMPNN, SmoothMPNN and AdaptMPNN. We show that SortMPNN is lower-Lipschitz in expectation, even with a width of 1. SortMPNN is thus the first MPNN with both upper and lower Lipschitz guarantees. In contrast, we show that ReluMPNN and SmoothMPNN are only lower-Hölder, with an exponent that deteriorates as the MPNN depth grows. The exponent bounds of both the graph and multiset models appear in table 1.

Our analysis provides an adversarial example illustrating the deterioration of the expected lower-Hölder exponent of sum-based MPNNs as the depth increases. We show in Table 2 that sum-based MPNNs fail to learn a simple binary classification task on such data, while SortMPNN and AdaptMPNN handle this task easily. We also provide experiments on several graph datasets which show that SortMPNN often outperforms standard MPNNs on non-adversarial data, and is more robust to reduction in model size.

#### 1.0.1   Notation

$S^{d-1}$ denotes the unit sphere $\{x\in\mathbb{R}^{d}|\quad\|x\|_{2}=1\}$. The notation ’$a\sim S^{d-1}$ and $b\sim[-B,B]$’ implies that the distribution on $a$ (respectively $b$) is taken to be uniform on $S^{d-1}$ (respectively $[-B,B]$), and that $a$ and $b$ are drawn independently. The inner product of two vectors $a,b$ is denoted by $a\cdot b$.

## 2 General framework

In this section, we present the framework for analyzing separation quality of parametric functions. We begin by extending the notion of Hölder stability to a family of parametric functions.

### 2.1 Hölder stability for parametric functions

Let $(X,d_{X})$ and $(Y,d_{Y})$ be psuedo-metric(^1^11A psuedo metric $d_{X}$ differs from a standard metric only in that the $d_{X}(x_{1},x_{2})=0$ is possible even if $x_{1}\neq x_{2}$.)spaces. , let $W$ be some set of parameters, and let $f(x;w):X\times W\to Y$ be a parametric function. we say that $f(x;w)$ is *uniformly* Lipschitz with constant $L>0$, if for all $w\in W$, the
function $x\mapsto f(x;w)$ is $L$-Lipschitz. This definition can naturally be extended to other upper-Hölder exponents, but this will not be necessary as the parametric functions we discuss in this paper will all end up being uniformly Lipschitz. Similarly, we could define a notion of uniform lower-Hölder, but as discussed in this introduction, this notion is too stringent for the problems we are discussing, and so we introduce the alternative notion of lower-Hölder *in expectation*:

###### Definition 2.1 .

Let $\alpha>0,p\geq 1$. Let $(X,d_{X})$ and $(Y,d_{Y})$ be psuedo-metric spaces, and let $(W,\mathcal{B},\mu)$ be a probability space. Let $f(x;w):X\times W\to Y$ be a parametric function, so that for every fixed $x\in X$ the function $f(x;\cdot)$ is measurable.

We say that $f$ is $\alpha$ lower-Hölder *in expectation* if there exists a strictly positive
$c$ such that

$$ $c^{p}\leq\mathbb{E}_{w\sim\mu}\left\{\frac{d_{Y}(f(x;w),f(x^{\prime};w))}{d_{X }(x,x^{\prime})^{\alpha}}\right\}^{p},\quad\forall(x,x^{\prime})\in X\times X, \text{ with }d(x,x^{\prime})>0$ $$

A few remarks on this definition are in order. Firstly, we note that if the set of $w$ for which $f(x;w)$ is $\alpha$ lower-Hölder has positive probability, then $f(x;w)$ will be $\alpha$ lower-Hölder in expectation. The opposite is not true: Relu-sum networks for multisets are never injective, and hence never lower-Hölder, but we will show that they are lower-Hölder in expectation. This is possible since it is defined in a pairwise sense, similar to the separation analysis in Morris et al. ().

Next, we note that our definition takes an expectation over the parameters, but requires uniform boundedness over $x,x^{\prime}$ pairs. This choice was made since data distribution is unknown, while the parameter distribution at initialization can be chosen by the algorithm. Furthermore, we conjecture that bad separation at initialization will be difficult to overcome during training. Supporting evidence for this claim will be shown later on (see table 2), where models with a large lower-Hölder in expectation exponent, fail to learn an adversarial binary classification task.

Thirdly, we note that as with the standard Hölder definitions, Hölder in expectation exponents that are closer to one will be regarded as having higher separation quality.

Fourthly, $p$ in the definition refers to the $\ell_{p}$ norm chosen for the feature spaces. Results stated in the introduction are for the default $p=2$.

Finally, while the expected value across random parameters is informative as is, the variance can potentially have a strong effect on the separation in practice. Luckily, our analysis will typically apply for models $f(x;w)$ of width $1$, and models with width $W$ can be viewed as stacking of $W$ multiple model instances. Wider models will have the same expected distortion as the width one model, but the variance will converge to zero, in a rate proportional to $1/W$.
For more details on this effect see Appendix B.1.

We now proceed to analyze Hölder stability for parametric models on multisets.

###### Definition 2.1 .

## 3 Multiset Hölder stability

Standing assumptions:
Throughout the rest of this section, we will say that $(n,d,p,\Omega,z)$ satisfy our standing assumptions, if $n,d$ are natural numbers,
$p\geq 1$, the set $\Omega$ is a compact subset of $\mathbb{R}^{d}$, and $z$ will be a point in $\mathbb{R}^{d}\setminus\Omega$.

A *multiset* $\{\!\!\{x_{1},\ldots,x_{k}\}\!\!\}$ is a collection of $n$ unordered elements where (unlike sets) repetitions are allowed. We denote by $\mathcal{S}_{\leq n}(\Omega)$ and $\mathcal{S}_{=n}(\Omega)$ the space of all multisets with at most $n$ elements (respectively, exactly $n$ elements), which reside in $\Omega$.

A popular choice of a metric on $\mathcal{S}_{=n}(\Omega)$ is the Wasserstein metric

$$ ${W}_{1}(\{\!\!\{x_{1},\ldots,x_{n}\}\!\!\},\{\!\!\{y_{1},\ldots,y_{n}\}\!\!\}) =\min_{\tau\in S_{n}}\sum_{j=1}^{n}\|x_{j}-y_{\tau(j)}\|_{p}$ $$

Following (Chuang & Jegelka, ), we extend the Wasserstein metric to the space $\mathcal{S}_{\leq n}(\Omega)$ by applying the augmentation map which adds the vector $z$ to a given multiset until it is of maximal size $n$:

$$ $\rho_{(z)}\left(\{\!\!\{x_{1},\ldots,x_{k}\}\!\!\}\right)=\{\!\!\{x_{1},\ldots ,x_{k},x_{k+1}=z,\ldots.x_{n}=z\}\!\!\}$ $$

This induces the augmented Wasserstein metric $W_{1}^{z}$ on $\mathcal{S}_{\leq n}(\Omega)$,

$$ $W_{1}^{z}(S_{1},S_{2})=W_{1}(\rho_{(z)}(S_{1}),\rho_{(z)}(S_{2})),\quad S_{1}, S_{2}\in\mathcal{S}_{\leq n}(\Omega)$ $$

Now that we have introduced the augmented Wasserstein as the metric we will use for multisets, we turn to analyze the multiset neural networks we are interested in.

A popular building block in architectures for multisets (Zaheer et al., ) ,as well as MPNNs (Xu et al., ; Schütt et al., ; Gilmer et al., ), is based on summation of element-wise neural network application. For a given activation $\sigma:\mathbb{R}\to\mathbb{R}$, we define a parametric function $m_{\sigma}$ on multisets in $\mathcal{S}_{\leq n}(\Omega)$, with parameters $(a,b)\in\mathbb{R}^{d}\oplus\mathbb{R}$, by

$$ $m_{\sigma}(\{\!\!\{x_{1},\ldots,x_{r}\}\!\!\};a,b)=\sum_{i=1}^{r}\sigma(a\cdot x _{i}-b).$ (1) $$

We note that we focus on $m_{\sigma}$ with scalar outputs, since the variant with vector output has the same Hölder properties, as discussed in the end of Section 2.

#### 3.0.1   ReLU summation

We set out to understand the expected Hölder stability of $m_{\sigma}$, starting from the case $\sigma=\mathrm{ReLU}$:

theoremReluMoments
For $(n,d,p,\Omega,z)$ satisfying our standing assumptions, assume that $a\sim S^{d-1}$ and $b\sim[-B,B]$. Then $m_{ReLU}$ is uniformly Lipschitz. Moreover,

- 1.
$m_{ReLU}$ is not $\alpha$ lower-Hölder in expectation for any $\alpha<\frac{p+1}{p}$.
- 2.
If $\|x\|<B$ for all $x\in\Omega$, then $m_{ReLU}$ is $\frac{p+1}{p}$ lower-Hölder in expectation.

##### Proof idea: the ± ϵ plus-or-minus italic-ϵ \pm\epsilon ± italic_ϵ example.

The following example, which we nickname the $\pm\epsilon$ example, gives a good intuition for the Hölder behavior of $m_{\mathrm{ReLU}}$. Denote $X_{\epsilon}=\{\!\!\{-\epsilon,\epsilon\}\!\!\}$ and $X_{2\epsilon}=\{\!\!\{-2\epsilon,2\epsilon\}\!\!\}$. The Wasserstein distance between these two multisets is proportional to $\epsilon$. Now, let us consider the images of these multisets under $m_{\mathrm{ReLU}}$, where we assume for simplicity that $d=1$ and $a=1$. Note that when $b>2\epsilon$ we will get $ReLU(x-b)=0$ for all $x$ in either sets, and so $m_{\mathrm{ReLU}}(X_{2\epsilon};1,b)=m_{\mathrm{ReLU}}(X_{\epsilon};1,b)$.
We will get a similar results when $b\leq-2\epsilon$. In this case, for every $x\geq-2\epsilon$ we have $\mathrm{ReLU}(x-b)=x-b$. Therefore, we will obtain that

$$ $m_{\mathrm{ReLU}}(X_{2\epsilon};1,b)=(-2\epsilon-b)+(2\epsilon-b)=(-\epsilon-b )+(\epsilon-b)=m_{\mathrm{ReLU}}(X_{\epsilon};1,b),\quad\forall b<-2\epsilon.$ $$

We conclude that $|m_{\mathrm{ReLU}}(X_{2\epsilon};1,b)-m_{\mathrm{ReLU}}(X_{\epsilon};1,b)|^{p}$ is zero for all $b$ outside the interval $(-2\epsilon,2\epsilon)$. Inside this interval of length $4\epsilon$, we will typically have $|m_{\mathrm{ReLU}}(X_{2\epsilon};1,b)-m_{\mathrm{ReLU}}(X_{\epsilon};1,b)|^{p}
\sim\epsilon^{p}$, so that the expectation over all $b$ will be proportional to $4\epsilon\cdot\epsilon^{p}\sim\epsilon^{p+1}$. To ensure that the ratio between $\epsilon^{p+1}$ and $W_{p}(X_{2\epsilon},X_{\epsilon})^{\alpha p}\sim\epsilon^{\alpha p}$ will
have a strictly positive lower bound as $\epsilon\rightarrow 0$
, we need to choose $\alpha\geq\frac{p+1}{p}$. The proof that $\alpha=\frac{p+1}{p}$ is actually enough, and several other details necessary to turning this argument into a full proof, are given in the appendix.

#### 3.0.2   Adaptive ReLU

To attain the $(p+1)/p$ exponent in theorem 3.0.1, we need to assume that the bias is drawn from an interval $[-B,B]$ which is large enough so that $\|x\|<B$ for all $x\in\Omega$. This assumption may be difficult to satisfy, especially in MPNNs where the features in intermediate layers are effected by previous parameter choices. Additionally, the $\pm\epsilon$ example shows that when $b$ is outside the range of the multiset features, $m_{\mathrm{ReLU}}$ is not effective. Motivated by these observations, we propose a novel parametric function for multisets, based on $\mathrm{ReLU}$, where the bias is automatically adapted to feature values. For a multiset $X=\{\!\!\{x_{1},\ldots,x_{r}\}\!\!\}$ and $a\in\mathbb{R}^{d},t\in[0,1]$, the adaptive $\mathrm{ReLU}$ function $m_{\mathrm{ReLU}}^{\mathrm{adapt}}(X;a,t)$ is defined using $m=\min\{a\cdot x_{1},\ldots,a\cdot x_{r}\},M=\max\{a\cdot x_{1},\ldots,a\cdot x
_{r}\},b=(1-t)m+tM$ and

$$ $m_{\mathrm{ReLU}}^{\mathrm{adapt}}(X;a,t)=[r,m,M,\frac{1}{r}\sum_{i=1}^{r} \mathrm{ReLU}(a\cdot x_{i}-b)]$ $$

The output of $m_{\mathrm{ReLU}}^{\mathrm{adapt}}$ is a four dimensional vector . The last coordinate of the vector is essentially $m_{\mathrm{ReLU}}$, where the bias $b$ now only varies between the minimal and maximal value of multiset features. The first three coordinates are simple invariant features of the multisets: its cardinality and minimal and maximal value. In the appendix we prove
{restatable}theoremAdaptiveRelu
For $(n,d,p,\Omega,z)$ satisfying our standing assumptions, assume that $a\sim S^{d-1}$ and $t\sim[0,1]$. Then the function $m_{\mathrm{ReLU}}^{\mathrm{adapt}}:\mathcal{S}_{\leq n}(\Omega)\times S^{d-1}
\times[0,1]\to\mathbb{R}^{4}$ is uniformly Lipschitz and $\frac{p+1}{p}$ lower-Hölder in expectation.

We note that while adaptive ReLU solves the $\pm\epsilon$ example, a similar issue can be created for adaptive ReLU with multisets such as $\{\!\!\{1,0,0,-1\}\!\!\}$ and $\{\!\!\{1,\epsilon,-\epsilon,-1\}\!\!\}$, which force the adaptive bias to cover the interval $[1,1]$. Nonetheless, the delicate bias range assumption is no longer an issue.

#### 3.0.3   Summation with smooth activation

We now consider the Hölder properties of $m_{\sigma}$ when the activation $\sigma$ is smooth (e.g., sigmoid, SiLU, tanh). To understand this scenario, it can be instructive to first return to the $\pm\epsilon$ example. In this example, since the elements in $X_{\epsilon}$ and $X_{2\epsilon}$ both sum to the same number, one can deduce that (for any choice of bias) the first order Taylor expansion of $m_{\sigma}(X_{2\epsilon})-m_{\sigma}(X_{\epsilon})$ will vanish, so that this difference will be of the order of $\epsilon^{2}$. Based on this example we can already see that for smooth activations $\alpha\geq 2$. However, it turns out that there are adversarial examples with much worse behavior. Specifically, note that $X_{\epsilon}$ and $X_{2\epsilon}$ are ’problematic’ because their first moments are identical. In the same spirit, by choosing a pair $X,X^{\prime}$ of distinct multisets with $n$ elements in $\mathbb{R}$, whose first $n-1$ moments are identical, we can deduce that $\alpha\geq n$:

Figure: Figure 2: $l_{2}$ vs. $W_{2}$ distance on multiple adversarial multiset-pairs. Results are in accordance with our theoretical results (see Table1)
Refer to caption: extracted/5914493/figs/multiset_exponents.png

theoremSmoothMoments
Assume $(n,d,p,\Omega,z)$ satisfy our standing assumptions,
$\sigma:\mathbb{R}\mapsto\mathbb{R}$ has $n$ continuous derivatives, and $a\sim S^{d-1},b\sim[-B,B]$. If the function $m_{\sigma}$ is $\alpha$ lower-Hölder in expectation, then $\alpha\geq n$.

#### 3.0.4   Sort based

We now present the fourth and final multiset parametric function which we would like to address. Based on ideas from (Balan et al., ; Dym & Gortler, ), we consider the parametric function

$$ $S_{z}(X;a,b)=b\cdot\mathrm{sort}(a\cdot\rho_{(z)}(X))$ (2) $$

where $X\in\mathbb{R}^{d\times n}$ and $a,b$ are in $S^{d-1}$ and $S^{n-1}$ respectively. Note $S_{z}$ explicitly includes the augmentation by $z$, while in our previous discussions this was not necessary for expected Hölder stability with respect to $W_{p}^{z}$. On the bright side, the sort-based parametric functions are lower-Lipschitz in expectation.

theoremSortHolder
For $(n,d,p,\Omega,z)$ satisfying our standing assumptions. Assume that $a\sim S^{d-1}$ and $b\sim S^{n-1}$. Then $S_{z}(\cdot;a,b)$ is uniformly upper Lipschitz and lower Lipschitz in expectation.

The proof of this claim is that such mappings were shown to be bi-Lipschitz for all parameters, when duplicated enough times (with independent parameters) (Balan et al., ). This implies bi-Lipschitzness in expectation for the duplicated function, which in turn implies bi-Lipschitzness in expectation for a single function since the relevant expectations are identical.

##### Summary

In Figure 2, we choose a pair of distinct multisets $X,X^{\prime}$ with $16$ scalar features each, which have $15$ identical moments (for details on how this was constructed see Subsection C.2 in the appendix). The figure plots the distance between embedding distances versus the Wasserstein distance, for pairs $\epsilon X,\epsilon X^{\prime}$, for varying values of $\epsilon$. The figure illustrates that the $\alpha=3/2$ Hölder exponent of ReLU (we take $p=2$) and the $\alpha\geq n$ exponent of sigmoid summation are indeed encountered in this example. The sort and adaptive ReLU methods display a linear plot in this case (recall that while sort is indeed lower-Lipschitz in expectation, adaptive ReLU is not, but showing for adaptive ReLU will require different adversarial examples).

## 4 MPNN Hölder Stability

We will now analyze the expected Hölder stability of MPNNs for graphs, using the various multiset embeddings discussed previously. We begin by defining MPNNs, WL tests, and the graph WL-metric of our choice.

### 4.1 MPNNs and WL

We denote a graph by $G=(V,E,X)$ where $V$ are the nodes, $E$ are the edges, and $X=(x_{v})_{v\in V}$ are the initial node features. Throughout this section we will assume, as in the previous section, that $\Omega$ is a compact subset of $\mathbb{R}^{d}$, and $z$ is some fixed point in $\mathbb{R}^{d}\setminus\Omega$. We denote the neighborhood of a node $v$ by $\mathcal{N}_{v}=\{\!\!\{u\in V|(u,v)\in E\}\!\!\}$. MPNNs iteratively use the graph structure to redefine node features in the following manner:
initialize
$x_{v}^{(0)}=x_{v}$.

For $\mathbf{k=1\ldots,K}$

$$ $\displaystyle\text{{AGGREGATE}: }c_{v}^{(k)}=\phi^{(k)}(\{\!\!\{x_{u}^{(k-1)}| u\in\mathcal{N}_{v}\}\!\!\})$ $\displaystyle\text{{COMBINE}: }x_{v}^{(k)}=\psi^{(k)}(x_{v}^{(k-1)},c_{v}^{(k)})$ $$

In order to achieve a final graph embedding, a final READOUT step is performed

$$ $\textit{READOUT}:c^{\mathrm{global}}=\eta(\{\!\!\{x_{v}^{(K)}|v\in V\}\!\!\})$ $$

We note that different instantiations of the COMBINE, AGGREGATE, and READOUT parametric functions will lead to different architectures. We will discuss several choices, and their expected Hölder stability properties, slightly later on. Once such an instantiation is chosen, an MPNN will be a parametric model of the form $f(G;w)=c^{\mathrm{global}}$, where $w$ denotes the concatenation of all parameters of the COMBINE, AGGREGATE, and READOUT functions.

The Weisfeiler-Lehman (WL) graph isomorphism test (Weisfeiler & Lehman, ; Huang & Villar, ) checks whether two graphs $G,H$ are isomorphic (identical up to node permutation) by applying an MPNN-like procedure while choosing the COMBINE, AGGREGATE and READOUT functions to be hash functions, and running the procedure for at most $K=max\{|V_{G}|,|V_{H}|\}$ (Morris et al., ) iterations. The WL test can separate many, but not all, pairs of non-isomorphic graphs.

Notably, MPNNs can only separate pairs of graphs which are WL separable, and therefore if we hope for MPNN Hölder stability in expectation, we need to choose a *WL metric*: a pseudo-metric on graphs which is zero if and only if the pair of graphs cannot be separated by WL. Of the several possible choices in the recent literature for such a metric
(Grohe, ; Böker, ; Chen et al., ), we choose the family of Tree Mover’s Distance ($TMD^{(K)}$) suggested in (Chuang & Jegelka, ), which is described in detail in Appendix E. $TMD^{(K)}$ is a distance which is non-zero if and only if two graphs can be separated by $K$ iterations of the WL test. It is based on recursive application of the augmented Wasserstein to multisets of WL computation trees. We note that when $K$ is large enough, $TMD^{(K)}$ is a WL metric as discussed above.

### 4.2 MPNN stability analysis

We now consider the Hölder stability of MPNNs. In this setting, the graph domain we consider are graphs in $G_{\leq n}(\Omega)$, the collection of graphs with up to $n$ nodes, and node features in a compact set $\Omega\subseteq\mathbb{R}^{d}$.

It is natural to expect that the stability properties of an MPNN will be closely related to the AGGREGATE, READOUT and COMBINE functions composing the MPNN.
As a COMBINE function, we will choose some parametric functions which is uniformly upper-Lipschtiz, and lower Lipschitz in expectation. As the COMBINE function is a vector-to-vector function, this is easy to achieve even via a linear function $x_{v}^{(k)}=Ax_{v}^{(k-1)}+Ba_{v}^{(k)}.$
In Appendix D we prove the Lipschitz properties of this and three other COMBINE choices.

The AGGREGATE and READOUT functions are multiset to vector functions, and we have discussed the Hölder properties of four different such choices in Section 3. In practice, at each iteration $k$ we take $W_{k}$ parallel copies of these functions with independent parameters ($W$ corresponds to the width of the MPNN). As discussed in the end of Subsection 2.1, this maintains the same results in expectation while reducing the variance. For example, the SortMPNN architecture will be defined using the sort-based parametric functions $S_{z}$ from (2) via

$$ $\displaystyle\text{{AGGREGATE}: }c_{v,i}^{(k)}=S_{z_{k}}\left(\{\!\!\{x_{u}^{( k-1)}|u\in\mathcal{N}_{v}\}\!\!\};a_{i}^{(k)},b_{i}^{(k)}\right),\quad i=1, \ldots W_{k}$ $\displaystyle\text{{COMBINE}: }x_{v}^{(k)}=A^{(k)}x_{v}^{(k-1)}+B^{(k)}c_{v}^{ (k)}$ $\displaystyle\textit{READOUT}:c^{\mathrm{global}}_{i}=S_{z_{K+1}}\left(\{\!\! \{x_{v}^{(K)}|v\in V\}\!\!\};a_{i}^{(K+1)},b_{i}^{(K+1)}\right),\quad i=1, \ldots,W_{K+1}$ $$

Similarly, we will use the terms AdaptMPNN, ReluMPNN and SmoothMPNN to denote the MPNN obtained by replacing $S_{z}$ with the appropriate functions $m_{\mathrm{ReLU}}^{\mathrm{adapt}}$, $m_{\mathrm{ReLU}}$ or $m_{\sigma}$ with a smooth $\sigma$.

Our first result regards the upper Lipschitz bound.

theoremLipschitzMPNN (Uniformly Lipschitz MPNN embeddings, informal version)
Let $f:G_{\leq n}(\Omega)\rightarrow\mathbb{R}^{m}$ be an MPNN with $K$ layers. If the functions used for the aggregation $\phi^{(k)}$, combine $\psi^{(k)}$, and readout $\eta$ are all uniformly upper Lipschitz, then $f$ is uniformly upper Lipschitz with respect to $TMD^{(K)}$. In particular, ReluMPNN, SmoothMPNN, AdaptMPNN and SortMPNN are all uniformly upper Lipschitz.

It would seem natural
to expect that similar results will hold for lower-Hölder in expectation guarantees: namely, that if the AGGREGATE, COMBINE and READOUT functions used in the MPNN are lower-Hölder in expectation, then the overall MPNN will be lower-Hölder in expectation as well. (^2^22Indeed, we made this claim in a previous version. We thank Soutrik Sarangi and Abir De for pointing out the error in the proof.)
Unfortunately, this isn’t always the case. In Appendix F.2 we give an example of a pair of graphs which cannot be separated by any choice of parameter of a Relu-MPNN with width $1$ and depth $2$, even though they are separated by $2$ iterations of WL.

The argument above does not rule out the possibility that a wider Relu-MPNN will be lower-Hölder in expectation. Indeed, we suspect that this will be the case for an appropriate width $W=W(d,n)$, but we leave a formal proof of this result for future work. We do prove that, for arbitrarily wide ReluMPNN and SmoothMPNN, we cannot obtain a very good Hölder exponent, with a bound that deteriorates as the depth increases:

theoremlb
Assume that ReluMPNN with depth $K$ is $\alpha$ lower-Hölder in expectation with respect to $TMD^{(K)}$, then $\alpha\geq 1+\frac{K+1}{p}$. If SmoothMPNN with depth $K$ is $\alpha$ lower-Hölder in expectation then $\alpha\geq 2^{K+1}$.
The proof of this theorem is based on an adversarial set of examples we call $\epsilon$-Trees. These are defined recursively, where the first set of trees are of height two, and the leaves contain the $\pm\epsilon$ multisets. Deeper trees are then constructed by building upon substructures from the trees in the previous step as depicted in figure 7 alongside a rigorous formulation and proof in Appendix F.3.

In contrast to the previously discussed methods, SortMPNN is lower Lipschitz in expectation, even with a width of $1$.

theoremThmSortMPNNInformal(informal)
For any given $W\geq 1,K\geq 0$, SortMPNN with width $W$ and depth $K$ is lower Lipschitz in expectation with respect to $TMD^{(K)}$.

While we don’t formally analyze the lower-Hölder properties of AdaptMPNN, we conjecture its worst-case behavior will be similar to ReluMPNN. However, in some settings it will have better stability. For example, for our adversarial $\epsilon$-Trees example AdaptMPNN has a Lipschitz-like behavior, as shown in Figure 3.

## 5 Experiments

In the following experiments, we evaluate the four architectures SortMPNN, AdaptMPNN, ReluMPNN and SmoothMPNN. As the last two architectures closely resemble standard MPNN like GIN (Xu et al., ) with ReLU/smooth activation, our focus is mainly on the SortMPNN and AdaptMPNN architectures.
In our experiments we consider several variations of these architectures which were omitted in the main text for brevity, and are described in appendix G alongside further experiment details.

##### ϵ italic-ϵ \epsilon italic_ϵ -Tree dataset

**Table 2: $\epsilon$-Tree binary classification results**
| Model | Accuracy |
| --- | --- |
| GIN(Xu et al., ) | 0.5 |
| GCN(Kipf & Welling, ) | 0.5 |
| GAT(Velickovic et al., ) | 0.5 |
| ReluMPNN | 0.5 |
| SmoothMPNN | 0.5 |
| SortMPNN | 1.0 |
| AdaptMPNN | 1.0 |

In order to show that the expected (lower) Hölder exponent is indeed a good indicator of separation quality, and to further validate the importance of separation quality analysis, we first focus on the adversarial $\epsilon$-tree construction used to prove Theorem 4.2.

To begin with, we show how randomly initialized MPNNs with the different multiset embeddings we analyzed distort the TMD metric on the $\epsilon$-trees.
In figure 3 we plot the distance between the embeddings provided by the MPNNs, as a function of the parameter $\epsilon$ determining the $\epsilon$ trees (this parameter is proportional to the TMD distance between the trees). We see how the lower-Hölder exponent increases with depth for ReluMPNN and SmoothMPNN, in a manner which is consistent with the $\alpha=1+\frac{K+1}{p}$ and $2^{K+1}$ bounds suggested by Theorem 4.2 (for $p=2$). SortMPNN, in contrast, displays virtually no distortion in this example, and its behavior is consisten with the bi-Lipschitness predicted by our theory. Similar results (without theoretical justification) can be observed for AdaptMPNN.

We then proceed to train a binary classifier using these MPNNs in an attempt to separate the tree pairs.

Figure: Figure 3: $l_{2}$ vs. TMD on $\epsilon$-Trees. The targeted ReluMPNN and SmoothMPNN exponents deteriorate with depth, in accordance with our theory (see Table 1)
Refer to caption: extracted/5914493/figs/epsilon_trees_scatter.png

The results in Table 2 show how SortMPNN and AdaptMPNN achieve perfect performance, while ReluMPNN and SmoothMPNN, in addition to several baseline MPNNs completely fail. This serves as proof to the importance of separation quality analysis, especially in light of the fact smooth activation MLP moments (Amir et al., ) and GIN (Xu et al., ) are in theory capable of separation.

##### TUDataset

While the $\epsilon$-Tree dataset emphasizes the importance of high separation quality, validating our analysis, it is not obvious what effect this has in real world datasets, where we don’t necessarily know how the exponent behaves.
Therefore, we test SortMPNN and AdaptMPNN on a subset of the TUDatasets (Morris et al., ), including Mutag, Proteins, PTC, NCI1 and NCI109. Table 3 shows that SortMPNN and AdaptMPNN outperform several baseline MPNNs (results are reported using the evaluation method from (Xu et al., )).
Note that in all experiments we only compare to Vanilla MPNNs (1-WL), and not more expressive architectures which often do reach better results, at a higher computational cost.

**Table 3: Classification accuracy on a variety of TUdatasets (Morris et al., ). Best in Bold, second underlined.**
| Dataset | Mutag | Proteins | PTC | NCI1 | NCI109 |
| --- | --- | --- | --- | --- | --- |
| GIN(Xu et al., ) | 89.4$\pm$5.6 | 76.2$\pm$2.8 | 64.6$\pm$7 | 82.7$\pm$1.7 | 82.2$\pm$1.6 |
| GCN(Kipf & Welling, ) | 85.6$\pm$5.8 | 76$\pm$3.2 | 64.2$\pm$4.3 | 80.2$\pm$2.0 | NA |
| GraphSage(Hamilton et al., ) | 85.1$\pm$7.6 | 75.9$\pm$3.2 | 63.9$\pm$7.7 | 77.7$\pm$1.5 | NA |
| SortMPNN | 90.99$\pm$6.2 | 76.46$\pm$3.68 | 66.31$\pm$6.73 | 83.55$\pm$1.82 | 82.75$\pm$1.60 |
| AdaptMPNN | 90.41$\pm$6.1 | 75.12$\pm$3.64 | 66.87$\pm$5.37 | 82.77$\pm$1.72 | 83.26$\pm$0.86 |

##### LRGB

We further evaluate our architectures on two tasks from LRGB (Dwivedi et al., ): peptides-func, which is a multi-label graph classification task, and the peptides-struct regression task. Results in table 4 show that SortMPNN and AdaptMPNN outperform other methods on peptides-func. For peptides-struct, GCN is the best performing method with SortMPNN a close second.

**Table 4: LRGB results. Best in bold, Second underlined.**
| Dataset | peptides-func (AP$\uparrow$) | peptides-struct (MAE$\downarrow$) |
| --- | --- | --- |
| GINE(Hu* et al., ) | 0.6621$\pm$0.0067 | 0.2473$\pm$0.0017 |
| GCN(Kipf & Welling, ) | 0.6860$\pm$0.0050 | 0.2460$\pm$0.0007 |
| GatedGCN(Bresson & Laurent, ) | 0.6765$\pm$0.0047 | 0.2477$\pm$0.0009 |
| SortMPNN | 0.6940$\pm$0.0049 | 0.2464$\pm$0.0024 |
| AdaptMPNN | 0.6934$\pm$0.0099 | 0.2484$\pm$0.0034 |

Both experiments above adhere to the 500K parameter constraint as in (Dwivedi et al., ). In addition, we reran the peptides-struct experiment with a 100K, 50K, 25K, 7K and 1K parameter budget. The results in figure 4 show SortMPNN and AdaptMPNN outperforming GCN for smaller models, with SortMPNN achieving the best results. We believe this is related to the fact that even a SortMPNN with width=1 is lower-Lipschitz in expectation, and accordingly when only a small number of features is available, its advantage on other methods is more substantial.

Figure: Figure 4: pep-struct, #params vs. MAE
Refer to caption: extracted/5914493/figs/pep_struct_small.png

##### Subgraph aggregation networks - Zinc12K

Finally, we experiment the use of our methods as the backbone MPNN in a more advanced equivariant subgraph aggregation network (ESAN) (Bevilacqua et al., ), which has greater separation power than MPNNs.
To this extent, we run the experiment from (Bevilacqua et al., ) on the ZINC12K dataset, where we swap the base encoder from GIN (Xu et al., ) to SortMPNN and AdaptMPNN.

**Table 5: ESAN (Bevilacqua et al., ) MAE per base encoder on ZINC12K. Best in bold, second underlined.**
| Method | GIN (Xu et al., ) | SortMPNN | AdaptMPNN |
| --- | --- | --- | --- |
| DS-GNN (ED) | 0.172$\pm$0.008 | 0.157$\pm$0.007 | 0.176$\pm$0.008 |
| DS-GNN (ND) | 0.171$\pm$0.010 | 0.152$\pm$0.009 | 0.168$\pm$0.008 |
| DS-GNN (EGO) | 0.126$\pm$0.006 | 0.104$\pm$0.004 | 0.127$\pm$0.007 |
| DS-GNN (EGO+) | 0.116$\pm$0.009 | 0.115$\pm$0.008 | 0.126$\pm$0.007 |
| DSS-GNN (ED) | 0.172$\pm$0.005 | 0.169$\pm$0.004 | 0.173$\pm$0.007 |
| DSS-GNN (ND) | 0.166$\pm$0.004 | 0.167$\pm$0.006 | 0.167$\pm$0.008 |
| DSS-GNN (EGO) | 0.107$\pm$0.005 | 0.115$\pm$0.010 | 0.126$\pm$0.007 |
| DSS-GNN (EGO+) | 0.102$\pm$0.003 | 0.121$\pm$0.005 | 0.131$\pm$0.006 |

The results shown in table 5 show that SortMPNN outperforms GIN in five of the eight different scenarios. For fair comparison, the models don’t surpass the 100K parameter budget.

##### Code and timing

Timing for the LRGB experiment are provided in table 7 in the appendix. SortMPNN is marginally slower than GCN, and AdaptMPNN is X1.87 times slower. Code is available at (^3^33https://drive.filen.io/d/9e745743-42b7-4cb7-ae3b-2b4653b612f9#bwcvyWoNTmbQU7493xzx8EhEXAy9vj5v) .

## 6 Related work

##### Sorting

Sorting based operations were used for permutation invariant networks on multisets in (Zhang et al., ), and as readout functions for MPNNs in (Balan et al., ; Zhang et al., ). To the best of our knowledge, our SortMPNN is the first network using sorting for aggregations, and the first MPNN with provable bi-Lipschitz (in expectation) guarantees.

##### Bi-Lipschitz stability

An alternative bi-Lipschitz embedding for multisets was suggested in (Cahill et al., ). This embedding has higher computational complexity then sort based embedding. Some additional examples of recent works on bi-Lipschitz embeddings include (Balan & Tsoukanis, ; Cahill et al., ; ).

We note that
(Böker et al., ) does provide lower and upper stability estimates for MPNNs with respect to a WL metric. These stability estimates are in an $\epsilon-\delta$ sense, which do not rule out arbitrarily bad Hölder exponents. Moreover, they only consider graphs without node features. On the other hand, their analysis is more general in that they consider graphs of arbitrary size, and their graphon limit.

## 7 Summary and Limitations

We presented expected Hölder stability analysis for functions based on ReLU summation, smooth activation summation, adaptive ReLU, and sorting. Our theoretical and empirical results suggest SortMPNN as a promising alternative to traditional sum-based MPNNs.

A computational limitation of SortMPNN is that it requires prior knowledge of maximal multiset sizes for augmentation. We believe that future work will reveal ways of achieving lower-Lipschitz architectures without augmentation. Other avenues of future work include analyzing Hölder properties with respect to other graph metrics other than TMD, and resolving some questions we left open such as upper bounds for the Hölder exponents of smooth activations.

Acknowledgements: N.D. and Y.D. are supported by Israeli Science Foundation grant no.
272/23. We thank Haggai Maron and Guy Bar-Shalom for their useful insights and assistance.

## Appendix A Definitions and notation

The following definitions and notation are used throughout the appendices.

- •
We use the function $F(x,x^{\prime},w)=F_{p,\alpha,f}(x,x^{\prime},w)$ to denote
$F(x,x^{\prime},w)=\left\{\frac{d_{Y}(f(x,w),f(x^{\prime},w))}{d_{X}(x,x^{
\prime})^{\alpha}}\right\}^{p}$
- •
We call a pair of multisets balanced, if the number of elements they contain are equal, and otherwise we call them unbalanced.

## Appendix B Hölder stability in expectation properties

In this section we fill in the details of some properties of lower-Hölder in expectation functions, which were discussed in Subsection 2.1.

### B.1 Reducing variance by averaging

Given $N\in\mathbb{N}$, we can extend $f(x;w):X\times W\to Y$ to a new parametric function $f_{N}:X\times W^{N}\to Y^{N}$ defined as

$$ $f_{N}(x;w_{1},\ldots,w_{N})=\left[f(x,w_{1}),\ldots,f(x,w_{N})\right]$ (3) $$

where the f measure on $W^{N}$ is the product measure, and the distance we take on $Y^{N}$ is

$$ $d_{Y_{N}}([y_{1},\ldots,y_{N}],[y_{1}^{\prime},\ldots,y_{N}^{\prime}])=\left\{ \frac{1}{N}\sum_{n=1}^{N}d(y_{n},y_{n}^{\prime})^{p}\right\}^{1/p}$ $$

The function $F_{N}$ corresponding to this choice of $f_{N}$ is

$$ $\displaystyle F_{N}(x,x^{\prime},w)$ $\displaystyle=\left\{\frac{d_{Y_{N}}(f_{N}(x,\mathbf{w}),f_{N}(x^{\prime}, \mathbf{w}))}{d_{X}(x,x^{\prime})^{\alpha}}\right\}^{p}$ $\displaystyle=\left(\frac{1}{d_{X}(x,x^{\prime})}\right)^{\alpha p}\frac{1}{N} \sum_{n=1}^{N}d(f(x,w_{n}),f(x^{\prime},w_{n}))^{p}$ $\displaystyle=\frac{1}{N}\sum_{n=1}^{N}F(x,x^{\prime},w_{n})$ $$

Thus $F_{N}$ is the average of $N$ independent copies of $F$, which implies that for all $x\neq x^{\prime}$ the random variable $F_{N}(x,x^{\prime},\mathbf{w})$ has the same expectation as $F$, and its standard deviation is uniformly bounded by $\sigma/\sqrt{N}$, where $\sigma$ is a bound on the standard deviation of $F(x,x^{\prime},w)$ which is uniform in $x,x^{\prime}$ (assuming that such a bound exists). In particular, this means that *for any* $x,x^{\prime}$, not only will the expectation of $F_{N}(x,x^{\prime},\mathbf{w})$ be bounded by from below by $m^{p}$, but also as $N$ goes to infinity the probability of $F_{N}(x,x^{\prime},\mathbf{w})$ being larger than $m^{p}-\epsilon$ will go to one.

### B.2 Boundedness assumptions

In the scenarios we consider in this paper the metric spaces are bounded. Assume $0<\alpha<\beta$, and the distances between any two elements $x,x^{\prime}$ in a metric space $X$ are bounded by some constant $B$, then

$$ $d(x,x^{\prime})^{\beta}=B^{\beta}\left(\frac{d(x,x^{\prime})}{B}\right)^{\beta }\leq B^{\beta}\left(\frac{d(x,x^{\prime})}{B}\right)^{\alpha}=B^{\beta-\alpha }d(x,x^{\prime})^{\alpha}.$ $$

In other words, in bounded metric spaces, up to a constant, increasing the exponent decreases the value. It follows that a parametric functions which is $\alpha$ lower Hölder in expectation is also $\beta$ lower Hölder in expectation.

### B.3 Composition

In this section we discuss what happens when we compose two functions $f(x;w):X\times W\mapsto Y$ and $g(y;v):Y\times V\mapsto Z$ to get a new parametric function

$$ $g\circ f(x;w,v)=g(f(x;w);v).$ $$

###### Lemma B.1 .

Given two parametric functions $f(x;w):X\times W\mapsto Y$ and $g(y;v):Y\times V\mapsto Z$ which are $\alpha$ and $\beta$ lower Hölder in expectation with $\alpha,\beta\geq 1$, then $g\circ f$ is $\alpha\cdot\beta$ lower Hölder in expectation. Similarly, if $f,g$ are uniformly upper Lipschitz, then $f\circ g$ is uniformly upper Lipschitz.

###### Proof.

We omit the proof for uniform upper Lipschitz, and prove only the more challenging case of lower-Hölder in expectation. We have that for strictly positive $m_{f},m_{g}$, $\mathbb{E}_{w\sim\mu}[\{\frac{d_{Y}(f(x;w),f(x^{\prime};w))}{d_{X}(x,x^{\prime
})^{\alpha}}\}^{p}]\geq m_{f}$ and $\mathbb{E}_{v\sim\gamma}[\{\frac{d_{Z}(g(y,v),g(y^{\prime},v))}{d_{Y}(y,y^{
\prime})^{\beta}}\}^{p}]\geq m_{g}$, from which we conclude

$$ $\mathbb{E}_{w,v\sim(\mu,\gamma)}[\{\frac{d_{Z}(g(f(x;w),v),g(f(x^{\prime};w),v ))}{d_{X}(x,x^{\prime})^{\alpha\cdot\beta}}\}^{p}]$ $$

$$ $\stackrel{{\scriptstyle\text{w,v are independent}}}{{=}}\int_{w}\int_{v}\{ \frac{d_{Z}(g(f(x;w),v),g(f(x^{\prime};w),v))}{d_{X}(x,x^{\prime})^{\alpha \cdot\beta}}\}^{p}d\gamma d\mu$ $$

$$ $=\int_{w}\{\frac{d_{Y}(f(x;w),f(x^{\prime};w))^{\beta}}{d_{X}(x,x^{\prime})^{ \alpha\cdot\beta}}\}^{p}\int_{v}\{\frac{d_{Z}(g(f(x;w),v),g(f(x^{\prime};w),v) )}{d_{Y}(f(x;w),f(x^{\prime};w))^{\beta}}\}^{p}d\gamma d\mu$ $$

$$ $=\int_{w}\{\frac{d_{Y}(f(x;w),f(x^{\prime};w))^{\beta}}{d_{X}(x,x^{\prime})^{ \alpha\cdot\beta}}\}^{p}\cdot\mathbb{E}_{v\sim\gamma}[\{\frac{d_{Z}(g(f(x;w),v ),g(f(x^{\prime};w),v))}{d_{Y}(f(x;w),f(x^{\prime};w))^{\beta}}\}^{p}]d\mu$ $$

$$ $\geq\int_{w}\{\frac{d_{Y}(f(x;w),f(x^{\prime};w))^{\beta}}{d_{X}(x,x^{\prime}) ^{\alpha\cdot\beta}}\}^{p}\cdot m_{g}d\mu$ $$

$$ $=m_{g}\cdot\mathbb{E}_{w\sim\mu}[(\{\frac{d_{Y}(f(x;w),f(x^{\prime};w))}{d_{X} (x,x^{\prime})^{\alpha}}\}^{p})^{\beta}]$ $$

$$ $\stackrel{{\scriptstyle(*)}}{{\geq}}m_{g}\cdot(\mathbb{E}_{w\sim\mu}[\{\frac{d _{Y}(f(x;w),f(x^{\prime};w))}{d_{X}(x,x^{\prime})^{\alpha}}\}^{p}])^{\beta}$ $$

$$ $\geq m_{g}\cdot m_{f}^{\beta}$ $$

Where $(*)$ is true by Jensen’s inequality, since $(\frac{d_{Y}(f(x;w),f(x^{\prime};w))}{d_{X}(x,x^{\prime})^{\alpha}})^{p}$ is non negative, and $\phi:\mathbb{R}^{+}\cup\{0\}\mapsto\mathbb{R}$, $\phi(t)=t^{\beta}$ is convex for $\beta\geq 1$
∎

###### Lemma B.1 .

###### Proof.

## Appendix C Multiset embeddings analysis proofs

The following are the proofs of the claims regarding the exponent of Hölder stability in expectation from section 3.

Before stating and proving the relevant claims, we will present some lemmas that will be used throughout this section and those that follow.

We first state the well known property of norm equivalence on finite dimensional normed spaces

###### Lemma C.1 .

Let $1\leq p,q\leq\infty$ and $d\in\mathbb{N}$. There exists $c_{p,q,d}>0$ s.t. for any $v\in\mathbb{R}^{d}$

$$ $\|v\|_{p}\geq c_{p,q,d}\cdot\|v\|_{q}$ $$

In particular, this means that

###### Lemma C.2 .

Let $1\leq p,q\leq\infty$ and $\Omega\subseteq\mathbb{R}^{d}$ There exists $C_{p,q,d}>0$ s.t. for any $X,Y\in\mathcal{S}_{\leq n}(\Omega)$

$$ $W_{p}^{(z)}(X,Y)\geq C_{p,q,d}\cdot W_{q}^{(z)}(X,Y)$ $$

###### Proof.

$$ $\displaystyle W_{p}^{(z)}(X,Y)$ $\displaystyle=\left\{\min_{\tau\in S_{n}}\sum_{j=1}^{n}\|\rho_{(z)}(x)_{j}- \rho_{(z)}(y)_{\tau(j)}\|_{p}^{p}\right\}^{1/p}$ $\displaystyle\stackrel{{\scriptstyle\ref{norm_equiv}}}{{\geq}}c_{p,1,d}\cdot \left\{\min_{\tau\in S_{n}}\sum_{j=1}^{n}\|\rho_{(z)}(x)_{j}-\rho_{(z)}(y)_{ \tau(j)}\|_{p}\right\}$ $\displaystyle\stackrel{{\scriptstyle\ref{norm_equiv}}}{{\geq}}c_{p,1,d}\cdot c _{p,q,d}\cdot\left\{\min_{\tau\in S_{n}}\sum_{j=1}^{n}\|\rho_{(z)}(x)_{j}-\rho _{(z)}(y)_{\tau(j)}\|_{q}\right\}$ $\displaystyle\stackrel{{\scriptstyle\ref{norm_equiv}}}{{\geq}}c_{p,1,d}\cdot c _{p,q,d}\cdot c_{q,1,d}\cdot\left\{\min_{\tau\in S_{n}}\sum_{j=1}^{n}\|\rho_{( z)}(x)_{j}-\rho_{(z)}(y)_{\tau(j)}\|_{q}^{q}\right\}^{\frac{1}{q}}$ $\displaystyle=c_{p,1,d}\cdot c_{p,q,d}\cdot c_{q,1,d}\cdot W_{q}^{(z)}(X,Y)$ $$

The upper bound can be proven in the same manner.
∎

Now that we are equipped with the above lemmas, we begin by proving for smooth functions.

###### Lemma C.1 .

###### Lemma C.2 .

###### Proof.

### C.1 Summing over a smooth activation

*

###### Proof.

It is sufficient to prove the claim in the case where $d=1$.

Since $\sigma$ is $n$ times continuously differentiable, there exists some constant $C$ such that

$$ $|\frac{\partial^{n}}{\partial t^{n}}\sigma(at+b)|\leq C,\forall(t,a,b)\in[-1,1 ]\times S^{d-1}\times[-B,B].$ $$

Now, let $f$ be any function with $n$ continuous derivatives and $|f^{(n)}(t)|<C$ for all $t\in[-1,1]$. For $x\in(-1,1)^{n}$ denote $F(x)=\sum_{i=1}^{n}f(x_{i})$. (when $f(t)=\sigma(at+b)$ we have $F(x)=m_{\sigma}(x;a,b)$). By Taylor’s approximation of $f$ around 0 0, there exists points $c_{1},\ldots,c_{n}\in(-1,1)$ such that

$$ $\displaystyle F(x)$ $\displaystyle=\sum_{i=1}^{n}\left(f(0)+f^{\prime}(0)x_{i}+\frac{f^{\prime \prime}(0)}{2}x_{i}^{2}+\ldots+\frac{f^{(n-1)}(0)}{(n-1)!}x_{i}^{n-1}+\frac{f^ {(n)}(c_{i})}{n!}x_{i}^{n}\right)$ $\displaystyle=nf(0)+f^{\prime}(0)\sum_{i=1}^{n}x_{i}+\frac{f^{\prime\prime}(0) }{2}\sum_{i=1}^{n}x_{i}^{2}+\ldots+\frac{f^{(n-1)}(0)}{(n-1)!}\sum_{i=1}^{n}x_ {i}^{n-1}+\sum_{i=1}^{n}x_{i}^{n}\frac{f^{(n)}(c_{i})}{n!}$ $$

It follows that if $x,y$ are two vectors which are not the same, even up to permutation, but their first $n-1$ moments are identical, then for every $\epsilon\in(0,1)$, the vectors $\epsilon x$ and $\epsilon y$ will also have the same $n-1$ moments, and so there exist $c_{1},\ldots,c_{n},d_{1},\ldots,d_{n}$ in $(-1,1)$ such that

$$ $|F(\epsilon x)-F(\epsilon y)|=\epsilon^{n}|\sum_{i=1}^{n}\frac{f^{(n)}(c_{i})} {n!}x_{i}^{n}+\sum_{i=1}^{n}\frac{f^{(n)}(d_{i})}{n!}y_{i}^{n}|\leq\frac{2n}{n !}C\epsilon^{n}$ $$

In contrast, the Wasserstein distance between $\epsilon x$ and $\epsilon y$ scales like $\epsilon$:

$$ $W_{p}^{z}(\epsilon x,\epsilon y)=W_{p}(\epsilon x,\epsilon y)=\epsilon W_{p}(x ,y)$ $$

It follows that for all $\alpha<n$, all parameters $w=(a,b)\in S^{d-1}\times[-B,B]$, and all $\epsilon\in(0,1)$,

$$ $\left|\frac{m_{\sigma}(\epsilon x;w)-m_{\sigma}(\epsilon y;w)}{W_{p}^{\alpha}( \epsilon x,\epsilon y)}\right|^{p}\leq CW_{p}(x,y)^{-\alpha}\frac{2n}{n!} \epsilon^{(n-\alpha)p}$ $$

Taking the expectation over the parameters $w$ on the left hand side, we get the same inequality. Taking the limit $\epsilon\rightarrow 0$, we obtain

$$ $\mathbb{E}_{w}\left|\frac{m_{\sigma}(\epsilon x;w)-m_{\sigma}(\epsilon y;w)}{W _{p}^{\alpha}(\epsilon x,\epsilon y)}\right|^{p}\leq CW_{p}(x,y)^{-\alpha} \frac{2n}{n!}\epsilon^{(n-\alpha)p}\rightarrow 0$ $$

from which we deduce that $m_{\sigma}$ is not $\alpha$ lower-Hölder in expectation.

To conclude this argument , we note that there always exist $x,y\in\mathbb{R}^{n}$ which are not identical up to permutation, but have the same first $n-1$ moments. This is because it is known that no mapping from $(-1,1)^{n}$ to $\mathbb{R}^{n}$, and in particular the mapping which takes a vector of length $n$ to its first $n-1$ moments, cannot be injective, up to permutations (Wagstaff et al., ). To ensure that $x,y$ are in $(-1,1)^{n}$, we can scale both vectors by a sufficiently small positive number. In the next subsection we discuss a constructive method for producing pairs of sets with $n-1$ equal moments.
∎

###### Proof.

### C.2 Creating sets with n − 1 𝑛 1 n-1 italic_n - 1 equal moments

In the proof of theorem 2 we
explained that for every natural $n$, there exists pairs
of vectors $x,y\in\mathbb{R}^{n}$ which are not identical, up
to permutation, but have the same first $n-1$ moments.
We now give a constructive algorithm to construct such
pairs for $n$ which is a power of $2$, that is, $n=2^{k}$
for some natural $k$. This algorithm was used to produce
figure 2 in the main text.

Our algorithm operates by recursion on $k$. For $k=1$, our goal is to find a pair of vectors $x^{(1)},y^{(1)}$ of length $2^{k}=2$, whose first $2^{k}-1=1$ moments are equal. We can choose for example

$$ $x^{(1)}=[-1,1],y^{(1)}=[-2,2].$ $$

Now, assume that for a given $k$ we have a pair of vectors $x=x^{(k)}$ and $y=y^{(k)}$ of length $n=2^{k}$ which are not identical, up to permutation, but whose first $2^{k}-1$ moments are identical. We then perform the following three steps:

step 1: We translate all coordinates of $x$ and $y$ by the same number $t$, which we take to be the minimum of all coordinates in $x$ and $y$. We thus get a new pair of vectors

$$ $x^{\prime}=[x_{1}-t,\ldots,x_{n}-t],\quad y^{\prime}=[x_{1}-t,\ldots,x_{n}-t]$ $$

whose coordinates are all non-negative. These vectors are still distinct, even up to permutation, and have the same first $n-1=2^{k}-1$ moments.

step 2:
We take the root of the non-negative entries in both vectors, to get new vectors

$$ $x^{\prime\prime}=[\sqrt{x_{1}^{\prime}},\ldots,\sqrt{x_{n}^{\prime}}],\quad y^ {\prime\prime}=[\sqrt{y_{1}^{\prime}},\ldots,\sqrt{y_{n}^{\prime}}].$ $$

These two vectors are distinct, even up to permutations, and they agree on the even moments $2,4,\ldots,2(n-1)$.

step 3:
Finally, we define a new pair of vectors of cardinality $2n=2^{k+1}$

$$ $x^{(k+1)}=[-x^{\prime\prime},x^{\prime\prime}],\quad y^{(k+1)}=[-y^{\prime \prime},y^{\prime\prime}]$ $$

These two new vectors are still distinct up to permutations. There even moments up to order $2(n-1)$ are still identical. Moreover, all odd moments of both vectors are zero, and hence identical. In particular, the two vectors have the same moments of order
$1,2,\ldots,2(n-1),2n-1$, which is what we wanted.

### C.3 Sorting

Our goal in this subsection is proving

*

Recall that $S_{z}$ is defined via

$$ $S_{z}(X;a,b)=b\cdot\mathrm{sort}(a\cdot\rho_{(z)}(X))$ $$

Once $\rho_{(z)}$ is applied, we can think of $S_{z}$ as operating on multisets with exactly $n$ elements, coming from the domain $\Omega\cup\{z\}$. The function $S_{z}$ is then a composition of parametric functions $L\circ s$, where

$$ $L(y;b)=b\cdot y$ (4) $$

and

$$ $s(X;a)=sort(a^{T}X)$ (5) $$

By the rules of composition (Lemma B.1), it is sufficient to show that both functions are uniformly Lipschitz, and lower-Lipschitz in expectation.

We first show for $L$

###### Lemma C.3 .

For every $p\geq 1$, the function $L:\mathbb{R}^{n}\times S^{n-1}\to\mathbb{R}$ is uniformly upper-Lipschitz, and lower Lipschitz in expectation.

###### Proof.

Uniformly upper Lipschitz
For every $b\in S^{n-1}$ and $x,x^{\prime}\in\mathbb{R}^{n}$ we have, using Cauchy-Schwarz,

$$ $|L(x;b)-L(x^{\prime};b)|=|b\cdot(x-x^{\prime}),|\leq\|x-x^{\prime}\|$ $$

so for all $b$ we have a Lipschitz constant of 1.

Lower Lipschitz in expectation
For every $x\neq x^{\prime}$ in $\mathbb{R}^{n}$ we have, due to lemma C.1, for the appropriate $c>0$

$$ $\displaystyle\mathbb{E}_{b\sim S^{n-1}}\left[\left\{\frac{|L(x;b)-L(x^{\prime} ;b)|}{\|x-x^{\prime}\|_{p}}\right\}^{p}\right]$ $\displaystyle\geq c^{p}\cdot\mathbb{E}_{b\sim S^{n-1}}\left[\left\{\frac{|L(x; b)-L(x^{\prime};b)|}{\|x-x^{\prime}\|_{2}}\right\}^{p}\right]$ $\displaystyle=c^{p}\cdot\mathbb{E}_{b\sim S^{n-1}}\left[|b\cdot\frac{x-x^{ \prime}}{\|x-x^{\prime}\|_{2}}|^{p}\right]$ $\displaystyle=c^{p}\cdot\mathbb{E}_{b\sim S^{n-1}}\left[|b\cdot e_{1}|^{p} \right]>0$ $$

where the last equality is because the distribution is rotational invariant.
∎

We now prove for $s$

###### Lemma C.4 .

Let $d,n,p$ be natural numbers. Let $s$ be as defined in (5), then $s$ is uniformly upper Lipschitz, and lower Lipschitz in expectation. (here, the domain of $s$ is the space $\mathcal{S}_{=n}(\Omega)$ of multisets with $n$ elements in a compact set $\Omega\subseteq\mathbb{R}^{d}$, endowed with the $W_{1}$ metric. $a$ is drawn uniformly from $S^{d-1}$, and the metric on the output of $s$ is the $\ell_{p}$ distance).

###### Proof.

Due to norm equivalence as stated in lemma C.1, it is sufficient to prove the claim when $p=2$.

Fix some balanced multisets $X\not\sim Y$. Let $\sigma$ be a permutation which minimizes $\sum_{j=1}^{n}\|x_{j}-y_{\tau(j)}\|_{2}$ over all $\tau\in S_{n}$. Then for every $a\in S^{d-1}$ we have, using Cauchy-Schwartz and norm equivalence C.1

$$ $\displaystyle\|s(X;a)-s(Y;a)\|_{2}^{2}$ $\displaystyle=\|sort(a^{T}X)-sort(a^{T}Y)\|_{2}^{2}$ $\displaystyle=\min_{\tau\in S_{n}}\sum_{j=1}^{n}|a^{T}x_{j}-a^{T}y_{\tau(j)}|^ {2}$ $\displaystyle\leq\sum_{j=1}^{n}|a^{T}x_{j}-a^{T}y_{\sigma(j)}|^{2}$ $\displaystyle\leq\sum_{j=1}^{n}\|x_{j}-y_{\sigma(j)}\|_{2}^{2}$ $\displaystyle\leq C\left(\sum_{j=1}^{n}\|x_{j}-y_{\sigma(j)}\|_{2}\right)^{2}$ $\displaystyle=W_{1}(X,Y)^{2}$ $$

Concluding $s$ is uniformly upper Lipschitz.

In addition, we also have that for large enough $N$, the function $s_{N}(\cdot;a)$ (concatenation of $s$ $N$ times, divided by $1/N$) will be lower Lipschitz for Lebesgue almost every $a\in S^{d-1}$, as proved in (Balan et al., ). In particular, it follows that $s_{N}$ is lower Lipschitz in expectation. As discussed in Subsection B.1, this implies that $s$ is lower Lipschitz in expectation as well.
∎

###### Lemma C.3 .

###### Proof.

###### Lemma C.4 .

###### Proof.

### C.4 Analyzing ReLU

*

###### Proof.

We divide the proof into three parts, in accordance with the three parts of the theorem.

Part 1: Uniform Lipschitz
We first prove a uniform Lipschitz bound for all balanced multisets $Y,Y^{\prime}$ of cardinality $k\leq n$. Denote $X=\rho_{(z)}(Y),X^{\prime}=\rho_{(z)}(Y^{\prime})$. Then for every permutation $\tau\in S_{n}$,

$$ $\displaystyle|m_{ReLU}(Y;a,b)-m_{ReLU}(Y^{\prime};a,b)|$ $\displaystyle=|m_{ReLU}(X;a,b)-m_{ReLU}(X^{\prime};a,b)|$ $\displaystyle=|\sum_{i=1}^{n}ReLU(ax_{i}-b)-ReLU(ax_{\tau(i)}^{\prime}-b)|$ $\displaystyle\stackrel{{\scriptstyle(*)}}{{\leq}}\sum_{i=1}^{n}|(ax_{i}-b)-(ax _{\tau(i)}^{\prime}-b)|$ $\displaystyle=\sum_{i=1}^{n}|a\cdot(x_{i}-x_{\tau(i)}^{\prime})|$ $\displaystyle\stackrel{{\scriptstyle(**)}}{{\leq}}\sum_{i=1}^{n}\|a\|_{2}\|x_{ i}-x_{\tau{(i)}}\|_{2}$ $\displaystyle=\sum_{i=1}^{n}\|x_{i}-x_{\tau{(i)}}\|_{2},$ $$

where (*) is because $ReLU$ is Lipschitz with constant 1, and (**) is from Cauchy-Schwartz.

Since the inequality we obtain holds for all permutations $\tau$, we can take the minimum over all permutations to obtain that

$$ $|m_{ReLU}(Y;a,b)-m_{ReLU}(Y^{\prime};a,b)|\leq W_{1}(X,X^{\prime})=W_{1}^{z}(Y ,Y^{\prime})$ $$

To address multisets $Y,Y^{\prime}$ of different sizes, we first note that since the elements of the multisets are in $\Omega$, and the parameters $a,b$ come from a compact set, there exists some constant $M>0$ such that

$$ $|m_{ReLU}(Y;a,b)-m_{ReLU}(Y^{\prime};a,b)|\leq M.$ $$

On the other hand, for all $Y,Y^{\prime}$ of different sizes, we will always have that

$$ $W_{1}^{z}(Y,Y^{\prime})\geq dist(z,\Omega)>0$ $$

therefore

$$ $|m_{ReLU}(Y;a,b)-m_{ReLU}(Y^{\prime};a,b)|\leq M\frac{W_{1}^{z}(Y,Y^{\prime})} {W_{1}^{z}(Y,Y^{\prime})}\leq\frac{M}{dist(z,\Omega)}W_{1}^{z}(Y,Y^{\prime})$ $$

Combining this with our bound for multisets of equal cardinality, we see that $m_{ReLU}$ is uniformly Lipschitz with constant $\max\{1,\frac{M}{dist(z,\Omega)}\}$.

Part 2: Lower bound on expected Hölder exponent

We show that $m_{ReLU}$ is not $\beta$ lower-Hölder in expectation for all $\beta<(p+1)/p$.

Let $X$ be some matrix in $\Omega^{n}$ whose first two columns are the same $x_{1}=x_{2}$. Let $q$ be a vector with unit norm. For every $\epsilon>0$ define

$$ $X_{\epsilon}=[x_{1}-\epsilon q,x_{1}+\epsilon q,x_{3},\ldots,x_{n}]$ $$

It is not difficult to see that for all small enough $\epsilon$ we have that $W_{1}(X_{\epsilon},X)=2\epsilon$. On the other hand, for every fixed $a\in S^{d-1}$ and $b\in[-B,B]$ we have that, denoting $y=a\cdot x_{1}$ and $\delta=|\epsilon a\cdot q|$, we have

$$ $\displaystyle m_{ReLU}(X;a,b)-m_{ReLU}(X_{\epsilon};a,b)$ $\displaystyle=2\mathrm{ReLU}(a\cdot x_{1}-b)-\mathrm{ReLU}(a\cdot(x_{1}- \epsilon q)-b)$ $\displaystyle\quad\quad-\mathrm{ReLU}(a(\cdot x_{1}+\epsilon q)-b)$ $\displaystyle=2\mathrm{ReLU}(y-b)-[\mathrm{ReLU}(y+\delta-b)+\mathrm{ReLU}(y- \delta-b)]$ $$

Note that if $b>y+\delta$ then the expression above will be zero because all arguments of the ReLUs will be negative, and if $b<y-\delta$ then the expression above will also be zero because the arguments of all ReLUs will be positive so that we obtain

$$ $2\mathrm{ReLU}(y-b)-[\mathrm{ReLU}(y+\delta-b)+\mathrm{ReLU}(y-\delta-b)]=2(y- b)-[y+\delta-b+y-\delta-b]=0$ $$

Thus this expression will not vanish only if $b\in[y-\delta,y+\delta]$ which is an interval of diameter $2\delta\leq 2\epsilon$. For each $b$ in this interval we have, since ReLU is 1-Lipschitz, that

$$ $\displaystyle|2\mathrm{ReLU}(y-b)-[\mathrm{ReLU}(y+\delta-b)+\mathrm{ReLU}(y- \delta-b)]|\leq$ $\displaystyle|\mathrm{ReLU}(y-b)-\mathrm{ReLU}(y+\delta-b)|+|\mathrm{ReLU}(y-b )-\mathrm{ReLU}(y-\delta-b)|\leq 2\delta$ $$

So that in total the expectation for fixed $a$, over all $b$, is bounded by

$$ $\mathbb{E}_{b\sim[-B,B]}|m_{ReLU}(X;a,b)-m_{ReLU}(X_{\epsilon};a,b)|^{p}\leq(2 \epsilon)^{p}\mu[y-\delta.y+\delta]\leq\frac{1}{2B}(2\epsilon)^{p+1}$ $$

Which implies the same bound when taking the expectation over $a$ and $b$. Thus overall we obtain for all $\beta<(p-1)/p$ that

$$ $\mathbb{E}_{a,b}\left\{\frac{|m_{ReLU}(X;a,b)-m_{ReLU}(X_{\epsilon};a,b)|}{W_{ p}(X,X_{\epsilon})^{\beta}}\right\}^{p}\leq\frac{1}{2B}\frac{(2\epsilon)^{p+1} }{2^{\beta/p}\epsilon^{\beta\cdot p}}=\frac{2^{p+1}}{2B2^{\beta/p}}\epsilon^{p +1-\beta\cdot p}\overset{\epsilon\rightarrow 0}{\rightarrow}0$ $$

While if $m_{ReLU}$ were $\beta$ lower-Hölder in expectation this expression should have been uniformly bounded away from zero.

Part 3: lower-Hölder in expectation
Next we show that $m_{ReLU}$ is $(p+1)/p$ lower-Hölder in expectation. We first consider the restriction of $m_{ReLU}$ to the subspace of $\mathcal{S}_{\leq n}(\Omega)$ which contains only multisets of cardinality exactly $n$. We denote this subspace by $\mathcal{S}_{=n}(\Omega)$. In this case, we realize $m_{ReLU}$ as the composition of two functions:
the functions
$s(X;a)=sort(a^{T}X)$ from Lemma C.4, and the function $q(x;b)=\sum_{i}\mathrm{ReLU}(x_{i}-b)$
Note that indeed $f(X;a,b)=(q\circ s)(a,b)$. Since we already know that $s$ is lower Lipschitz in expectation, it is sufficient to show that $q$ is $(p+1)/p$ lower-Hölder in expectation, due to the theorem on composition B.3.

Let us denote $\hat{B}=\max\{\|x\||\quad x\in\Omega\}$. By assumption $\hat{B}<B$.

Note that the domain of $q$ is contained in

$$ $\Omega_{[-B,B]}=\{x\in\mathbb{R}^{n}|-B\leq x_{1}\leq x_{2}\leq\ldots\leq x_{n }\leq B\}.$ $$

Since the $2$-norm and $\infty$-norm are equivalent on $\mathbb{R}^{n}$ C.1, to address the case of balanced multisets it is sufficient to prove

###### Lemma C.5 .

Let $p>0,B>0$ and $n$ be a natural number. Let $q_{b}(x)=q(x;b)$ be the function described previously, defined on the domain

$$ $\Omega_{[A,B]}={x\in\mathbb{R}^{n}|\quad A\leq x_{1}\leq\ldots\leq x_{n}\leq B}$ $$

There is a constant

$$ $C_{n,p}=\frac{1}{8n(B4-A)^{p}}$ $$

such that for all $x,y\in\Omega_{[A,B]}$,

$$ $\mathbb{E}_{b\sim[A,B]}|q_{b}(x)-q_{b}(y)|^{p}\geq C_{n,p}\|x-y\|_{\infty}^{p+ 1}.$ $$

###### Proof.

Let $x\neq y$ be a pair in $\Omega_{[A,B]}$. Let $s$ be an index for which $|y_{s}-x_{s}|=\|y-x\|_{\infty}$.
without loss of generality assume that $y_{s}>x_{s}$. For every $b\in\mathbb{R}$ we denote $\Delta(b)=q_{b}(y)-q_{b}(x)$. Let $t$ be the smallest integer such that $x_{t}\geq y_{s}$. Note that $t>s$. We now have

$$ $\displaystyle\Delta(y_{s})$ $\displaystyle=\sum_{j=1}^{n}\left[\mathrm{ReLU}(y_{j}-y_{s})-\mathrm{ReLU}(x_{ j}-y_{s})\right]=\sum_{j>s}(y_{j}-y_{s})-\sum_{j\geq t}(x_{j}-y_{s})$ $\displaystyle=\sum_{s<k<t}(y_{k}-y_{s})+\sum_{j\geq t}(y_{j}-x_{j})$ $$

Now

$$ $\displaystyle\Delta(x_{s})$ $\displaystyle=\sum_{i<s}\mathrm{ReLU}(y_{i}-x_{s})+(y_{s}-x_{s})+\sum_{j>s} \left[(y_{j}-x_{s})-(x_{j}-x_{s})\right]$ $\displaystyle\geq(y_{s}-x_{s})+\sum_{j>s}(y_{j}-x_{j})$ $\displaystyle=(y_{s}-x_{s})+\sum_{s<k<t}(y_{k}-y_{s})+\sum_{s<k<t}(y_{s}-x_{k} )+\sum_{j\geq t}(y_{j}-x_{j})$ $\displaystyle\geq(y_{s}-x_{s})+\sum_{s<k<t}(y_{k}-y_{s})+\sum_{j\geq t}(y_{j}- x_{j})$ $\displaystyle=(y_{s}-x_{s})+\Delta(y_{s})$ $$

We deduce that $\Delta(x_{s})-\Delta(y_{s})\geq y_{s}-x_{s}$, and therefore at least one of $|\Delta(x_{s})|$ and $|\Delta(y_{s})|$ is larger than $\frac{y_{s}-x_{s}}{2}$, so we found some $b_{0}$ for which $|\Delta(b_{0})|\geq\frac{y_{s}-x_{s}}{2}$. Next, we note that since $\Delta$ is a sum of $2n$ ReLU functions which are all $1$-Lipschitz, $\Delta$ is $(2n)$ Lipschitz. Therefore if $b$ is such that $|b-b_{0}|\leq\delta:=\frac{1}{8n}(y_{s}-x_{s})$ then

$$ $||\Delta(b)|-|\Delta(b_{0})||\leq|\Delta(b)-\Delta(b_{0})|\leq 2n\delta=\frac{ y_{s}-x_{s}}{4}$ $$

implying that $|\Delta(b)|\geq\frac{y_{s}-x_{s}}{4}$. Thus

$$ $\displaystyle\mathbb{E}_{b}|q_{b}(y)-q_{b}(x)|^{p}$ $\displaystyle=\mathbb{E}_{b}|\Delta(b)|^{p}\geq\mu\{b||b-b_{0}|<\delta\}\cdot \left[\frac{y_{s}-x_{s}}{4}\right]^{p}\geq\frac{\delta}{(B-A)4^{p}}|y_{s}-x_{s }|^{p}$ $\displaystyle=\frac{1}{8n(B-A)4^{p}}|y_{s}-x_{s}|^{p+1}=\frac{1}{8n(B-A)4^{p}} \|y-x\|_{\infty}^{p+1}$ $$

∎

Now let us consider the case of unbalanced multisets $Y,Y^{\prime}$ in $\mathcal{S}_{\leq n}(\Omega)$. Our goal will be to show that the distance between all unbalanced multisets is uniformly bounded from below away from zero.

For fixed $a\in S^{d-1}$ denote $y=a\cdot Y$ and $y^{\prime}=a\cdot Y^{\prime}$. Denote $x=\rho_{(-B)}(y)$, that is, $x$ is the multiset obtained from adding elements with value $-B$ to $y$ until it has $Y$ elements. Similarly, denote $x^{\prime}=\rho_{(-B)}(y^{\prime})$. Note that since $y$ and $y^{\prime}$ don’t have the same number of elements, and all entries of $y$ are in $[-B^{\prime},B^{\prime}]$, we have that

$$ $\|x-x^{\prime}\|_{\infty}\geq B-B^{\prime}.$ $$

For all $b\in[-B,B]$ we have that $\mathrm{ReLU}(-B-b)=0$, and therefore, according to Lemma C.5, we have

$$ $\displaystyle\mathbb{E}_{b}|m_{ReLU}(Y;a,b)$ $\displaystyle-m_{ReLU}(Y^{\prime};a,b)|^{p}=\mathbb{E}_{b}|q_{b}(y)-q_{b}(y^{ \prime})|^{p}$ $\displaystyle=\mathbb{E}_{b}|q_{b}(x)-q_{b}(x^{\prime})|^{p}\geq C_{n,p}\|x-x^ {\prime}\|_{\infty}^{p+1}\geq C_{n,p}(B-B^{\prime})^{p+1}$ $$

We deduce that

$$ $\displaystyle\mathbb{E}_{a,b}|m_{ReLU}(Y;a,b)$ $\displaystyle-m_{ReLU}(Y^{\prime};a,b)|^{p}\geq C_{n,p}^{p}(B-B^{\prime})^{p}$ $\displaystyle=C_{n,p}^{p}(B-B^{\prime})^{p}\left(\frac{W_{1}^{z}(Y,Y^{\prime}) }{W_{1}^{z}(Y,Y^{\prime})}\right)^{p+1}\geq C\left(W_{1}^{z}(Y,Y^{\prime}) \right)^{p+1}$ $$

For an appropriate constant $C$, where we use the fact that $W_{1}^{z}$ is bounded from above. We have obtained a $(p+1)/p$ lower Hölder bound for both balanced and unbalanced multisets, and so we are done.
∎

###### Proof.

###### Lemma C.5 .

###### Proof.

##### Adaptive ReLU

*

###### Proof.

To prove this theorem, we first recall the definition of $m_{\mathrm{ReLU}}^{\mathrm{adapt}}$

$$ $\displaystyle m=\min\{a\cdot x_{1},\ldots,a\cdot x_{r}\},\quad M=\max\{a\cdot x _{1},\ldots,a\cdot x_{r}\},\quad b=(1-t)m+tM$ $\displaystyle m_{\mathrm{ReLU}}^{\mathrm{adapt}}(X;a,t)=[r,m,M,\frac{1}{r}\sum _{i=1}^{r}\mathrm{ReLU}(a\cdot x_{i}-b)]$ $$

To begin with, we note that the case of unbalanced multisets is easy to deal with. There exists constants $0<c=1<C$ such that, for every pair of unbalanced multisets $Y,Y^{\prime}$, and every choice of parameters $a,t$,

$$ $1\leq\|m_{\mathrm{ReLU}}^{\mathrm{adapt}}(X;a,t)-m_{\mathrm{ReLU}}^{\mathrm{ adapt}}(X;a,t)\|_{p}\leq C$ $$

The lower bound follows from the fact that the first coordinate of $m_{\mathrm{ReLU}}^{\mathrm{adapt}}$ is the cardinality of the sets. The upper bound follows from compactness. Similarly, the augmented Wasserstein distance between all unbalanced multisets in $\mathcal{S}_{\leq n}(\Omega)$ is uniformly bounded from above and below. This can be used to obtain both uniform upper Lispchitz bounds, and lower Hölder bounds, as discussed in previous proofs.

Thus, it is sufficient to prove uniform upper Lispchitz bounds, and lower Hölder bounds in expectation, for balanced multisets. Without loss of generality we can assume the balanced multisets both have maximal cardinality $n$. So we need to prove the claim on the space $\mathcal{S}_{=n}(\Omega)$.

We can write $m_{\mathrm{ReLU}}^{\mathrm{adapt}}$, restricted to $\mathcal{S}_{=n}(\Omega)$, as a composition

$$ $m_{\mathrm{ReLU}}^{\mathrm{adapt}}(X;a,t)=Q\circ s(X;a.t)$ $$

where $s(X;a)=\mathrm{sort}(a^{T}X)$, which is the uniformly upper Lipschitz, and lower Lipschitz in expectation function defined in Lemma C.4, and
$Q(x;t)$ is defined via

$$ $m=\mathrm{min}\{x_{1},\ldots,x_{n}\},\quad M=\mathrm{max}\{x_{1},\ldots,x_{n} \},\quad b=(1-t)m+tM$ $$

and

$$ $Q(x;t)=[n,m,M,\frac{1}{n}\sum_{i=1}^{n}\mathrm{ReLU}(x_{i}-b)]$ $$

As $s$ is is uniformly upper Lipschitz, and lower Lipschitz in expectation, it is sufficient to show that $Q$ is upper Lipchitz, and lower Hölder in expectation.

We note that since $\Omega$ is bounded, there exists some $B>0$ such that $\|x\|\leq B,\forall x\in\Omega$, and for this $B$ we have that the image of $s$ is contained in

$$ $\Omega_{[-B,B]}=\{x\in\mathbb{R}^{n}|\quad-B\leq x_{1}\leq\ldots\leq x_{n}\leq B\}$ $$

We can therefore think of $\Omega_{[-B,B]}$ as the domain of $Q$.

We begin with the upper Lipschitz bound. Let $x,y$ be vectors in $\Omega_{[-B,B]}$, which we identify with multisets with $n$ elements in $\mathbb{R}$. Denote the maximum and minimum of $x$ by $M_{x}$ and $m_{x}$. Define the maximum and minimum of $y$ by $M_{y}$ and $m_{y}$. Denote

$$ $\phi(s,t,m,M)=\mathrm{ReLU}(s-[(1-t)m+tM])$ $$

Then for all
$t\in[0,1]$.

$$ $\displaystyle||Q(X;t)-Q(Y;t)||_{p}\leq C\cdot||Q(X;t)-Q(Y;t)||_{1}$ $\displaystyle=C\cdot(|m_{x}-m_{y}|+|M_{x}-M_{y}|+\frac{1}{n}|\sum_{i=1}^{n} \phi(x_{i},t,m_{x},M_{x})-\sum_{i=1}^{n}\phi(y_{i},t,m_{y},M_{y})|)$ $\displaystyle\leq C\cdot(|m_{x}-m_{y}|+|M_{x}-M_{y}|+\frac{1}{n}\sum_{i=1}^{n} |x_{i}-y_{i}+t(M_{x}-M_{x})+(1-t)(m_{x}-m_{y})|)$ $\displaystyle\leq C\cdot(|m_{x}-m_{y}|+|M_{x}-M_{y}|+\frac{1}{n}\sum_{i=1}^{n} (|x_{i}-y_{i}|+|(M_{x}-M_{y})|+|(m_{x}-m_{y}))|)$ $\displaystyle\leq C\cdot(||x-y||_{\infty}+||x-y||_{\infty}+\frac{1}{n}\sum_{i= 1}^{n}3\cdot||x-y||_{\infty})$ $\displaystyle=5C\cdot||x-y||_{\infty}\leq 5CC^{\prime}\cdot||x-y||_{p}$ $$

Where $C,C^{\prime}$ are the constants obtained from norm equivalence over $\mathbb{R}^{n}$.

To obtain a Hölder lower bound, the idea of the proof is that for given balanced multisets $x,y\in\mathbb{R}^{n}$, we know from the analysis of $m_{\mathrm{ReLU}}$ that we can get a lower-Hölder bound when considering biases going between

$$ $m_{x,y}=\min\{m_{x},m_{y}\}\text{ and }M_{x,y}=\max\{M_{x},M_{y}\}$ $$

and then showing that the difference between this case and the function $Q$ where the bias range depends on the maximum and minimum of the individual multisets $x,y$, is proportional to the different between the minimum and maximum of $x$ and $y$, which also appear in $Q$. Indeed, since $\mathrm{ReLU}$ is Lipschitz we have for every $s,m,M,\hat{m},\hat{M}$ in $\mathbb{R}$ and $t\in[0,1]$ that

$$ $\displaystyle|\phi(s,t,m,M)-\phi(s,t,\hat{m},\hat{M})|$ $\displaystyle\leq|-t(M-\hat{M})-(1-t)(m-\hat{m})|$ (6) $\displaystyle\leq t|M-\hat{M}|+(1-t)|m-\hat{m}|\leq|M-\hat{M}|+|m-\hat{m}|$ $$

Next, to bound $\|Q(x;t)-Q(y;t)\|_{p}$, due to equivalence of norms, it is sufficient to bound $\|Q(x;t)-Q(y;t)\|_{1}$. We obtain

$$ $\displaystyle\|Q(x;t)-Q(y;t)\|_{1}=$ $\displaystyle|m_{x}-m_{y}|+|M_{x}-M_{y}|$ $\displaystyle+\frac{1}{n}|\sum_{i=1}^{n}\phi(x_{i},t,m_{x},M_{x})-\sum_{i=1}^{ n}\phi(y_{i},t,m_{y},M_{y})|$ $\displaystyle=|m_{x}-m_{xy}|+|M_{x}-M_{xy}|+|m_{xy}-m_{y}|+|M_{xy}-M_{y}|$ $\displaystyle+\frac{1}{n}|\sum_{i=1}^{n}\phi(x_{i},t,m_{x},M_{x})-\sum_{i=1}^{ n}\phi(y_{i},t,m_{y},M_{y})|$ $\displaystyle\stackrel{{\scriptstyle\eqref{eq:helper}}}{{\geq}}\frac{1}{n}| \sum_{i=1}^{n}\phi(x_{i},t,m_{xy},M_{xy})-\sum_{i=1}^{n}\phi(x_{i},t,m_{x},M_{ x})|$ $\displaystyle+\frac{1}{n}|\sum_{i=1}^{n}\phi(y_{i},t,m_{y},M_{y})-\sum_{i=1}^{ n}\phi(y_{i},t,m_{xy},M_{xy})|$ $\displaystyle+\frac{1}{n}|\sum_{i=1}^{n}\phi(x_{i},t,m_{x},M_{x})-\sum_{i=1}^{ n}\phi(y_{i},t,m_{y},M_{y})|$ $\displaystyle\stackrel{{\scriptstyle\text{triangle ineq.}}}{{\geq}}\frac{1}{n} |\sum_{i=1}^{n}\phi(x_{i},t,m_{xy},M_{xy})-\sum_{i=1}^{n}\phi(x_{i},t,m_{x},M_ {x})$ $\displaystyle+\sum_{i=1}^{n}\phi(y_{i},t,m_{y},M_{y})-\sum_{i=1}^{n}\phi(y_{i} ,t,m_{xy},M_{xy})$ $\displaystyle+\sum_{i=1}^{n}\phi(x_{i},t,m_{x},M_{x})-\sum_{i=1}^{n}\phi(y_{i} ,t,m_{y},M_{y})|$ $\displaystyle=\frac{1}{n}|\sum_{i=1}^{n}\phi(x_{i},t,m_{xy},M_{xy})-\sum_{i=1} ^{n}\phi(y_{i},t,m_{xy},M_{xy})|$ $\displaystyle=|q(x;(1-t)m_{xy}+tM_{xy})-q(y;(1-t)m_{xy}+tM_{xy})|$ $$

By Lemma C.5, we deduce that

$$ $\displaystyle\mathbb{E}_{t\sim[0,1]}\|Q(x;t)-Q(y;t)\|_{1}^{p}$ $\displaystyle\geq\mathbb{E}_{b\sim[m_{xy},M_{xy}]}|q(x;b)-q(y;b)|^{p}$ $\displaystyle\geq\frac{1}{8n4^{p}(M_{xy}-m_{xy})}\|y-x\|_{\infty}^{p+1}\geq \frac{1}{16nB4^{p}}\|y-x\|_{\infty}^{p+1}$ $$

by invoking the equivalence of the infinity norm and $p$ norm we are done.
∎

###### Proof.

## Appendix D Lipschitz COMBINE operations

In this section we describe how to construct COMBINE functions which are both uniformly Lipschitz, and lower-Lipschitz in expectation. The input to these functions are a pair of vectors $x_{1},y_{1}$ with the metric $\|x_{1}\|_{p}+\|y_{1}\|_{p}$. The output will be a vector (or scalar), and the metric on the output space will again be the $p$ norm.

Our analysis covers four 2-tuple embeddings, all of which are uniformly upper Lipschitz, and lower Lipschitz in expectation.

###### Theorem D.1 .

Let $p>0$, the following 2-tuple embeddings are all bi-Lipschitz in expectation with respect to $l_{p}$ and the 2-tuple metric previously defined.

- 1.
Linear combination: Given a 2-tuple $(x,y)\in\mathbb{R}^{d}\times\mathbb{R}^{d}$ and $\alpha\sim U[-D,D]$ we define
$f(x,y;\alpha)=\alpha\cdot x+y$
- 2.
Linear transform and sum: Given a 2-tuple $(x,y)\in\mathbb{R}^{k}\times\mathbb{R}^{l}$, and $A\in\mathbb{R}^{l\times k}$, whose rows are independently and uniformly sampled from $S^{k-1}$, we define
$f(x,y;A)=Ax+y$
- 3.
Concatenation:
This embedding concatenates the tuple entries, and has no parameters
$f(x,y)=\begin{bmatrix}x\\
y\end{bmatrix}$
- 4.
Concat and project: Given $(x,y)\in\mathbb{R}^{k}\times\mathbb{R}^{l}$ and $\theta\sim S^{k+l-1}$
$f(x,y;\theta)=\theta\cdot\begin{bmatrix}x\\
y\end{bmatrix}$

We note that the last method, concat and project, is the method used in the definition of SortMPNN in the main text. We also note that the first method, linear combination, is the method used in (Xu et al., ). Unlike the other methods, it requires the vectors $x,y$ to be of the same dimension.

In order to compare the proposed COMBINE functions which all share the bi-Lipschitz in expectation property, one could further explore separation quality by comparing the distortion defined by $\frac{M}{m}$ where $M$ ($m$) is the upper (lower) Lipschitz in expectation bound. A lower distortion would indicate that the function doesn’t change the input metric as much (possibly up to multiplication by some constant), resulting in high quality separation.

In figure 5, we plot the empirical distortion of the different COMBINE functions with varying input and embedding dimensions. The embedding dimension is controlled by stacking multiple instances of the functions with independent parameters. Note that concatenation isn’t parametric and therefore we only have a single output dimension for it.

The experiment is run on $1,000$ different random tuple pairs $\{(x,y),(z,w)\}$ for each of the four input dimensions experimented with, where $x,y,z,w\in\mathbb{R}^{\text{in_dim}}$. All random data vectors are sampled with entries drawn from the normal distribution. The empirical distortion is then computed by taking the largest empirical ratio $M$ between the input and output distances and dividing it by the smallest ratio $m$ between the input and output distances.

As we see in the figure, all proposed functions stabilize around a constant value as the embedding dimension increases. This value is the expected distortion on the experiment data. We can see that concat and project maintains low expected distortion across all settings, while linear combination has a higher expected distortion but not by much. Linear transform and sum seems to have higher expected distortion for lower input dimensions, but it improves with large input dimensions, matching concat and project distortion.

Another interesting result is that we can see that concat and project seems to have the highest variance of the three for most input dimensions, since the empirical distortion for low embedding dimensions is higher than other functions in most cases.

Figure: Figure 5: Distortion of the 2-tuple embeddings as a function of input and embedding dimension. LC=linear combination, LTSum=linear transform and sum, CP=concat project, C=concatenation.
Refer to caption: extracted/5914493/figs/combine_experiment.png

###### Proof of Theorem D.1 .

Linear combination: Let $x,z\in\mathbb{R}^{k},y,w\in\mathbb{R}^{l}$
We begin by proving that $f$ is lower Lipschitz in expectation.

$$ $\mathbb{E}_{\alpha\sim U[-D,D]}[||f(x,y;\alpha)-f(z,w;\alpha)||_{p}^{p}]= \mathbb{E}_{\alpha\sim U[-D,D]}[||(\alpha\cdot x+y)-(\alpha\cdot w+z)||_{p}^{p}]$ $$

$$ $=\mathbb{E}_{\alpha\sim U[-D,D]}[||\alpha\cdot(x-z)+(y-w)||_{p}^{p}]$ $$

From norm equivalence C.1, $\exists c_{1}$ s.t.

$$ $\geq c_{1}\cdot\mathbb{E}_{\alpha\sim U[-D,D]}[||\alpha\cdot(x-z)+(y-w)||_{1}^ {p}]$ $$

$$ $\stackrel{{\scriptstyle\text{Jensen's ineq.}}}{{\geq}}c_{1}\cdot\mathbb{E}_{ \alpha\sim U[-D,D]}[||\alpha\cdot(x-z)+(y-w)||_{1}]^{p}$ $$

$$ $=c_{1}\cdot{E}_{\alpha\sim U[-D,D]}[||\alpha\cdot(x-z)||_{1}+||(y-w)||_{1}]^{p}$ $$

From norm equivalence C.1, $\exists c_{2}$ s.t.

$$ $\geq c_{1}\cdot c_{2}\cdot({E}_{\alpha\sim U[-D,D]}[|\alpha|]\cdot||(x-z)||_{p }+||(y-w)||_{p})^{p}$ $$

$$ $=c_{1}\cdot c_{2}\cdot(\frac{D}{2}\cdot||(x-z)||_{p}+||(y-w)||_{p})^{p}$ $$

$$ $\geq c_{1}\cdot c_{2}\cdot min(1,\frac{D}{2}^{p})\cdot(||(x-z)||_{p}+||(y-w)|| _{p})^{p}$ $$

For uniform upper Lipschitz bounds we have

$$ $\displaystyle||f(x,y;\alpha)-f(z,w;\alpha)||_{p}$ $\displaystyle=||(\alpha\cdot x+y)-(\alpha\cdot w+z)||_{p}$ $\displaystyle\leq\|\alpha\cdot x-\alpha\cdot w\|_{p}+\|y-z\|_{p}$ $\displaystyle\leq\mathrm{max}(1,D)^{p}\cdot(||(x-z)||_{p}+||(y-w)||_{p})$ $$

Linear transform and sum: We prove lower-Lipschitzness in expectation. Let $a_{i}$ denote the i’th row of $A$. Recall that for any $1\leq i\leq l$, $a_{i}$ is drawn uniformly from $S^{k-1}$.
Let $x,z\in\mathbb{R}^{k},y,w\in\mathbb{R}^{l}$

$$ $\mathbb{E}_{A}[||f(x,y;A)-f(z,w;A)||_{p}^{p}]=\mathbb{E}_{A}[||A(x-z)+(y-w)||_ {p}^{p}]$ $$

$$ $=\mathbb{E}_{A}[\sum_{i=1}^{l}|a_{i}\cdot(x-z)|^{p}+\sum_{j=1}^{l}|y_{j}-w_{j} |^{p}]=\sum_{i=1}^{l}\mathbb{E}_{a_{i}\sim S^{k-1}}[|a_{i}\cdot(x-z)|^{p}]+ \sum_{j=1}^{l}|y_{j}-w_{j}|^{p}$ $$

From lemma C.3 $\exists b>0$ s.t.

$$ $\displaystyle\mathbb{E}_{A}[||f(x,y;A)-f(z,w;A)||_{p}^{p}]$ $\displaystyle\geq n\cdot b\cdot||x-z||_{p}^{p}+||y-w||_{p}^{p}$ $$

Let us denote $t=\begin{bmatrix}||x-z||_{p}\\
||y-w||_{p}\end{bmatrix}\in\mathbb{R}^{2}$.
Then

$$ $n\cdot b\cdot||x-z||_{p}^{p}+||y-w||_{p}^{p}\geq min(1,n\cdot b)\cdot(||x-z||_ {p}^{p}+||y-w||_{p}^{p})=min(1,n\cdot b)\cdot||t||_{p}^{p}$ $$

From norm equivalence C.1, $\exists c>0$

$$ $\geq min(1,n\cdot b)\cdot c\cdot||t||_{1}^{p}=min(1,n\cdot b)\cdot c\cdot(||x- z||_{p}+||y-w||_{p})^{p}$ $$

Uniform upper Lipschitzness is straightforward to prove.

Concatenation: We show that this non-parametric function is bi-Lipschitz.
Let $x,z\in\mathbb{R}^{k},y,w\in\mathbb{R}^{l}$. Denote as before
$t=\begin{bmatrix}||x-z||_{p}\\
||y-w||_{p}\end{bmatrix}\in\mathbb{R}^{2}$. Then
the distance difference in the output space is given by $||f(x,y)-f(z,w)||_{p}^{p}=\|t\|_{p}$, while the difference in the input space is given by

$$ $\|x-z\|_{p}+\|y-w\|_{p}=\|t\|_{1}.$ $$

The claim follows from equivalence of $p$ norm and $1$ norm on $\mathbb{R}^{2}$ C.1.

Concat and project: The claim follows from the fact that the concatenation operation is bi-Lipschitz, and Lemma C.3.
∎

###### Theorem D.1 .

###### Proof of Theorem D.1 .

## Appendix E Tree Mover’s Distance

In this appendix section we review the definition of the Tree Mover’s Distance (TMD) from (Chuang & Jegelka, ), which is the way we measures distances between graphs in this paper.

We first review Wasserstein distances. Recall that if $(X,D)$ is a metric space, $\Omega\subseteq X$ is a subset, and $z$ is some point in $X$ with $dist(z,\Omega)>0$, then we can define the Wasserstein distance on the space of multisets consisting of $n$ elements in $\Omega$ via

$$ ${W}_{1}(\{\!\!\{x_{1},\ldots,x_{n}\}\!\!\},\{\!\!\{y_{1},\ldots,y_{n}\}\!\!\}) =\min_{\tau\in S_{n}}\sum_{j=1}^{n}D(x_{j},y_{\tau(j)})$ $$

The augmentation map on multisets of size $r\leq n$ is defined as

$$ $\rho_{(z)}\left(\{\!\!\{x_{1},\ldots,x_{r}\}\!\!\}\right)=\{\!\!\{x_{1},\ldots ,x_{r},x_{r+1}=z,\ldots.x_{n}=z\}\!\!\}$ $$

and the augmented distance on multisets of size up to $n$ is defined via

$$ $W_{D,1}^{(z)}(X,\hat{X})=W_{D,1}(\rho_{(z)}X,\rho_{(z)}\hat{X})$ $$

We now return to define the TMD. We consider the space of graphs $G_{\leq n}(\Omega)$, consisting of graphs with $\leq n$ nodes, with node features coming from a compact domain $\Omega\subseteq\mathbb{R}^{d}$. We also fix some $z\in\mathbb{R}^{d}\setminus\Omega$. The TMD is defined using the notion of computation trees:

###### Definition E.1 .

(Computation Trees). Given a graph $G=(V,E)$ with node features $\{\!\!\{x_{v}\}\!\!\}_{v\in V}$, let $T_{v}^{(0)}$ be the rooted tree with a single node $v$, which is also the root of the tree, and node features $x_{v}$. For $K\in\mathbb{N}^{+}$ let $T_{v}^{(K)}$ be the depth-$K$ computation tree of node $v$ constructed by connecting the neighbors of the leaf nodes of $T_{v}^{(K-1)}$ to the tree. Each node is assigned the same node feature it had in the original graph $G$. The multiset of depth-$K$ computation trees defined by $G$ is denoted by $\mathcal{T}_{G}^{(K)}:=\{\!\!\{T_{v}^{(K)}\}\!\!\}_{v\in V}$.
Additionally, for a tree $T_{r}$ with root $r$, we denote by $\mathcal{T}_{r}$ the multiset of subtrees that root at the descendants of $r$.

###### Definition E.2 .

(Blank Tree). A blank tree $\bar{T}_{z}$ is a tree (graph) that contains a single node and no edge, where the node feature is the blank vector $z$.

Recall that by assumption, all node features will come from the compact set $\Omega$, and $z\not\in\Omega$.

We can now define the tree distance:

###### Definition E.3 .

(Tree Distance).(^4^44Note the difference from the original definition in (Chuang & Jegelka, ) is due to our choice to set the depth weight to 1 and using the 1-Wasserstein which is equivalent to optimal transport)
The distance between two trees $T_{a},T_{b}$ with features from $\Omega$ and $z\not\in\Omega$, is defined recursively as

$$ $TD(T_{a},T_{b}):=\begin{cases}\|x_{a}-x_{b}\|_{p}+W_{TD,1}^{(\bar{T}_{z})}( \mathcal{T}_{a},\mathcal{T}_{b})&\text{if $K>0$}\\ \|x_{a}-x_{b}\|_{p}&\text{otherwise}\end{cases}$ $$

where $K$ denotes the maximal depth of the trees $T_{a}$ and $T_{b}$.

###### Definition E.4 .

(Tree Mover’s Distance). Given two graphs, $G,H$ and $w,K\geq 0$, the tree mover’s distance is defined as

$$ $TMD^{(K)}(G,H)=W^{(\bar{T}_{z})}_{TD,1}(\mathcal{T}_{G}^{(K)},\mathcal{T}_{H}^ {(K)})$ $$

where $\mathcal{T}_{G}^{(K)}$ and $\mathcal{T}_{H}^{(K)}$ denote the multiset of all depth $K$ computational trees arising from the graphs $G$ and $H$, respectively.
We refer the reader to (Chuang & Jegelka, ) where they prove this is a pseudo-metric that fails to distinguish only graphs which cannot be separated by $K$ iterations of the WL test.

###### Definition E.1 .

###### Definition E.2 .

###### Definition E.3 .

###### Definition E.4 .

## Appendix F MPNN Hölder Proofs

In the main text we informally stated the following theorem
\LipschitzMPNN*

Before stating the formal theorem and proof, we will present an important lemma alongside some assumptions and notation that will be used throughout this section.

###### Lemma F.1 .

Let $\Omega\subseteq\mathbb{R}^{d}$ be a compact set and let $f:G_{\leq n}(\Omega)\rightarrow\mathbb{R}^{m}$ be a depth $K$ MPNN such that all the aggregation $\phi^{(k)}$, combine $\psi^{(k)}$ and readout $\eta$ are continuous. Let $z_{k}\not\in\Omega_{k}$, where $\Omega_{k}$ is the set of all possible feature vectors after running $k\leq K$ message passing layers. Then, for any $G\in G_{\leq n}(\Omega),v\in V_{G}$ and for any $p\geq 1$, there exist $C\geq c>0$ such that

$$ $C\cdot TD(T_{v}^{(k)},\bar{T}_{z})\geq||z_{k}-x_{v}^{(k)}||_{p}\geq c\cdot TD( T_{v}^{(k)},\bar{T}_{z})$ $$

###### Proof.

We begin by showing that $\Omega_{k}$ is compact.

Since the number of nodes in graphs from $G_{\leq n}(\Omega)$ is bounded, the set of all possible graph topologies, i.e. pairs $(V,E)$ with at most $n$ nodes, is finite. Let $(V,E)$ be such a fixed topology. For every $k$, the function $f^{(k)}_{V,E}\left((x_{v})_{v\in\Omega},\theta\right)$ mapping features $\Omega$ and parameters $\theta$ to new parameter $x_{v}^{(k)}$ is continuous since the aggregation and combine functions are continuous. Accordingly the image $f{(k)}_{V,E}(\Omega)$ is compact, and therefore

$$ $\Omega_{k}=\bigcup_{(V,E),|V|\leq n}f^{(k)}_{V,E}(\Omega)$ $$

is compact as well.

Therefore:

$$ $a:=\inf_{x\in\Omega_{k}}||x-z_{k}||_{p}>0$ $$

$$ $d:=\sup_{x\in\Omega_{k}}||x-z_{k}||_{p}<\infty$ $$

Additionally, since $\Omega$ is compact and $z\not\in\Omega$, we can show recursively, starting from $k=1$, that for all $k\in\{0\ldots,K\}$ and all trees $T_{v}^{(k)}\in\mathcal{T}^{(k)}_{G_{\leq n}(\Omega)}$,

$$ $b:=\inf_{T_{v}^{(k)}\in\mathcal{T}^{(k)}_{G_{\leq n}(\Omega)}}TD(\bar{T}_{z},T _{v}^{(k)})>0$ $$

$$ $e:=\sup_{T_{v}^{(k)}\in\mathcal{T}^{(k)}_{G_{\leq n}(\Omega)}}TD(\bar{T}_{z},T _{v}^{(k)})<\infty$ $$

Where $\mathcal{T}^{(k)}_{G_{\leq n}(\Omega)}$ denotes the set of all possible height-$k$ computation trees from $G_{\leq n}(\Omega)$.
Therefore, $\forall G\in G_{\leq n}(\Omega),v\in V_{G}$:

$$ $\frac{d}{b}\cdot TD(\bar{T}_{z},T_{v}^{k})\geq d\cdot\frac{b}{b}=d\geq||z_{k}- x_{v}^{(k)}||_{p}\geq a=a\cdot\frac{e}{e}\geq\frac{a}{e}\cdot TD(\bar{T}_{z},T _{v}^{k})$ $$

To conclude we define $C:=\frac{d}{b},c:=\frac{a}{e}$.

∎

###### Lemma F.1 .

###### Proof.

##### Assumptions and notations

As stated in the main text, we consider graphs with up to $n$ nodes. We will also make the disjointness assumption: We assume as previously that the initial features of the graphs all reside in a compact set $\Omega=\Omega_{0}\subseteq\mathbb{R}^{d}$. We denote by $\Omega_{k}$ the space of all possible features which can be obtained by the MPNN at question after $k$ iteration, with any choice of parameters. This set is also compact (see proof of F.1). The disjointness assumptions is that the ’augmentation vector’ $z_{k}$ is not in $\Omega_{k}$, for all $k=0,1,\ldots,K$.

We define for $k=0,\ldots,K-1$

$$ $x_{v,s}^{(k)}=\left\{\begin{array}[]{ll}x_{s}^{(k)}&\mbox{if }(v,s)\in E\\ z_{k}&\mbox{if }otherwise\end{array}\right.,\text{ and }T_{v,s}^{(k)}=\left\{ \begin{array}[]{ll}T_{s}^{(k)}&\mbox{if }(v,s)\in E\\ \bar{T}_{z}&\mbox{if }otherwise\end{array}\right.$ $$

For $K$ we denote

$$ $x_{s}^{(K)}=\left\{\begin{array}[]{ll}x_{s}^{(K)}&\mbox{if }s\in V\\ z_{K}&\mbox{if }otherwise\end{array}\right.,\text{ and }T_{s}^{(K)}=\left\{ \begin{array}[]{ll}T_{s}^{(K)}&\mbox{if }s\in V\\ \bar{T}_{z}&\mbox{if }otherwise\end{array}\right.$ $$

We will also denote

$$ $x_{v,\bullet}^{(k)}=[x_{v,1}^{(k)},\ldots,x_{v,n}^{(k)}]\text{ and }x_{\bullet }^{(K)}=[x_{1}^{(K)},\ldots,x_{n}^{(K)}]$ $$

We can now state the full formal statement of the theorem:

###### Theorem F.2 .

(Uniformly Lipschitz MPNN embeddings, full version)
Let $p\geq 1$ and $K\in\mathbb{N}$,
Let $f:G_{\leq n}(\Omega)\rightarrow\mathbb{R}^{m}$ be a continuous MPNN with $K$ message passing layers. Under the above assumptions F and given the following holds for all $1\leq k\leq K$:

- 1.
The aggregation, $\phi^{(k)}$, is uniformly upper Lipschitz w.r.t. the augmented Wasserstein distance $W_{1}^{(z_{k})}$ on $\mathcal{S}_{\leq n}(\Omega_{k})$.
- 2.
The combine function $\psi^{(k)}$, is uniformly upper Lipschitz.
- 3.
the readout function $\eta$, is uniformly upper Lipschitz w.r.t. the augmented Wasserstein distance $W_{1}^{z_{(K+1)}}$ on $\mathcal{S}_{\leq n}(\Omega_{K+1})$..

Then, there exist constants $C_{1},C_{2}>0$ such that

- 1.
The node embeddings after $k$ message passing layers are uniformly upper Lipschitz w.r.t. $TD$ on their depth $k$ computation trees
$\|x^{(k)}_{v}-x^{(k)}_{u}\|_{p}^{p}\leq C_{1}\cdot TD^{p}(T^{(k)}_{v},T^{(k)}_
{u})$
- 2.
The graph embeddings are uniformly upper Lipschitz w.r.t. $TMD^{(K)}$
$\|c^{\mathrm{global}}-\hat{c}^{\mathrm{global}}\|_{p}^{p}\leq C_{2}\cdot TMD^{
(K)}(G,\hat{G})^{p}$

###### Proof.

Since we are proving for uniformly Lipschitz, we will allow ourselves to omit the parameters for ease of notation, as the proof holds for any set of parameters.

###### Theorem F.2 .

###### Proof.

##### Proof of First Claim

Let $u,v\in G$. The proof will be by induction on $k$. for $k=0$, by definition

$$ $\displaystyle\|x_{v}^{(0)}-x_{u}^{(0)}\|_{p}^{p}=TD^{p}(T_{v}^{(0)},T_{u}^{(0)})$ $$

We now assume correctness for $k-1$ and prove for $k$.

- 1.
First, using the properties of $\phi^{(k)}$ and the induction hypothesis we get:
$\displaystyle\|\phi^{(k)}(\mathcal{N}_{v}^{(k-1)})-\phi^{(k)}(\mathcal{N}_{u}^
{(k-1)})\|_{p}^{p}\leq c_{\phi}W_{1}^{p}(x^{(k-1)}_{v,\bullet},x^{(k-1)}_{u,
\bullet})$
$\displaystyle\stackrel{{\scriptstyle\ref{wass_equiv}}}{{\leq}}c_{\phi}\cdot c_
{1}W_{p}^{p}(x^{(k-1)}_{v,\bullet},x^{(k-1)}_{u,\bullet})=c_{\phi}\cdot c_{1}(
\min_{\tau\in S_{n}}\sum_{s=1}^{n}\|x^{(k-1)}_{v,s}-x^{(k-1)}_{u,\tau(s)}\|_{p
})^{p}$
$\displaystyle\stackrel{{\scriptstyle\ref{norm_equiv}}}{{\leq}}c_{\phi}\cdot c_
{1}\cdot c_{2}(\min_{\tau\in S_{n}}\sum_{s=1}^{n}\|x^{(k-1)}_{v,s}-x^{(k-1)}_{
u,\tau(s)}\|_{p}^{p})$
$\displaystyle\stackrel{{\scriptstyle\text{ind. hypothesis}+\ref{lemma:blank_
tree_equiv}}}{{\leq}}c_{\phi}\cdot c_{1}\cdot c_{2}\cdot c_{(k-1)}\sum_{s=1}^{
n}TD^{p}(T^{(k-1)}_{v,s},T^{(k-1)}_{u,\tau^{*}(s)})$
$\displaystyle\stackrel{{\scriptstyle\ref{norm_equiv}}}{{\leq}}c_{\phi}\cdot c_
{1}\cdot c_{2}\cdot c_{(k-1)}\cdot c_{3}(\sum_{s=1}^{n}TD(T^{(k-1)}_{v,s},T^{(
k-1)}_{u,\tau^{*}(s)}))^{p}$
$\displaystyle=cW^{p}_{TD,1}(T^{(k-1)}_{v,\bullet},T^{(k-1)}_{u,\bullet})$
Where $\tau^{*}:=\mathrm{argmin}_{\tau\in S_{n}}\sum_{s=1}^{n}TD(T^{(k-1)}_{v,s},T^{(
k-1)}_{u,\tau(s)})$, $c_{1},c_{2},c_{3}$ are the relevant constants from the lemmas we used and $c$ is the multiplication of all previous constants.
- 2.
Next, we note that
$\displaystyle TD(T_{v}^{(k-1)},T_{u}^{(k-1)})\leq TD(T_{v}^{(k)},T_{u}^{(k)})$
Since the depth $k-1$ computation trees are subtrees of the depth $k$ computation trees.
- 3.
Next, using the properties of $\psi^{(k)}$, the induction hypothesis and the above we get:
$\displaystyle\|x_{v}^{(k)}-x_{u}^{(k)}\|_{p}^{p}=\|\psi^{(k)}(x_{v}^{(k-1)},
\phi^{(k)}(\mathcal{N}_{v}^{(k-1)}))-\psi^{(k)}(x_{u}^{(k-1)},\phi^{(k)}(
\mathcal{N}_{u}^{(k-1)}))\|_{p}^{p}$
$\displaystyle\leq c_{\psi}\cdot(\|x_{v}^{(k-1)}-x_{u}^{(k-1)}\|_{p}^{p}+\|\phi
^{(k)}(\mathcal{N}_{v}^{(k-1)})-\phi^{(k)}(\mathcal{N}_{u}^{(k-1)})\|_{p}^{p})$
$\displaystyle\stackrel{{\scriptstyle\text{ind. hypothesis}}}{{\leq}}c_{\psi}
\cdot(c_{(k-1)}TD^{p}(T^{(k-1)}_{v},T^{(k-1)}_{u})+\|\phi^{(k)}(\mathcal{N}_{v
}^{(k-1)})-\phi^{(k)}(\mathcal{N}_{u}^{(k-1)})\|_{p}^{p})$
$\displaystyle\stackrel{{\scriptstyle\ref{neighborhood_prop}}}{{\leq}}c_{\psi}
\cdot(c_{(k-1)}TD^{p}(T^{(k-1)}_{v},T^{(k-1)}_{u})+cW^{p}_{TD,1}(T^{(k-1)}_{v,
\bullet},T^{(k-1)}_{u,\bullet}))$
$\displaystyle\stackrel{{\scriptstyle\ref{TD_monotonicity}}}{{\leq}}c_{\psi}
\cdot(c_{(k-1)}TD^{p}(T^{(k)}_{v},T^{(k)}_{u})+cW^{p}_{TD,1}(T^{(k-1)}_{v,
\bullet},T^{(k-1)}_{u,\bullet}))$
$\displaystyle\stackrel{{\scriptstyle\ref{norm_equiv}}}{{\leq}}c_{\psi}\cdot max
(c_{(k-1)},c)\cdot c_{4}(TD(T^{(k)}_{v},T^{(k)}_{u})+W_{TD,1}(T^{(k-1)}_{v,
\bullet},T^{(k-1)}_{u,\bullet}))^{p}$
$\displaystyle\leq c_{\psi}\cdot max(c_{(k-1)},c_{1})\cdot c_{4}\cdot 2^{p}
\cdot TD^{p}(T^{(k)}_{v},T^{(k)}_{u})$

Where $c_{4}$ is the relevant constant from norm equivalence,
concluding the proof of the first claim.

##### Proof of Second Claim

Let us denote the final message passing features corresponding to the nodes of G by $x_{v}^{(K)}$, and the depth-$K$ computation trees corresponding to each node of G by $T_{v}^{(K)}$. Similarly, we denote the features corresponding to the nodes of $\hat{G}$ by $\hat{x}_{v}^{(K)}$ and the trees by $\hat{T}_{v}^{(K)}$. Then, using the properties of $\eta$ and first part of the theorem we get

$$ $\displaystyle\|c^{\mathrm{global}}-\hat{c}^{\mathrm{global}}\|_{p}^{p}=\|\eta( x^{(K)}_{\bullet})-\eta(\hat{x}^{(K)}_{\bullet})\|_{p}^{p}$ $\displaystyle\leq c_{\eta}W_{1}^{p}(x^{(K)}_{\bullet},\hat{x}^{(K)}_{\bullet}) \stackrel{{\scriptstyle\ref{wass_equiv}}}{{\leq}}c_{\eta}\cdot c_{1}W_{p}^{p}( x^{(K)}_{\bullet},\hat{x}^{(K)}_{\bullet})=c_{\eta}\cdot c_{1}(\min_{\tau\in S _{n}}\sum_{s=1}^{n}\|x^{(K)}_{s}-\hat{x}^{(K)}_{\tau(s)}\|_{p}^{p})$ $\displaystyle\stackrel{{\scriptstyle(*)}}{{\leq}}c_{\eta}\cdot c_{1}\cdot C \cdot\sum_{s=1}^{n}TD^{p}(T^{(K)}_{s},\hat{T}^{(K)}_{\tau^{*}(s)})\stackrel{{ \scriptstyle\ref{norm_equiv}}}{{\leq}}c_{\eta}\cdot c_{1}\cdot C\cdot c_{2} \cdot(\sum_{s=1}^{n}TD(T^{(K)}_{s},\hat{T}^{(K)}_{\tau^{*}(s)}))^{p}$ $\displaystyle=c_{\eta}\cdot c_{1}\cdot C\cdot c_{2}\cdot TMD^{(K)}(G,\hat{G})^ {p}$ $$

Where (*) is the first claim and lemma F.1, $c_{1},c_{2}$ are the relevant constants from norm equivalence and $\tau^{*}=\mathrm{argmin}_{\tau\in S_{n}}\sum_{s=1}^{n}TD(T^{(K)}_{s},\hat{T}^{
(K)}_{\tau(s)})$, concluding the proof.

∎

### F.1 SortMPNN

In the main text we stated the informal theorem
\ThmSortMPNNInformal*

We next state the theorem formally and present the proof.

Our results in this section hold for any number of MPNN iterations $K$, and any choice of ’widths’ $W_{1},\ldots,W_{K+1}$. For simplicity we prove this for the apriori hardest case, $W_{1}=1=\ldots=W_{K+1}$. We will call this version of SortMPNN thin-SortMPNN.

Explicitly, thin-SortMPNN is defined by:

SortMPNN: For $\mathbf{k=1\ldots,K}$

$$ $\displaystyle\text{{AGGREGATE}: }c_{v}^{(k)}=\mathrm{sort}\left(a^{(k)}\cdot x _{v,\bullet}^{(k-1)}\right)$ $\displaystyle\text{{COMBINE}: }x_{v}^{(k)}=d^{(k)}\cdot concat(x_{v}^{(k-1)},c _{v}^{(k)})$ $$

Note that we omit the second inner product in the definition of $S_{z}$ as it is superfluous: this inner-product would be subsumed by the following inner product with $d^{(k)}$ in the COMBINE function.

The READOUT function is given by

$$ $\textit{READOUT}:c^{\mathrm{global}}=b^{({K+1})}\cdot\mathrm{sort}\left(a^{(K+ 1)}\cdot x_{\bullet}^{(K)}\right),$ $$

We will denote by $\theta^{(k)}$ the concatenation of all network parameters up to the creation of the node features $x_{v}^{(k)}$, where $\theta^{(0)}$ is just an ’empty vector’. We denote all network parameters, including those used by the readout function by $\theta^{(K+1)}$. The distribution on each of the parameter vectors $a^{(k)},b^{(k)},d^{(k)}$ is taken to be uniform on the unit sphere of the relevant dimension, as discussed in the main text.

We will now state our theorem on lower Lipschitzness of thin-SortMPNN.

###### Theorem F.3 .

(lower Lipschitz SortMPNN, formal)
Under the above assumptions F, and given $a^{(k)},b^{(k)}$ are distributed uniformly on the appropriate unit sphere, then for thin-SortMPNN with depth $K$ the following holds:

- 1.
The node embeddings after $k$ message passing layers are lower Lipschitz in expectation w.r.t. $TD$ on their depth $k$ computation trees
$\mathbb{E}_{\phi^{(k)}}[\|x^{(k)}_{v}-x^{(k)}_{u}\|_{p}^{p}]\geq c_{1}\cdot TD
^{p}(T^{(k)}_{v},T^{(k)}_{u})$
- 2.
The graph embeddings are lower Lipschitz in expectation w.r.t. $TMD^{(K)}$
$\mathbb{E}_{\phi^{(K+1)}}[\|c^{\mathrm{global}}-\hat{c}^{\mathrm{global}}\|_{p
}^{p}]\geq c_{2}\cdot TMD^{(K)}(G,\hat{G})^{p}$

For the proof of the theorem we will need the following simple but useful lemma

###### Lemma F.4 .

Let $D,N$ be natural numbers and $p>1$. Then there exists a positive $\delta=\delta(D,N,p)$ such that, for all $N$ fixed vectors $x_{1},\ldots,x_{N}$ in $\mathbb{R}^{D}$

$$ $\mathbb{P}\left\{a\in S^{D-1}|\quad|a\cdot x_{i}|\geq\delta\|x_{i}\|_{p}, \forall i=1,\ldots,N\right\}\geq 1/2$ $$

###### Proof of Lemma F.4 .

Due to equivalence of norms C.1, it is sufficient to prove the claim when $p=2$.

For every $y\in S^{D-1}$ and $\delta>0$, denote

$$ $B(y,\delta)=\{a\in S^{D-1}|\quad|a\cdot y|<\delta\}$ $$

We note that for any fixed positive $\delta$ and $y,y^{\prime}\in S^{D-1}$, the probability of
$B(y,\delta)$ and $B(y^{\prime},\delta)$ will be the same, due to the rotation invariance of the uniform measure on $S^{D-1}$, and the fact that if $R$ is a rotation taking $y$ to $y^{\prime}$, then

$$ $a\in B(y,\delta)\text{ iff }|a\cdot y|<\delta\text{ iff }|Ra\cdot Ry|<\delta \text{ iff }Ra\in B(y^{\prime},\delta)$ $$

We can therefore denote the probability of $B(y,\delta)$ by $p_{\delta}$, and this definition does not depend on the choice of $y$. Next, note that $p_{\delta}$ is a non-negative seuence converging monotonely to 0 0 as $\delta\rightarrow 0$. Accordingly,we can choose some $\delta_{0}$ such that $p_{\delta}<\frac{1}{2N}$.

Now, assume we are given $N$ fixed points in $\mathbb{R}^{D}$. Without loss of generality we can assume the first $M$ points are non-zero, and the last $N-M$ points are zero. We then have

$$ $\displaystyle\mathbb{P}\left\{a\in S^{D-1}|\quad|a\cdot x_{i}|\geq\delta\|x_{i }\|_{2},\forall i=1,\ldots,N\right\}$ $\displaystyle=\mathbb{P}\left\{a\in S^{D-1}|\quad|a\cdot x_{i}|\geq\delta\|x_{ i}\|_{2},\forall i=1,\ldots,M\right\}$ $\displaystyle=\mathbb{P}\left\{a\in S^{D-1}|\quad|a\cdot\frac{x_{i}}{\|x_{i}\| _{2}}|\geq\delta,\forall i=1,\ldots,M\right\}$ $\displaystyle=1-\mathbb{P}\left(\cup_{i=1}^{M}B(\frac{x_{i}}{\|x_{i}\|_{2}}, \delta_{0})\right)$ $\displaystyle\geq 1-\sum_{i=1}^{M}\mathbb{P}\left(B(\frac{x_{i}}{\|x_{i}\|_{2} },\delta_{0})\right)$ $\displaystyle\geq 1-\frac{M}{2N}$ $\displaystyle\geq 1-\frac{N}{2N}=\frac{1}{2}.$ $$

∎

Using the lemma, we can now prove the theorem:

###### Proof of Theorem F.3 .

In general, to show that a function $f:X\times W\to Y$ is lower-Lipschitz in expectation, it suffices to show that there exist $\delta,\epsilon>0$ such that for all $x_{1},x_{2}\in X$, the probability of the set $\{w\in W:\quad\frac{|f(x_{1},w)-f(x_{2},w)|}{d_{X}(x_{1},x_{2})}>\epsilon\}$ is larger than $\delta$. We will use this alternative requirement in the proof of both parts of the theorem.

###### Theorem F.3 .

###### Lemma F.4 .

###### Proof of Lemma F.4 .

###### Proof of Theorem F.3 .

##### Proof of First Claim

For ease of notation, we will use $W$ throughout this proof to denote $W_{TD,1}^{\bar{T}_{z}{}}$, and $W^{p}(\cdot,\cdot)=W_{TD,1}^{\bar{T}_{z}{}}(\cdot,\cdot)^{p}$.

We prove the claim be induction on $k$. For $k=0$ we have equality

$$ $\|x_{v}^{(0)}-x_{u}^{(0)}\|^{p}=TD^{p}(T_{v}^{0},T_{u}^{0})$ $$

We now assume correctness for $k-1$ and prove for $k$. Now note that

- 1.
With probability of at least $p_{k-1}$ on $\theta^{(k-1)}$, we know that for all nodes $u,v$
$\|x_{v}^{(k-1)}-x_{u}^{(k-1)}\|^{p}\geq c_{k}TD^{p}(T_{v}^{(k-1)},T_{u}^{(k-1)})$
- 2.
Once $\theta^{(k-1)}$ is fixed, and all features $x_{v}^{(k-1)}$ are determined, we can use Lemma F.4 to show that with an appropriate $\delta>0$ and probability of at least $1/2$ on $a^{(k)}$,
$|a^{(k)}\cdot(x-y)|\geq\delta\|x-y\|$
for all $x,y$ in the set $\{x_{v}^{(k)}\}_{v\in V}\cup\{z_{k}\}$. It follows that, with probability $1/2p_{k-1}$ on $(a^{(k)},\theta^{(k-1)})$
$\displaystyle\|c_{v}^{(k)}-c_{u}^{(k)}\|_{p}^{p}$
$\displaystyle=\|\mathrm{sort}\left(a^{(k)}\cdot x_{v,\bullet}^{(k-1)}\right)-
\mathrm{sort}\left(a^{(k)}\cdot x_{u,\bullet}^{(k-1)}\right)\|_{p}^{p}$
$\displaystyle=\min_{\tau\in S_{n}}\sum_{s=1}^{n}|a^{(k)}\cdot(x^{(k-1)}_{v,s}-
x^{(k-1)}_{u,\tau(s)})|^{p}$
$\displaystyle\geq\min_{\tau\in S_{n}}\sum_{s=1}^{n}\delta^{p}\|x^{(k-1)}_{v,s}
-x^{(k-1)}_{u,\tau(s)}\|_{2}^{p}$
$\displaystyle\geq\delta^{p}C_{1}^{p}\min_{\tau\in S_{n}}\sum_{s=1}^{n}\|x^{(k-
1)}_{v,s}-x^{(k-1)}_{u,\tau(s)}\|_{p}^{p}$
$\displaystyle\stackrel{{\scriptstyle(*)}}{{\geq}}\delta^{p}C_{1}^{p}\tilde{c}_
{k}\min_{\tau\in S_{n}}TD^{p}(T_{v,s}^{(k-1)},T_{u,\tau(s)}^{(k-1)})$
$\displaystyle\geq\delta^{p}C_{1}^{p}\tilde{c}_{k}C_{2}\left[\min_{\tau\in S_{n
}}TD(T_{v,s}^{(k-1)},T_{u,\tau(s)}^{(k-1)})\right]^{p}$
$\displaystyle=c_{k,1}W^{p}(\mathcal{T}_{v}^{(k-1)},\mathcal{T}_{u}^{(k-1)})$
where in the last equation we used $c_{k,1}$ to denote the product of all constants appearing previously. (*) is the induction hypothesis and lemma F.1.
- 3.
Once $(a^{(k)},\theta^{(k-1)})$ are determined, we know that $c_{v}^{(k)}$ and $x_{v}^{(k-1)}$ are determined. Thus, for an appropriate positive $\delta^{\prime}$, we know that
with probability of at least $1/2$ on $d^{(k)}$, and thus with probability of at least $p_{k}:=p_{k-1}/4$ on $\theta^{(k)}=(d^{(k)},a^{(k)},\theta^{(k-1)})$, we have
$\displaystyle|x_{v}^{(k)}-x_{u}^{(k)}|^{p}$
$\displaystyle=|d^{(k)}\cdot concat(x_{v}^{(k-1)}-x_{u}^{(k-1)},c_{v}^{(k)}-c_{
u}^{(k)})|^{p}$
$\displaystyle\geq(\delta^{\prime})^{p}\|concat(x_{v}^{(k-1)}-x_{u}^{(k-1)},c_{
v}^{(k)}-c_{u}^{(k)})\|_{p}^{p}$
$\displaystyle\geq(\delta^{\prime})^{p}\left(\|x_{v}^{(k-1)}-x_{u}^{(k-1)}\|_{p
}^{p}+\|c_{v}^{(k)}-c_{u}^{(k)}\|_{p}^{p}\right)$
$\displaystyle\geq(\delta^{\prime})^{p}c_{k-1}TD^{p}(T_{v}^{(k-1)},T_{u}^{(k-1)
})+(\delta^{\prime})^{p}c_{k,1}W^{p}(\mathcal{T}_{v}^{(k-1)},\mathcal{T}_{u}^{
(k-1)})$
$\displaystyle\geq(\delta^{\prime})^{p}\min\{c_{k-1},c_{k,1}\}\|x_{v}^{(0)}-x_{
u}^{(0)}\|_{p}^{p}+W^{p}(\mathcal{T}_{v}^{(k-1)},\mathcal{T}_{u}^{(k-1)})$
$\displaystyle\stackrel{{\scriptstyle(*)}}{{\geq}}C(\delta^{\prime})^{p}\min\{c
_{k-1},c_{k,1}\}\left(\|x_{v}^{(0)}-x_{u}^{(0)}\|_{p}+W(\mathcal{T}_{v}^{(k-1)
},\mathcal{T}_{u}^{(k-1)})\right)^{p}$
$\displaystyle=c_{k}TD^{p}(T_{v}^{(k)},T_{u}^{(k)})$
where in the last equation we used $c_{k}$ to denote all constants incurred up to this step. In the inequality (*) we used equivalence of norms in euclidean space, with an appropriate constant $C$.

This concludes the proof of the first part of the theorem.

##### Proof of Second Claim

Let us denote the final message passing features corresponding to the nodes of $G$ by $x_{v}^{(K)}$, and the depth-K computation trees corresponding to each node of $G$ by $T_{v}^{(K)}$. Similarly, we denote the features corresponding to the nodes of $\hat{G}$ by $\hat{x}_{v}^{(K)}$ and the trees by $\hat{T}_{v}^{(K)}$. By applying the first part of the theorem (to the disjoint union of the graphs $G$ and $\hat{G}$), we know that with probablity of at least $p_{K}$ on the parameters $\theta_{K}$, we have that

$$ $|x_{v}^{(K)}-\hat{x}_{u}^{(K)}|^{p}\geq c_{k}TD^{p}(T_{v}^{(K)},\hat{T}_{u}^{( K)})$ $$

Once $\theta_{K}$ is fixed, all node features $x_{v}^{(K)}$ and $\hat{x}_{u}^{(K)}$ are determined. By Lemma F.4, we have for an appropriate $\delta>0$ that, with probablity of at least $1/2$ on $a^{(K+1)}$,

$$ $|a^{(K+1)}\cdot(x_{v}^{(K)}-\hat{x}_{u}^{(K)})|\geq\delta\|x_{v}^{(K)}-\hat{x} _{u}^{(K)}\|_{2},\quad\forall u,v\in[n].$ $$

It follows that for an appropriate $\delta^{\prime}$, with probability of at least $p_{k}/4$ on $\theta_{K+1}=(b^{(K+1)},a^{(K+1)},\theta_{K})$

$$ $\displaystyle|c^{\mathrm{global}}-\hat{c}^{\mathrm{global}}|^{p}$ $\displaystyle=|b^{(K-1)}\cdot[\mathrm{sort}\left(a^{(K+1)}\cdot x_{\bullet}^{( K)}\right)-\mathrm{sort}\left(a^{(K+1)}\cdot\hat{x}_{\bullet}^{(K)}\right)]|^{p}$ $\displaystyle\geq\delta^{\prime p}\|\mathrm{sort}\left(a^{(K+1)}\cdot x_{ \bullet}^{(K)}\right)-\mathrm{sort}\left(a^{(K+1)}\cdot\hat{x}_{\bullet}^{(K)} \right)\|_{p}^{p}$ $\displaystyle=\delta^{\prime p}\min_{\tau\in S_{n}}\sum_{v=1}^{n}|a^{(K+1)} \cdot x_{v}^{(K)}-a^{(K+1)}\cdot\hat{x}_{\tau(v)}^{(K)}|^{p}$ $\displaystyle\geq(\delta\cdot\delta^{\prime})^{p}\min_{\tau\in S_{n}}\sum_{v=1 }^{n}\|x_{v}^{(K)}-x_{\tau(u)}^{(K)}\|_{p}^{p}$ $\displaystyle\stackrel{{\scriptstyle(*)}}{{\geq}}(C\cdot\delta\cdot\delta^{ \prime})^{p}\left(\min_{\tau\in S_{n}}\sum_{v=1}^{n}\|x_{v}^{(K)}-x_{\tau(u)}^ {(K)}\|_{p}\right)^{p}$ $\displaystyle\geq(C\cdot\delta\cdot\delta^{\prime})^{p}TMD^{(K)}(G,\hat{G})$ $$

where for $(*)$ we used equivalence of norms in Euclidean spaces, and the last inequality uses the first claim and lemma F.1.

∎

### F.2 MPNN with lower-Hölder components isn’t necessarily lower-Hölder

Consider the following setting: we consider graphs in $G_{\leq n}(\Omega)$ where $\Omega=[-2,2]\subseteq\mathbb{R}$. Assume we have an MPNN $f$ with two layers, where the aggregation functions are the one dimensional functions $q:(x_{1},\ldots,x_{n})\mapsto\sum_{i=1}^{n}\mathrm{ReLU}(x_{i}-b)$ where $b\sim[-3,3]$. We showed this type of multiset aggregation is Lower-Hölder in expectation C.4. We claim that with these one-dimensional aggregations and for any COMBINE function, there are graphs $G_{1},G_{2}\in G_{\leq n}(\Omega)$ which can be separated by two iterations of 1-WL, but cannot be separated by two of our MPNN iterations, for any choice of network parameters. An illustration of these graphs is given in Figure 6. The parameter $\epsilon$ can be taken to be any fixed number in $(0,1/2)$.
As can be seen, the depth-2 computation trees of the nodes at level 2 (children of the root, filled in red) differ between the graphs. Therefore, the 1-WL test will succeed in determining the graphs are non-isomorphic after two iterations.

We now explain why the aforementioned MPNN won’t succeed in separating the graphs in two iterations. In our analysis, we use the names of the nodes defined in the figure: $a_{1},\ldots,e_{1}$ for the nodes of $G_{1}$ and $a_{2},\ldots,e_{2}$ for the nodes of $G_{2}$. We denote that feature vector at (e.g.,) $a_{1}$ after the $k$-th iteration by $a_{1}^{(k)}$.

The core reason for the failure of two iterations of the MPNN to separate $G_{1},G_{2}$ is that, depending on the value of $b$, the first aggregation will not be able to simultaneously separate the multisets $\{\!\!\{-\epsilon,\epsilon\}\!\!\}$ and $\{\!\!\{0,0\}\!\!\}$ and the multisets $\{\!\!\{1-\epsilon,1+\epsilon\}\!\!\}$ and $\{\!\!\{1,1\}\!\!\}$. In general, we can consider three options, depending on the value of $b$:

- 1.
$b\in[-3,-\epsilon]\cup[\epsilon,1-\epsilon]\cup[1+\epsilon,3]$: In this case, the aggregation won’t separate the multisets $\{\!\!\{-\epsilon,\epsilon\}\!\!\}$ and $\{\!\!\{0,0\}\!\!\}$, and it also won’t separate the multisets $\{\!\!\{1-\epsilon,1+\epsilon\}\!\!\}$ and $\{\!\!\{1,1\}\!\!\}$.
In this case, after the first message passing layer, we will get equality between the corresponding features of nodes $a_{1}^{(1)}=a_{2}^{(1)}=c_{1}^{(1)}=c_{2}^{(1)},b_{1}^{(1)}=b_{2}^{(1)}=d_{1}^
{(1)}=d_{2}^{(1)}$. This means the multisets $\{\!\!\{a_{1}^{(1)},b_{1}^{(1)}\}\!\!\}$ and $\{\!\!\{a_{2}^{(1)},b_{2}^{(1)}\}\!\!\}$ will be equal, and so will the multisets $\{\!\!\{c_{1}^{(1)},d_{1}^{(1)}\}\!\!\}$ and $\{\!\!\{c_{2}^{(1)},d_{2}^{(1)}\}\!\!\}$. Therefore, the second message passing iteration will result in the equality between the corresponding features of nodes $e_{1}^{(2)}=e_{2}^{(2)},f_{1}^{(2)}=f_{2}^{(2)}$, and the graphs won’t be separated.
- 2.
$b\in(-\epsilon,\epsilon)$: In this case the aggregation can separate the multisets $\{\!\!\{-\epsilon,\epsilon\}\!\!\}$ and $\{\!\!\{0,0\}\!\!\}$, but not the multisets $\{\!\!\{1-\epsilon,1+\epsilon\}\!\!\}$ and $\{\!\!\{1,1\}\!\!\}$.
In this case, after a single message passing layer, we will get equality between features of the nodes $a_{1}^{(1)}=a_{2}^{(1)},c_{1}^{(1)}=c_{2}^{(1)},d_{1}^{(1)}=b_{2}^{(1)}=b_{1}^
{(1)}=d_{2}^{(1)}$.
This means the multisets $\{\!\!\{a_{1}^{(1)},b_{1}^{(1)}\}\!\!\}$ and $\{\!\!\{a_{2}^{(1)},b_{2}^{(1)}\}\!\!\}$ will be equal, and so will the multisets $\{\!\!\{c_{1}^{(1)},d_{1}^{(1)}\}\!\!\}$ and $\{\!\!\{c_{2}^{(1)},d_{2}^{(1)}\}\!\!\}$. Therefore, the second message passing iteration will result in the equality between the corresponding features of nodes $e_{1}^{(2)}=f_{2}^{(2)},f_{1}^{(2)}=e_{2}^{(2)}$, and the graphs won’t be separated.
- 3.
$b\in(1-\epsilon,1+\epsilon)$: In this case the aggregation can separate the multisets $\{\!\!\{1-\epsilon,1+\epsilon\}\!\!\}$ and $\{\!\!\{1,1\}\!\!\}$, but not the multisets $\{\!\!\{-\epsilon,\epsilon\}\!\!\}$ and $\{\!\!\{0,0\}\!\!\}$.
This will lead to a similar result as the previous case, just with a slight change to the equivalence classes of the nodes.

Overall, we have shown that the proposed MPNN, which is composed of lower-Hölder in expectation components, fails to separate the above graphs for any choice of parameters. Therefore, it cannot be lower-Hölder in expectation w.r.t. $TMD^{(2)}$.

We note that the graphs can be separated by a similar ReLUMPNN which is wider, or has additional message passing layers. Our point here is just that there exists an MPNN with $K$ message passing layers, composed of lower-Hölder in expectation components, which isn’t lower Hölder in expectation w.r.t. $TMD^{K}$.

Figure: Figure 6: A pair of labeled graphs $G_{1}$ (top) and $G_{2}$ (bottom) used to prove that MPNNs with Hölder building blocks aren’t necessarily Hölder
Refer to caption: extracted/5914493/figs/non_holder_mpnn_example.png

### F.3 ReluMPNN and SmoothMPNN

In the main text of the paper we stated the following informal theorem
\lb*
We note that the distribution on the parameters taken in this theorem is the same as described throughout the text. The aggregation used for both MPNNs at hand is of the form

$$ $\{\!\!\{x_{1},\ldots,x_{s}\}\!\!\}\mapsto\sum_{i=1}^{s}\rho(a\cdot x_{i}-b)$ $$

where $a$ and $b$ are drawn uniformly from the unit sphere and some interval $[-B,B]$ respectively. Of course, in ReluMPNN $\rho=\mathrm{ReLU}$ while with smoothMPNN $\rho$ is a smooth function. The parameters of the linear COMBINE function are also selected uniformly in the unit sphere.

##### ϵ italic-ϵ \epsilon italic_ϵ -Tree dataset

The proof of the theorem is based on a family of adversarial example which we denote by $\epsilon$-trees. The recursive construction depends on two parameters: a real number $\epsilon$, and an integer $T=0,1,\ldots$. This will give us a pair of trees $(G,\hat{G})$ where $G=G(T,\epsilon),\hat{G}=\hat{G}(T,\epsilon)$.

The recursive construction of $\epsilon$ trees is depicted in Figure 7. In the first step $T=0$, we define ’trees’ $a_{0},b_{0},c_{0},d_{0}$ which consist of a single node, with feature value of $\epsilon,-\epsilon,2\epsilon$ and $-2\epsilon$, respectively. We connect $a_{0},b_{0}$ via a new root with feature value 0 0 to create the tree $G_{0}=G(T=0,\epsilon)$, as depicted in the left hand side of Figure 7. We do the same thing with $c_{0},d_{0}$ to obtain the tree $\hat{G}_{0}=\hat{G}(T=0,\epsilon)$.

We now proceed recursively, again as shown in Figure 7. Assuming that for a given $T$ we have defined the trees $a_{T},b_{T},c_{T},d_{T},G_{T},\hat{G}_{T}$, we define the corresponding trees for $T+1$ as shown in the figure:

For example, $a_{T+1}$ is constructed by joining together two copies of $a_{T}$ and two copies of $b_{T}$ at a common root with feature value 0 0 as shown in the figure on the right hand side. From the figure one can also infer the construction of $b_{T+1},c_{T+1},d_{T+1}$ from $a_{T},b_{T},c_{T},c_{T}$. Finally, the tree $G_{T+1}$ is constructed by joining together $c_{T+1},d_{T+1}$ while $\hat{G}_{T+1}$ is constructed by joining together $a_{T+1},b_{T+1}$

###### Proof of Theorem 4.2 .

Denoting an MPNN of depth $T$ by $f_{T}(G,w)$, where $w$ denotes the network’s parameters, our goal is to provide bounds of the form

$$ $\mathbb{E}_{w}\|f_{T}(G(T,\epsilon),w)-f_{T}(\hat{G}(T,\epsilon),w)\|^{p}\leq C \epsilon^{p\beta}$ (7) $$

where $\beta$ depends on $T$ and the activation used, and on the other had show that $TMD^{T}(G(T,\epsilon),\hat{G}(T,\epsilon))$ scales linearly in $\epsilon$. This will prove that the MPNN cannot be $\alpha$-lower Hölder with any $\alpha<\beta$.

We can see that $TMD^{T}(G(T,\epsilon),\hat{G}(T,\epsilon))$ scales linearly in $\epsilon$ because the TMD is
homogenous. Namely, given any $G=(V,E,X)$ and
$\epsilon>0$, let $\epsilon G=(V,E,\epsilon X)$.
Then for any $T$ we have

$$ $TMD^{T}(\epsilon G,\epsilon\hat{G})=\epsilon\cdot TMD^{T}(G,\hat{G}).$ $$

We now turn to estimate the exponent $\beta$ in (7), for smooth or ReLU activations, and different values of $T$. To avoid cluttered notation we will explain what occurs for low values $T=0$ and $T=1$ in detail, and outline the argument for larger $T$.

When $T=0$, we apply an MPNN of depth $T=0$ to the graphs $G_{0},\hat{G}_{0}$, which means that we just apply a readout function to the initial features. This is exactly the $\pm\epsilon$ example discussed in the main text, and so (7) will hold with $\beta=(p+1)/p$ with ReLU activations, and $\beta=2$ with smooth activations.

We now consider the case $T=1$. Here we apply we an MPNN of depth $T=1$ to the graphs $G_{1},\hat{G}_{1}$. Let us denote the nodes of $G_{1}$ by $V_{1}$, using the natural correspondence defined by the figure (middle) we can think of $V_{1}$ as the nodes of $\hat{G}_{1}$ as well. We denote the node feature values at node $v$ after a single MPNN iteration by $x_{v,w}$ and $\hat{x}_{v,w}$, where $w$ denotes the network parameters as before.

Note that for the root $r$ we have that $x_{r,w}=\hat{x}_{r,w}$ for any choice of parameters $w$.

At initialization, all leaves of the trees $G_{0},G_{1}$
are assigned one of the values $\epsilon,-\epsilon,2\epsilon,-2\epsilon$ (according to the labeling $a_{0},b_{0},c_{0}$ or $d_{0}$). If two leaves are assigned the same value (both are denoted by, say, $a_{0}$), then they will have the same node feature after a single message passing iteration, since they are only connected to their father, and all fathers have the same initial node feature.

We now consider the two remaining nodes. We denote the node connected to the root from the left by $v_{1}$ and the node connected to the root from the right by $v_{2}$.

Let us first focus on the smooth activation case. We note that the children of $v_{1}$ and $v_{2}$, in both graphs, have initial features of order $\epsilon$, and the all sum to the same value zero. By considering the Taylor approximation of the activation, as discussed in the proof of Theorem 2, we see that for an appropriate constant $C$ we have that for all small enough $\epsilon$,

$$ $|c_{v,w}-\hat{c}_{u,w}|\leq C\epsilon^{2},\text{ for all parameters }w$ (8) $$

where $v$ and $u$ could either be $v_{1}$ or $v_{2}$, and we recall that $c_{v,w}$ is the output of the aggregation function, applied to the initial features of the neighbors of $v$ in $G$ (and $\hat{c}_{u,w}$ is defined analogously for the neighbors of $u$ in $\hat{G}$).

Additionally, we note that for every parameter $w$ we have that

$$ $c_{v_{1},w}+c_{v_{2},w}=2\left[\psi_{w}(\epsilon)+\psi_{w}(-\epsilon)+\psi_{w} (2\epsilon)+\psi_{w}(-2\epsilon)\right]=\hat{c}_{v_{1},w}+\hat{c}_{v_{2},w}$ (9) $$

where $\psi_{w}$ denotes the neural network used for aggregation. The next step in the MPNN procedure is the COMBINE step, after which we will have a bound of the form (9), namely

$$ $|x_{v,w}-\hat{x}_{u,w}|\leq C\epsilon^{2},\text{ for all parameters }w$ (10) $$

since the COMBINE functions we use are uniformly Lipschitz, and we will also have
or every parameter $w$ we have that

$$ $x_{v_{1},w}+x_{v_{2},w}=\hat{x}_{v_{1},w}+\hat{x}_{v_{2},w}$ (11) $$

since our COMBINE function is linear.

Now, when applying the readout function, we have eleven node features in the multiset of nodes of $G_{1}$ and $\hat{G}_{1}$. For any parameter vector $w$, if we remove the nodes $v_{1},v_{2}$, we get two identical multisets. Therefore, we only need to bound the difference

$$ $\|\left[\eta_{w}(x_{v_{1},w})+\eta_{w}(x_{v_{2},w})\right]-\left[\eta_{w}(\hat {x}_{v_{1},w})+\eta_{w}(\hat{x}_{v_{2},w}\right]\|$ $$

Since by (11) the features sum to the same value, and all are the same up to $\epsilon^{2}$, this difference goes like $\epsilon^{4}$. This is what we wanted. For general $T$, we will be able to apply this same argument recursively $T+1$ times to get the lower bound of $2^{T+1}$ for the smooth exponent.

Now, let us consider what happens when we replace the smooth activation with ReLU activation. In This case too, the only features which really matter are those corresponding to $v_{1},v_{2}$. Recalling our analysis of the $\pm\epsilon$ for ReLU activations, we note that if the bias of the network falls outside the interval $[-2\epsilon,2\epsilon]$, then all node features $x_{v_{1},w},x_{v_{2},w},\hat{x}_{v_{2},w},\hat{x}_{v_{2},w}$ will all be identical. If the bias does fall in that interval, then the difference between these node features can be bounded from above by $C\epsilon$ for an appropriate constant $C$. When applying readout, there is a probability of $\sim\epsilon^{2}$ that the final global features will be distinct (if the parameters of aggregation, and the parameters of readout, are both well behaved). The difference in this case can again be bounded by $C\epsilon$. It follows that

$$ $\mathbb{E}_{w}\|f_{T}(G(T,\epsilon),w)-f_{T}(\hat{G}(T,\epsilon),w)\|^{p}\leq C \epsilon^{2}\epsilon^{p}$ $$

which will lead to a lower bound on the Hölder constant $\alpha\geq(1+\frac{2}{p})$.

For general $T$, we will get that there is a probability of $\sim\epsilon^{T+1}$ to obtain any separation, and that, if this separation occurs, the difference between features will be $\leq C\cdot\epsilon$. As a result, the expectation

$$ $\mathbb{E}_{w}\|f_{T}(G(T,\epsilon),w)-f_{T}(\hat{G}(T,\epsilon),w)\|^{p}$ $$

will scale like $\epsilon^{T+1+p}$, which leads to our lower bound on the exponent for ReLU activations.
∎

Figure: Figure 7: $\epsilon$-tree dataset construction: The trees are built using basic elements $a_{t},b_{t},c_{t},d_{t}$. At step $T$, the basic elements are defined using subtrees from the previous step.
Refer to caption: extracted/5914493/figs/eps_tree_build.png

###### Proof of Theorem 4.2 .

## Appendix G Experiments

### G.1 Hardware

All experiments were executed on an a cluster of 8 nvidia A40 49GB GPUs.

### G.2 Architecture variations

In this subsection we will lay out the architectural variations that were considered for SortMPNN and AdaptMPNN, that were omitted from the main text for brevity. We include explanations of these variations with the connection between them and the theoretical analysis.

##### SortMPNN

First off, we must choose how to model of the blank vector, since SortMPNN uses the blank vector explicitly in (2). We propose three options for the (blank method):

- 1.
Learnable: For every multiset aggregation function there is a unique learnable blank vector that is used for augmentation. Since we assume the node features are from a compact subset of $\mathbb{R}^{d}$, there exist valid augmentation vectors, and upon training we can expect the learnable vector to converge to a valid augmentation vector. Note that this method assumes training takes place.
- 2.
Iterative update: In this method, the blank vector of each layer is the output of the previous MPNN layers on the blank tree, where the blank tree node feature is set to zero. This works under two assumptions: (1): The nodes in the dataset don’t have the feature zero. (2): The intermediate layer outputs on computation trees from the data don’t clash with the same layer output on the blank tree. This seems relatively reasonable since SortMPNN is bi-Lipschitz in expectation, which in expectation should lead to injectivity.
- 3.
Zero: simply using the zero vector as the blank vector. This works under the assumption that the zero vector is neither a valid initial feature, nor a valid output of intermediate layers. This is the weakest of the three, but was experimented with nonetheless.

Given a multiset $X$, we experimented with two variations for the aggregation implementation. Both share step (1), but differ in step (2):

- 1.
projection: $X\mapsto colsort(A\rho_{(z)}(X)):=Y\in\mathbb{R}^{m\times n}$, where $A\in\mathbb{R}^{m\times d}$.
- 2.
collapse: In this step, we either perform (1) a matrix collapse: $rowsum(B\odot Y)$ where $\odot$ is the element-wise product and $B\in\mathbb{R}^{m\times n}$ or (2) a vector collapse: $Yb$, for $b\in\mathbb{R}^{n}$.

Note that when choosing to use the matrix collapse, under the assumption that the blank method provides a valid blank vector, this is equivalent to running $m$ separate copies of $S_{z}$ (2) with independent parameters (at least upon initialization), just like we proposed in the main text. When using the vector collapse, the aggregation takes the form of multiple instances of $S_{z}$ except that that the parameter $b$ used in the inner product is shared across instances. Note the Hölder expectation doesn’t change in this case, but the variance is likely higher. We explored this option since it reduces the number of parameters needed for the aggregation, potentially easing the optimization process.

In addition, we experimented with adding bias to the projection, and to the vector collapse.

##### AdaptMPNN

The main detail we left out from the main text is the fact that the output of $m_{\mathrm{ReLU}}^{\mathrm{adapt}}$ is in $\mathbb{R}^{4}$, whereas we would want it to be in $\mathbb{R}$ in order to be able to stack $m$ instances and get an output in $\mathbb{R}^{d}$. In order to get the desired output dimension, we compose $m_{\mathrm{ReLU}}^{\mathrm{adapt}}$ with a projection, and the aggregation on a multiset $X\in\mathbb{R}^{d\times r}$ as follows:

- 1.
Stacked Adaptive ReLU:
$Y:=\begin{bmatrix}m_{\mathrm{ReLU}}^{\mathrm{adapt}}(X;a_{1},t_{1})\\
\vdots\\
m_{\mathrm{ReLU}}^{\mathrm{adapt}}(X;a_{m},t_{m})\end{bmatrix}\in\mathbb{R}^{m
\times 4}$
- 2.
Project: $rowsum(B\odot Y)$ where $\odot$ is the element-wise product and $B\in\mathbb{R}^{m\times 4}$ has rows that are drawn uniformly independently from $S^{3}$.

As we proved previously C.3, an inner product with a vector drawn uniformly from $S^{d-1}$ is bi-Hölder in expectation. From composition properties B.1, this means that the above stacking of multiple independent instances of $project\circ m_{\mathrm{ReLU}}^{\mathrm{adapt}}$ keeps the Hölder properties of $m_{\mathrm{ReLU}}^{\mathrm{adapt}}$.

In addition, we experimented with four optional changes to the aggregation:

- 1.
Adding the sum $\sum a\cdot x_{j}$ a 5’th coordinate to the output of $m_{\mathrm{ReLU}}^{\mathrm{adapt}}$. The idea behind this choice is that the sum is an informative feature which is encountered when using the standard ReLU activations and bias smaller than the minimal feature in the multiset. For adaptive ReLU, this summation will only occur when $t=0$.
- 2.
Note that the adaptive relu uses the parameter $t\in[0,1]$ to choose a bias within the minimal and maximal multiset values. When training, we experimented with either clamping $t$ to $[0,1]$ or otherwise lifting this constraint.
- 3.
We experimented with adding a bias term to the projection step, since this can potentially lead to stronger expressivity.
- 4.
the initialization of $t$ was optionally chosen to be linspace between $[0,1]$.

### G.3 Adversarial Multiset dataset experiments

Figure 1 and figure 2 both made use of the adversarial multiset pairs which are constructed as is described in detail in section C.2.

To produce figure 1, we ran $m_{relu},m_{\sigma}$ on a single adversarial pair of multisets that have $8$ scalar elements, with $7$ equal moments between the multisets, with $\epsilon=0.1$.

To produce figure 2, we used $200$ adversarial pairs, each containing $16$ scalar elements, where each pair shares the first $15$ moments. The pairs differ in the value of $\epsilon$, where $\epsilon\in[0,1]$

### G.4 ϵ italic-ϵ \epsilon italic_ϵ -tree dataset

##### Distortion experiment

The experiment depicted in 3 was run on a set of $1,000$ pairs of trees with $\epsilon$ going from $0.05$ to $0.4$. The first $500$ trees are of height $h=3$ and were input to MPNNs of depth $d=1$, and the other $500$ trees were of height $h=4$ and were input to MPNNs of depth $d=2$. The plots in figure 3 were separated by MPNN depth in order for the phenomenon to be clearly visible. We note that the feature of all non-leaf nodes which appear with 0 0 in figure 7, was set to be $1$. All the models in the experiment used linear combination as the combine function, had an embedding dimension of $45,000$, except SortMPNN that used an embedding dimension of $2,048$ due to resource constraints. The large width was used to show the behavior of the MPNNs as we approach the expected value over the parameters. In order to plot the curve in each subplot, the minimal ratio $r_{min}=l_{2}/TMD^{\alpha}$ was computed for each graph pair where $\alpha$ is the theoretical lower-Hölder exponent. Then, $r_{min}\cdot x^{\alpha}$ was plotted. Note that $r_{min}$ was only computed for half of the pairs with the largest TMD values in order to avoid the case where the embedding distance is zero leading to $r_{min}=0$.

Throughout the experiment, double precision was used in order to capture the small differences as needed. The measured embedding distance was taken over the output of the readout function, without further processing.

##### Separation experiment

For the trained graph separation experiment, $100$ pairs of $\epsilon$ trees of height $h=4$ with $\epsilon\in[0.1,1]$ were used. Trees corresponding to the first row in figure 7 were given the label 0 0, and trees from the second row were given the label $1$. The models were trained by inputting a single pair as a batch per training iteration. All models were composed of 2 message passing layers, followed by readout and a single linear layer to obtain a logit. All of $\{64,128,256,512\}$ were tested for the embedding dimension with results consistent regardless of the embedding dimension. The models were trained for 10 epochs using the adamw optimizer and a learning rate of 0.01.

### G.5 TUDataset

Results on datasets that appear in table 3 for GIN, GCN, GraphSage are from (Xu et al., ), except for GIN on NCI109 which is from (Chuang & Jegelka, ). The reported results in 3 for the fully trained SortMPNN and AdaptMPNN were chosen out to be the best single result out of 30 random hyper-parameter choices from the following hyperparameters:
batch size$\in\{32,64\}$, depth$\in[2-5]$, embed dim$\in\{16,32,64\}$, combine$\in\{ConcatProject,LinearCombination,LTSum\}$, dropout$\in\{0,0.1,0.2\}$,output mlp depth $\in\{1,2,3\}$, weight decay $\in\{0,0.01,0.1\}$, lr$\in[0.0001,0.01]$, optimizer $\in[adam,adamw]$, #epochs$=500$. In addition, for SortMPNN the blank method was taken to be iterative update, and the collapse method was taken to be from {matrix, vector} without bias begin used at any stage. For AdaptMPNN, adding the sum of the multiset was a hyper-parameter, as was the choice if to clamp the parameter $t$ to be in $[0,1]$. The choice of the 30 random configurations was chosen with wandb (Biewald, ) sweeps.

As stated in the main text, results are reported as in (Xu et al., ). Namely, we report the mean and standard deviation of the validation performance over stratified 10-folds, where the validation curves of all 10 folds are aggregated, and the best epoch is chosen.

### G.6 LRGB

In both experiments from table 4, the models adhere to the 500K parameter limit. The hyper-parameters were tuned for a single seed, and then the configuration with the best validation results was rerun for 4 different seeds. We report the mean and standard deviation on the test set in table 2. The code from the official LRGB github repo was used to run these experiments, with the only change being the addition of our architectures. The reported results for other models are taken from (Tönshoff et al., ).

Note that any theoretical result we presented that made use of a uniform distribution on the unit sphere holds also for any $t$-normalized unit sphere, i.e. ${x|\|x\|=t}$, where $t\in R$, where only the multiplicative constants in the proofs will change. In this experiment, for SortMPNN we normalized the parameters sampled from the unit sphere as follows: (1) the parameter $a$ was sampled from a $\frac{1}{d}$-normalized sphere, (2) the parameter $b$ was sampled from a $\frac{1}{n}$-normalized sphere. This normalization was chosen through hyper-parameter tuning on the validation set.

##### peptides-func

For SortMPNN, the following hyper-parameters were used to achieve the results shown in the table 4: blank method - learnable, collapse method - vector, bias used both in projection and collapse, combine - ConcatProject, embedding dimension - 165, message passing layers - 4, final mlp layers - 3, dropout - 0.05, weight decay - 0.01, positional encoding - RWSE (Dwivedi et al., ), learning rate - 0.005, scheduler - cosine with warmup, optimizer - adamW. Number of parameters: 494865.

For AdaptMPNN - the following hyper-parameters were used: the bias parameter $t$ was clamped to stay within $[0,1]$, the sum of a multiset was added as a 5’th coordinate, bias added in project, collapse method - LinearCombination, embedding dimension - 220, message passing layers - 6, final mlp layers - 3, dropout - 0.05,weight decay - 0.01, positional encoding - RWSE, learning rate - 0.001, scheduler - cosine with warmup, optimizer - adamW. Number of parameters: 480915.

##### peptides-struct

For SortMPNN: blank method - learnable, collapse method - vector, bias not used both in projection or collapse, combine - LinearCombination, embedding dimension - 200, message passing layers - 6, final mlp layers - 3, dropout - 0, weight decay - 0.01, positional encoding - LapPE (Dwivedi & Bresson, ), learning rate - 0.005, scheduler - step, optimizer - adamW. Number of parameters: 493858.

For AdaptMPNN:the bias parameter $t$ was clamped to stay within $[0,1]$, the sum of a multiset was added as a 5’th coordinate, bias added in project, collapse method - LinearCombination, embedding dimension - 220, message passing layers - 6, final mlp layers - 3, dropout - 0,weight decay - 0.1, positional encoding - LapPE, learning rate - 0.001, scheduler - step, optimizer - adamW. Number of parameters: 480973.

##### small models

Additional experiments run on peptides struct were performed by running models with the exact same hyper-parameters aside from the width. Results are shown in figure 4. Exact width and number of parameters (including for the 500K budget) appear in table 6.

**Table 6: [width | #params] per model for each parameter budget**
| Budget | SortMPNN | AdaptMPNN | GCN |
| --- | --- | --- | --- |
| 500K | [200 | 493.8K] | [220 | 478.8K] | [235 | 488K] |
| 100K | [70 | 95.9K] | [95 | 98.6K] | [100 | 98.4K] |
| 50K | [48 | 47.3] | [65 | 49.3K] | [68 | 48.8K] |
| 25K | [34 | 24.9K] | [44 | 24.3K] | [46 | 24.2K] |
| 7K | [17 | 6.7K] | [22 | 6.7K] | [23 | 6.8K] |
| 1K | [3 | 922] | [5 | 915] | [5 | 982] |

### G.7 Subgraph aggregation networks

In this experiment, we made use of edge features. The incorporation of edge features was done via the concatenation of the edge and relevant node features, followed by a linear layer projecting the concatenated vector to the desired dimension. The projected vector was then passed through a ReLU activation.

Test set mean and standard deviation were computed over 10 random seeds. Hyper-parameters were chosen by randomly running 40 configurations with the DS edge deleted configuration, and choosing the best. the validation set was done using a single seed. The code provided with (Bevilacqua et al., ) was used to run these experiments, with the only change being the addition of our architectures.

##### SortMPNN

Both for the DS and DSS experiments, the model consisted of 5 message passing layers, used sum Jumping Knowledge (JK) from (Xu et al., ), Combine - ConcatProject, matrix collapse, bias was added to projection, blank method - learnable. The model was trained for 400 epochs with lr=0.01 and batch size of 128.

##### AdaptMPNN

Both for the DS and DSS experiments, the model consisted of 5 message passing layers, sum JK for DS and last for DSS, combine - ConcatProject, sum wasn’t added, $t$ was clapmed to $[0,1]$, bias was added throughout aggregation steps. The model was trained for 400 epochs, with lr of 0.005 for DS and 0.01 for DSS.

## Appendix H Runtime

**Table 7: Average time per epoch of various models on the peptides-struct dataset. Average taken over 250 epochs.**
| Model | GINE | GCN | SortMPNN | AdaptMPNN |
| --- | --- | --- | --- | --- |
| Avg. time per epoch (s) | 12.16 | 14.25 | 16.08 | 26.77 |

## Peer Reviews

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: Holder continuity is usually used for exponents in the range [0,1]. Functions that are Holder continuous with exponents larger than 1 are trivial constant functions. It is not clear to me why the main results with \alpha > 1 are interesting. Is it because it is only the lower bound and not the upper? 

The choice of the probability measure under which Holder in expectation is taken seems to be essential. However, I don't see this appearing in the discussion or even in the proof sketches of most of the results. This was confusing to me. Even in Euclidean space, it seems it would be challenging to choose the right measure. If the data satisfied the manifold-hypothesis, for example, then a function that well-separates the data could be constant almost everywhere with respect to Lebesgue measure and thus have no lower bound with Lebesgue measure, but would have some appropriate lower bound with a measure tuned to the data.

In the definition of the Wasserstein metric on line 197 there is a p-norm on the right hand side, but it is W^1. Is this correct?

The authors give a lot of pathological adversarial examples of graphs that are difficult to separate. How practical are these examples?
soundness: 3
strengths: Deep graph neural networks is an important area of research, and there are an enormous number of different architectures that have been proposed. This paper provides some valuable insights into which networks give better separation properties. The lower Holder estimates in expectation are novel contributions to the literature, as far as I am aware. The new MPNN architectures with better separation properties are a nice contribution.
summary: The authors study stability and separability properties of multi-set and graph neural networks. They do this by establishing Holder lower bounds on the networks, which give some guarantees on separability. They argue that uniform Holder lower bounds are too restrictive, and instead propose a notion of lower Holder in expectation, which is averaged over a probability distribution on the parameters of the network. The main results in the paper are bounds on the lower Holder exponent for different types of graph and multiset neural networks. They find that common architectures have poor exponents and separability properties and propose two new message passing neural networks with better separation properties. They provide a number of numerical experiments with their newly proposed architectures, and compare against existing methods. 

Some typos I noticed: "decays rapidly with the network’s depth .", "paces. , let W be some set of parameters", "function. we say that", "multisets (Zaheer et al., 2017) ,as well as MPNNs"
weaknesses: The lower bounds may be pessimistic, relying on pathological adversarial examples. It is not clear how practical the advances are, in that the worst case graphs may not be encountered in practice.

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 2
questions: Comments and questions

1. Definition 2.1 depends on choices of $\alpha$ and $p \geq 1$. Shouldn't we be talking about whether functions are $\alpha,p$ lower-Holder in expectation instead of just $\alpha$ lower-Holder in expectation?
2. The remarks on line 162-165 don't make sense to me. Consider: "Firstly, we note that if the set of $w$ for which $f (x; w)$ is $\alpha$ lower-Hölder has positive probability, then $f(x; w)$ will be $\alpha$ lower-Hölder in expectation." I don't see why this is true. If we take $p=1$ and move $c^p$ in Definition 2.1 to the denominator of the RHS, then the $w$ for which $f(x; w)$ is $\alpha$ lower-Holder produce a value inside the expectation which is $\geq 1$. Then all we can say is that the RHS expectation is atleast the probability mass of such $w$'s that satisfy the original lower-Holder condition -- this bound cannot even be greater than $1$, which is what the LHS would become after dividing by $c^p$.
3. Is the implication of the experiments in Figure 3 that the presented bounds are tight?
4. Does the adaptive ReLU activation provide any experimental benefits?
5. Am I correct in concluding from Theorem 3.3 that any smooth activation (e.g. sigmoid) cannot be $\alpha$ lower-Holder in expectation (otherwise $\alpha$ would have to be infinity)? This seems counterintuitive; consider an activation like $a^{-1} \log(1 + e^{ax})$. This function is smooth and thus cannot be $\alpha$ lower-Holder in expectation, but in the limit $a \to \infty$ approaches the ReLU function which is $\alpha=3/2$ lower-Holder in expectation.

Minor remarks:
* Reading just theorem 3.1, it is unclear how $a$ and $b$ are related to the rest of the thoerem. I suggest writing $m_{ReLU}(\ \cdot\ ; a,b)$ in the theorem statement.
* Figures should be exported as pdfs or pgfs in order to avoid text blurring.
* There is a line of work on the abillity of GNNs to discriminate signals from a spectral theory perspective which might be worth discussing:

[a] Gama, Fernando, Joan Bruna, and Alejandro Ribeiro. "Stability properties of graph neural networks." IEEE Transactions on Signal Processing 68 (2020): 5680-5695.

[b] Pfrommer, Samuel, Alejandro Ribeiro, and Fernando Gama. "Discriminability of single-layer graph neural networks." ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2021.
soundness: 4
strengths: Overall, this is a heavily theoretical work that I am not well-equipped to confidently review. I have given it my best effort, although I will likely defer to other reviewers' judgements.

1. I really appreciate the $\epsilon$-tree experiment. I originally had my doubts that an architecture being $\alpha$ lower-Holder in expectation over *random weight parameters* is that meaningful for the practical ability of networks to distinguish inputs. Random-weight networks hardly reflect how trained networks ultimately behave. This experiment provides examples of networks which are and aren't $\alpha$ lower-Holder in expectation and shows that this property is critical in a synthetic task.
2. The theory is extensive but generally well-presented and the headline results are clear (although I did not check the proofs in detail). Proof sketches are included where appropriate.
3. The considered problem is interesting, and I appreciate how the authors generalized their framework to both multiset analysis framework to graph neural networks.
4. It's a surprising result to me that relu sum networks are $\alpha=3/2$ lower-Holder in expectation. I would have originally guessed that they would be just $\alpha=1$ lower-Holder (Lipschitz) in expectation.
summary: This work analyzes the abillity for multiset and graph neural networks to separate inputs. Due to the non-injectivity of ReLU networks, the authors introduce the property of a network being Holder in expectation over random weights (Holderness is just an extension of Lipschitz continuity with an exponential factor $\alpha$). The authors then derive the lower-Holder exponents for a variety of models, and provide preliminary empirical evidence that lower Holderness is related to network performance on some synthetic tasks.
weaknesses: 1. Some of the experimental benefits are pretty marginal (especially in table 4, where the gains are smaller than the exerimental variance).
2. While the $\epsilon$-Tree experiment provides some evidence for the importance of $\alpha$ lower-Holderness in expectation over random weights, it just includes a few anecdotal examples of architectures which do and don't work. A really convincing experiment would consist of a scalar-parameterized collection of function classes (which are in turn parameterized). Say this scalar smoothly varied the resulting $\alpha$; I'd then hope to see that as $\alpha$ increased, performance increases.
3. The paper seems rushed, and has many typos. Consider lines 145-147 for example (there are many more such mistakes). The authors should thoroughly proofread their work.
4. Related, the paper presentation is somewhat haphazard, with many text-wrapped figures and tables; these really should be combined appropriately and top-aligned to the page with no text wrapping.

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: - Regarding the weaknesses: Are there any major difficulties in a Hölder stability analysis of AdaptMPNN compared to ReluMPNN? If so, a short paragraph explaining this might be interesting.
- I am a bit puzzled about the results for the adaptive ReLU architectures: Both in Figure 2 and Figure 3, the corresponding architectures seem to be almost perfectly lower Lipschitz on the data, behaving similarly to the sorted architectures. If I am not mistaken, however, this is not guaranteed by the theory that is derived. I am assuming that this is not a consequence of the bound from Theorem 3.2 not being tight, since it is stated in lines 308-309 that the adaptive ReLU architecture for multisets is *not* lower Lipschitz in expectation, as can be shown by some adversarial graphs. However, I couldn't find any further justification for this statement in the manuscript, as Theorem 3.2 only shows that $(p+1)/p$-lower Hölder continuity in expectation (which evaluates to $\alpha=3/2$ for the $\ell^2$ case), while Theorem 3.1 only refers to standard ReLU summation. Could you maybe point me to a proof/justification or add one, and generally comment a bit more on the performance of the adaptive ReLU methods in Figure 2 and 3?
- As also mentioned in the paper (lines 344-346), there are multiple WL distances that have been suggested so far. As expected lower-Hölderness also crucially depends on the chosen distance functions, is there any particular reason why you chose the TMD?
- An interesting question would be to analyze more quantitatively how the expected Hölder stability evolves over training, as the theory is only laid out for the initial parameters (and depends crucially on the initialization distributions). While the results of Table 2 suggest that the stability is preserved during training, I would find more quantitative empirical results concerning the evolvement of the Hölder stability with the NN weights during training very interesting to see. However, as this is very open-ended and potentially costly, I would like to emphasize that this is merely as a suggestion.
soundness: 3
strengths: - The problem is well-motivated and I enjoyed reading the paper. To the best of my knowledge, the idea is novel, and as part of the ongoing effort to understand the geometry of MPNN embeddings, this work is a notable step.
- The theory is made quite accessible, and I particularly found the $\pm \epsilon$ example, which is used throughout the theory section, to be very instructive for a high-level understanding of the theorems.
- The presented empirical results are strong and insightful, showing that AdaptMPNN and SortMPNN can perfectly distinguish certain pairs of adversarial WL-distinguishable graphs, a task which standard MPNNs fail entirely.
summary: The paper introduces a framework for quantifying the separation quality of multiset and graph neural networks (GNNs) via Hölder stability. This is motivated by recent observations that certain graph pairs, despite being theoretically distinguishable by Weisfeiler-Leman (WL)-expressive MPNNs (message passing GNNs), may have their embeddings "squashed" by GNNs, to the point that they are practically indistinguishable. The authors argue that strict lower Hölder continuity, which is only fulfilled for injective embeddings, is too restrictive as a concept and instead consider lower Hölder continuity in expectation over fixed parameter distributions. Starting with the multiset case, which appears as building blocks of MPNNs, they analyze the expected Hölder stability of popular architectures and introduce two new approaches, one based on bias scaling and the other on sorting, that improve the theoretical guarantees. They then employ these new multiset architectures in MPNNs and analyze their expected Hölder stabilities w.r.t. a WL distance, the tree mover's distance (TMD) (Chuang & Jegelka, 2022). Notably, the sort-based technique is shown to be lower Lipschitz in expectation w.r.t. the TMD. The new architectures AdaptMPNN and SortMPNN are empirically analyzed on toy and real-world datasets, demonstrating superior separation compared to popular MPNN baselines.
weaknesses: In my opinion, the main weakness of the paper in its current state is that there is no theoretical analysis of the (expected) lower-Hölder properties of AdaptMPNN. While AdaptMPNN is one of the main methodological contributions of the work, only the lower-Hölder stability of the two "baselines" ReluMPNN and SmoothMPNN as well as SortMPNN are analyzed in section 4. Ideally, an outline of a worst case analysis as mentioned at the end of section 4, or at least some justification if or why there is something inherently harder about the analysis of AdaptMPNN compared to ReluMPNN might be an insightful addition to the manuscript.

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 3
questions: 1. The authors claimed that "SortMPNN is the first MPNN with both upper and lower Lipschitz guarantees" (Line 120-121). However, BLIS-Net proposed by Xu et al. [1] is shown to be Bi-Lipschitz with respect to a weighted inner product space. Can the authors compare and comment the differences?

2. Theorem 3.4 (SortMPNN is uniformly upper Lipschitz
and lower Lipschitz in expectation): This result seems straightforward given results in Balan et al. (2022). Can the authors discuss this more explicitly?

3. Connections between SortMPNNs and other related models: While sort-based aggregation function is rarely used in practice (and thus the novelty of the proposed SortMPNN), we can still see its "shades" in some related works. For example, in [2] the authors experimented with a sorting-based aggregation function; in [3] the authors proposed to learn a recursive aggregation function, which could in theory mimic a recursive sorting algorithm. Can the authors comment more of these related models?


References:

[1] Xu, Charles, et al. "BLIS-Net: Classifying and Analyzing Signals on Graphs." International Conference on Artificial Intelligence and Statistics. PMLR, 2024.


[2] Duvenaud, David K., et al. "Convolutional networks on graphs for learning molecular fingerprints." Advances in neural information processing systems 28 (2015).


[3] Ong, Euan, and Petar Veličković. "Learnable commutative monoids for graph neural networks." Learning on Graphs Conference. PMLR, 2022.
soundness: 3
strengths: 1. The study of separation power and stability of functions on multisets and graphs is of high interest to the community.

2. The proposed notion of lower-Hölder in expectation is interesting and can explain practical performance.

3. The proposed adversarial examples (i.e., $\pm \epsilon$-multisets, $\epsilon$-trees) are instructive and very useful to illustrative the theory.
summary: This paper proposed a pairwise separation quality analysis framework to study the Lipschitz and Hölder stability of functions on multisets and graphs. The authors proposed the notion of lower-Hölder in expectation (Definition 2.1), relaxing the uniform Hölder stability to the pairwise fashion. The authors showed that sum-based models on multisets and graphs are lower-Hölder in expectation with different exponents, depending on the activation choices (Table 1). Motivated from the theoretical analysis, the authors proposed AdaptMPNN and SortMPNN, which can separate the adversarial graph pairs and demonstrated improved empirical performance to existing MPNNs on standard benchmark datasets.
weaknesses: 1. The proposed SortMPNN deserves more remarks and discussion; see details in the Question section.

2. Empirical Comparison between AdaptMPNN, SortMPNN versus ReluMPNN and SmoothMPNN: While the authors compare these MPNN variants in Table 2 with the synthetic $\epsilon$-Tree dataset, the authors do not include the results of ReluMPNN and SmoothMPNN on the standard benchmark datasets (e.g. TUDataset, LRGB). Including such comparison could better bridge the theoretical insights in the paper to practice.

Meta-review: This paper studies the separation power and stability of functions on graphs and multisets. In addition to the theoretical stability bounds, which hold in expectation, the paper provides explicit adversarial examples that illustrate the boundary cases of the theoretical results. It also empirically shows the connection between the Holder stability bounds and the model performance. The paper finds that common architectures have poor Holder exponents and separability properties and proposes two new message-passing neural networks with better separation properties All reviewers stated that this is a good paper.

Decision: Accept (Oral)