Deo, Murthy and Sarker

# COMBINING RETROSPECTIVE APPROXIMATION WITH IMPORTANCE SAMPLING FOR
OPTIMISING CONDITIONAL VALUE AT RISK

Anand Deo
Karthyek Murthy
Tirtho Sarker
Singapore University of Technology and Design
8 Somapah Rd  SINGAPORE 487372

## ABSTRACT

This paper investigates the use of retrospective approximation solution
paradigm in solving risk-averse optimization problems effectively via
importance sampling (IS). While IS serves as a prominent means for tackling
the large sample requirements in estimating tail risk measures such as
Conditional Value at Risk (CVaR), its use in optimization problems driven by
CVaR is complicated by the need to tailor the IS change of measure differently
to different optimization iterates and the circularity which arises as a
consequence. The proposed algorithm overcomes these challenges by employing a
univariate IS transformation offering uniform variance reduction in a
retrospective approximation procedure well-suited for tuning the IS parameter
choice. The resulting simulation based approximation scheme enjoys both the
computational efficiency bestowed by retrospective approximation and
logarithmically efficient variance reduction offered by importance sampling.

##  1 INTRODUCTION

Conditional value at risk (CVaR) serves as a widely used risk measure towards
assessing tail risks in quantitative risk management and operations research
(see mcneil2015quantitative,rockafellar2000optimization). For a loss
$\ell(\bm{X},\bm{\theta})$ associated with a decision choice $\bm{\theta}$
under a random realization $\bm{X},$ let $v_{\beta}{(\bm{\bm{\theta}})}$
denote the $(1-\beta)$-th quantile of $\ell(\bm{X},\bm{\theta}).$ Then its
CVaR at the tail-level $\beta\in(0,1)$ is given by,

$C_{\beta}(\bm{\theta})=E\left[\ell(\bm{X},\bm{\theta})\ |\
\ell(\bm{X},\bm{\theta})\geq v_{\beta}{(\bm{\bm{\theta}})}\right],$

which measures the average loss over the worst $\beta$-fraction of the
realizations. Under the assumption that $\ell(\bm{X},\cdot)$ is convex,
$C_{\beta}{(\bm{\bm{\theta}})}$ has been shown to possess favourable
properties such as convexity, subadditivity and coherence (see
acerbi2002coherence). Minimizing CVaR $C_{\beta}(\bm{\theta})$ over a compact
set $\Theta$ enjoys the following variational representation (see
rockafellar2000optimization),

$c_{\beta}=\inf_{u\in\mathbb{R},\bm{\theta}\in\Theta}\left[u+\beta^{-1}E\left(\ell(\bm{X},\bm{\theta})-u\right)^{+}\right]=\inf_{u\in\mathbb{R},\bm{\theta}\in\Theta}f(u,\bm{\theta}),$
(1)

which is a convenient starting point for solving optimization problems.
Consequently, CVaR has served as one of the primary vehicle for introducing
risk aversion in a variety of planning and resource allocation problems in
operations research and quantitative finance.

Since (1) is rarely solvable in closed form, Monte Carlo methods are often
deployed to obtain an approximate solution. Given $N$ i.i.d. samples of data
$\bm{X}_{1},\ldots,\bm{X}_{N}$ from the distribution of $\bm{X}$, define the
sample averaged objective,

$\hat{f}_{n}(u,\bm{\theta})=\left[u+\frac{1}{N\beta}\sum_{i=1}^{N}(\ell(\bm{X}_{i},\bm{\theta})-u)^{+}\right]$

Then the sample average approximation (SAA) to the optimisation problem (1)
may be constructed as

$\hat{c}_{n}=\inf_{u,\bm{\theta}}\hat{f}_{n}(u,\bm{\theta}).$ (2)

Desirable large sample properties such as consistency and asymptotic normality
are well-known for SAA estimators (see [Theorem 3.2]shapiro1991asymptotic) and
their limiting variances scale inversely with the tail level $\beta$ of
interest. This scaling is consistent with the understanding that one would
need approximately $\tilde{O}(\beta^{-1})$ to witness loss scenarios exceeding
the $(1-\beta)$-quantile $v_{\beta}(\bm{\theta}),$ as $\beta\rightarrow 0.$
The computational effort in solving the resulting optimization problems
becomes large as a consequence and it becomes necessary to resort to
specialized sampling schemes such as importance sampling (see, for example,
he2021adaptive and references therein) or aggregation sampling
(Fairbrother_2019).

In the evaluation of objective $C_{\beta}(\bm{\theta})$ at a fixed decision
choice $\bm{\theta},$ Importance Sampling (IS) is helpful if one can identify
an alternate distribution under which the rare excess loss event
$\\{\ell(\bm{X},\bm{\theta})\geq v_{\beta}(\bm{\theta})\\}$ occurs more
commonly. Indeed, IS has a rich literature on such change of measure
prescriptions offering efficient variance reduction in the estimation of rare
event probabilities (see, for example, Heidelberger,JUNEJA2006291, and more
recently, arief2021deep,deo2021achieving,ahn2021efficient). For instance with
a good change of measure suited for tail probabilities of
$\ell(\bm{X},\bm{\theta}),$ glynn1996importance,glasserman2000variance
demonstrate how the variance reduction offered by the change of measure can be
translated to efficient estimation of tail quantiles.
bardou2008computation,egloff2010quantile,he2021adaptive develop adaptive
algorithms which sequentially search for a good sampler choice within a IS
distribution family in the estimation of CVaR.

In contrast to the prolific use of IS in tail estimation, its use in solving
CVaR-driven optimization problems is limited by the availability of effective
IS distribution families suited for complex objectives which arise in planning
problems, and more severely, by the need to tailor the alternate sampling
distribution choice differently for different decision choices: A change of
measure which offers variance reduction for a decision choice
$\bm{\theta}\in\Theta$ may end up being inappropriate for a different decision
choice $\bm{\theta}^{\prime}\in\Theta$ (see Example 1 in Section 2 and Figure
1 for an illustration). In a setting where identifying a good change of
measure itself may require solving a non-trivial optimization problem (see,
for eg., bai2020rare), this dependence introduces a circular conceptual
difficulty: A sampler choice which attains a small limiting variance depends,
in turn, on the optimal decision choice which is the goal of our estimation to
begin with. See he2021adaptive for a detailed description of this circularity
issue and how adaptive approaches in lemaire2010unconstrained,he2021adaptive
can be useful in overcoming this challenge. The cross-entropy method by
RUBINSTEIN199789, which iteratively updates the IS distribution choice, by
minimizing Kullback-Liebler divergence between a proposed IS distribution and
the theoretically optimal IS distributions for increasingly rare instances,
remains the most widely adopted approach.

The effectiveness of these adaptive schemes rely on working with a IS
distribution family expressive enough to mirror the properties of
theoretically optimal samplers at different decision choices
$\bm{\theta}\in\Theta.$ Their use, on the contrary, becomes practical if the
chosen IS distribution family has a simple parametric representation (such as
an exponential family) and is easy to sample from. Verifiably efficient
prescriptions of IS distribution families have however remained elusive except
in elementary instances involving piece-wise linear objectives and elliptical
distributions modeling uncertainty in the problem.

In this paper, we propose a simulation based approximation scheme which embeds
importance sampling naturally in the well-known retrospective approximation
solution paradigm for solving optimization problems
(CHEN_2001,pasupathy2010choosing). The retrospective approximation framework
optimally balances the errors due to sampling and optimization approximations
(see pasupathy2010choosing) and lends itself naturally to the adaptive
selection of IS parameter choices.

For tackling the earlier highlighted challenges pertaining to the proposal and
sampling of IS distribution family, we employ IS transformations which are
structured sufficiently to induce distributions approximating the
theoretically optimal zero variance measures. This deviates from the
conventional approach of directly selecting a suitable IS distribution family,
which has proven to be model and distribution specific and face scalability
challenges. Specifically, we employ the single parameter family of self-
structuring transformations $\\{\bm{T}_{h}(\bm{X}):h>0\\}$ proposed in
deo2021achieving, where for every $h>0,$ the mapping
$\bm{T}_{h}:\mathbb{R}^{d}\to\mathbb{R}^{d}$ is a deterministic bijective
function. Inheriting variance reduction properties exhibited in the context of
tail probability estimation in deo2021achieving, we exhibit uniform variance
reductions offered by this IS transformation family in the evaluation of the
objective $f(u,\bm{\theta})$ in (1) over compact subsets of the decision set
$\Theta.$ This ability to achieve versatile variance reduction for a wide
range of optimization iterates renders the resulting IS estimators as a
natural choice towards untying the circularity issue highlighted earlier.

Our proposal to use retrospective approximation together with IS deviates from
the existing adaptive IS literature in which Robbins-Monro stochastic
approximation remains the preferred solution paradigm. The predominance of
stochastic approximation in adaptive IS can be understood from the observation
that samples from previous iterates are less suited to be used for the
gradient evaluation at the current iterate due the differing changes of
measure adopted in each iteration. Since our sampler is based on
transformations of the underlying vector $\bm{X}$ and enjoys robust variance
reduction properties, it frees us to explore computationally attractive
alternative solution paradigms. While our exposition treats the specific
sampler family to be introduced shortly, we note that the enhanced
retrospective approximation scheme developed in Section 4 can be used in
conjunction with other distribution families as well.

Notation: We use $\xrightarrow{{P}}$ to denote convergence in probability and
$\Rightarrow$ to denote convergence in distribution. Boldface letters denote
vectors. Likewise for a function $\bm{f}:\mathbb{R}^{d}\to\mathbb{R}^{k}$,
$\bm{f}(\bm{x})=(f_{1}(\bm{x}),\ldots,f_{k}(\bm{x}))$. We let
$N(\mu,\sigma^{2})$ denote a normal variable with mean $\mu$ and variance
$\sigma^{2}$. Let $\|\bm{x}\|_{p}$ denote the $\ell_{p}$ norm of a vector
$\bm{x}\in\mathbb{R}^{d}$ and $B_{r}(\bm{x})$ denote the $l_{\infty}$-metric
ball of radius $r$ centred at $\bm{x}$. For an increasing function
$f:\mathbb{R}\rightarrow\mathbb{R}$, we let $f^{-1}$ denote its left-inverse.
For real valued functions $f$ and $g$, we say that $f(x)=O(g(x))$ as
$x\to\infty$ if there exist positive constants $M,x_{0}$ such that for all
$x>x_{0}$, $|f(x)|\leq M|g(x)|$. We say that $f(x)=\tilde{O}(g(x))$ if
$f(x)=O(g(x)\log^{k}(x)),$ for some $k>0$.

##  2 VARIANCE REDUCTION WITH IS TRANSFORMATIONS

In this section we present an IS estimator for approximating CVaR objective.
Recall our proposal to induce a suitable IS distribution via a transformation
of the random vector $\bm{X}$. To accomplish this in our context, define the
$\mathbb{R}^{d}$-valued function
$\bm{T}_{h}(\bm{x}):=\bm{x}[s_{h}]^{\bm{\kappa}(\bm{x})}$, where
$s_{h}=h\log\log(1/\beta),$ with $h>0,$ is a stretch factor which stretches
the coordinates suitably via the exponent,

$\bm{\kappa}(\bm{x}):=\frac{\log(1+|\bm{x}|)}{\rho\|\log(1+|\bm{x}|)\|_{\infty}},$

in order to generate more samples in the extreme risk regions. The stretch
factor $s_{h},$ when viewed as a function of tail level $\beta,$ is increasing
when the problem is made rarer by letting $\beta$ smaller (we henceforth drop
the dependence on $\beta$ in our notation). Exponentiation is done component-
wise in the above expression for $\bm{T}_{h}(\bm{x})$ as in,
$\bm{T}_{h}(\bm{x})=(x_{1}s_{h}^{\kappa_{1}(\bm{x})},\ldots,x_{d}s_{h}^{\kappa_{d}(\bm{x})}).$
The map $\bm{T}_{h}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ can be shown to
be invertible almost everywhere on $\mathbb{R}^{d}$ (see [Proposition
1]deo2021achieving) and the resulting vector $\bm{Z}$ has a probability
density if $\bm{X}$ has a density. Letting $f_{\bm{X}}$ and $f_{\bm{Z}}$
denote the respective densities of $\bm{X}$ and $\bm{Z},$ the likelihood ratio
resulting from this change-of measure is given by,

$\displaystyle\mathcal{L}_{h}=\frac{f_{\bm{X}}(\bm{Z})}{f_{\bm{Z}}(\bm{Z})}=\frac{f_{\bm{X}}(\bm{Z})}{f_{\bm{X}}(\bm{X})}J_{h}(\bm{X})$
(3)

An explicit expression of the Jacobian,
$J_{h}(\bm{x})=\partial\bm{T}_{h}(\bm{x})/\partial\bm{x}$ in (3), given in (A)
in the Appendix, can be obtained by replacing $(u/l)$ in [Algorithm
1]deo2021achieving by $s_{h}.$ With this change-of-measure, we have the
following unbiased estimator for the objective function in (1):

$\hat{f}_{is,n}(u,\bm{\theta})=\left[u+\frac{1}{n\beta}\sum_{i=1}^{n}(\ell(\bm{Z}_{i},\bm{\theta})-u)^{+}\mathcal{L}_{h,i}\right],$
(4)

where ${\bm{X}}_{1},\ldots,{\bm{X}}_{n}$ are drawn i.i.d. from $\bm{X}$,
$\bm{Z}_{i}=\bm{T}_{h}(\bm{X}_{i})$ and $\mathcal{L}_{h,i}$ denotes the
likelihood (3) evaluated at $\bm{Z}_{i}$. Subsequently, one may optimise over
the IS weighted objective (4). The IS based estimator of minimum CVaR
therefore becomes

$\hat{c}_{is,n}=\inf_{u,\bm{\theta}}\hat{f}_{is,n}(u,\bm{\theta})$ (5)

###  2.1 A General Distribution Tail Model for Obtaining Variance Reduction
Guarantees

We now outline the tail modelling framework under which the IS scheme
described above provides substantial variance reduction. We say that
$f:\mathbb{R}\rightarrow\mathbb{R}$ is regularly varying if for all
$x\in\mathbb{R}_{+}$, $\lim_{t\rightarrow\infty}{f(tx)}/{f(t)}=x^{p}$, for
some $p\in\mathbb{R}$ (see [Definition B.1.1]de2007extreme). In this case, we
write $f\in\mathscr{RV}(p).$ We say that a function
$f:\mathbb{R}^{d}_{+}\rightarrow\mathbb{R}_{+}$ is multivariate regularly
varying if for any sequence $\bm{x}_{n}$ of $\mathbb{R}^{d}_{+}$ satisfying
$\bm{x}_{n}\rightarrow\bm{x}\neq\bm{0},$

$\displaystyle\lim_{n\rightarrow\infty}n^{-1}f(\bm{h}(n)\bm{x}_{n})=f^{\ast}(\bm{x}),$
(6)

for some limiting $f^{\ast}:\mathbb{R}^{d}_{+}\rightarrow(0,\infty)$ and a
component-wise increasing $\bm{h}(t)=(h_{1}(t),\ldots,h_{d}(t))$ satisfying
$h_{i}\in\mathscr{RV}(1/\rho_{i}),\ \rho_{i}>0,i=1,\ldots,d.$.

######  Assumption 1.

The marginal distribution of $\bm{X}=(X_{1},\ldots,X_{d})$ is such that the
hazard functions of $\\{X_{i}:i=1\ldots,d$, $\\{\Lambda_{i}:i=1.\ldots,d\\}$
are eventually strictly increasing and
$\Lambda_{i}\in\mathscr{RV}(\alpha_{i})$ for some $\alpha_{i}>0$. Further the
density of $\bm{X}$ when written in the form,

$\displaystyle f_{\bm{X}}(\bm{x})=\exp(-\psi(\bm{x}))\
\bm{x}\in\mathbb{R}_{d}^{+},\text{ satisfies
$\psi\in\textnormal{$\mathscr{MRV}$}(\psi^{\ast},\bm{h}).$}$ (7)

A wide variety of parametric and nonparametric multivariate distributions,
including normal, exponential family, elliptical, log-concave distributions
and Archimedian copula models satisfy Assumption 1. Marginal distributions
which satisfy $\Lambda_{i}\in\mathscr{RV}(\alpha_{i})$ include all
distributions that are either Weibull-type heavy-tailed or possess lighter
tails (such as exponential, normal, etc.).

######  Assumption 2.

There exists a limiting function $\ell^{*}(\cdot;\cdot)$ and $\rho>0$, such
that whenever $\bm{x}_{n}\to\bm{x}\neq\bm{0}$ and
$\bm{\theta}_{n}\to\bm{\theta}$,

$n^{-\rho}\ell(n\bm{x}_{n};\bm{\theta}_{n})\to\ell^{*}(\bm{x};\bm{\theta}).$
(8)

such that for every $\bm{\theta}$, $\ell^{*}(\cdot,\bm{\theta})$ is a
homogeneous function of $\bm{x}$.

Assumption 2 is satisfied for example, when
$\ell(\bm{x},\bm{\theta})=c(\bm{\theta}^{\intercal}\bm{x})$, such that for
$x\neq 0$, $c(tx)/t^{\rho}\to c^{*}(x)$ as $t\to\infty$. It is also satisfied
in more complicated examples, such as two-stage linear optimisation problems
where

$\ell(\bm{x},\bm{\theta})=\bm{c}^{\intercal}\bm{x}+\inf_{\bm{A}\bm{y}+\bm{b}=\bm{x}}\bm{\theta}^{\intercal}\bm{y}.$
(9)

We refer the reader to [Section 2 and Appendix B]deo2021achieving for a more
elaborate discussion on the conditions under which Assumptions 1 and 2 are
satisfied.

###  2.2 Uniform Variance Reduction Guarantees

One way to quantify the benefits gained by the use of IS transforms is to
analyse the variance of the functions $\hat{f}_{1}$ and $\hat{f}_{is,1}$ as
defined in (2) and (4). Note that these are simply the sample variances of the
objective functions in SAA and IS-SAA respectively. Proposition 1 below
quantifies these.

######  Proposition 1.

Suppose Assumption 1 holds. Define the set
$S_{r}=\\{(u,\bm{\theta}):u^{*}(1-r)\leq u\leq u^{*}(1+r),\
\theta_{i}^{*}(1-r)\leq\theta_{i}\leq\theta^{*}_{i}(1+r)\\}$. Then for any
$r,h_{\min},h_{\max}\in(0,\infty)$,

$\lim_{\beta\to 0}\sup_{(u,\bm{\theta})\in
S_{r},h\in(h_{\min},h_{\max})}\left|\frac{\log\text{var}\left[[(\ell(\bm{Z};\bm{\theta})-u)^{+}]\mathcal{L}_{h}\right]}{\log\text{var}([(\ell(\bm{X},\bm{\theta})-u)^{+}])}-2\right|=0\text{
and }$ (10)

$\sup_{(u,\bm{\theta})\in
S_{r},h\in(h_{\min},h_{\max})}\frac{\text{var}(\hat{f}_{is,1}(u,\bm{\theta}))}{\text{var}(\hat{f}_{1}(u,\bm{\theta}))}=o\left(\beta^{\frac{1}{1+\rho(r)}}\right),$
(11)

where $\rho(r)$ decreases to $0$ as $r\to 0$.

Proposition 1 quantifies the scale of variance reduction in the objective due
to the IS transformation by comparing it with the sample average objective.
The first part shows that the asymptotic variance reduction is optimal when
viewed in the logarithmic scale (see BJZWSC). Further, since $r$ appears
multiplicatively in $S_{r}$, equations (10) implies that a uniform reduction
in the variance of the objective function is obtained over a substantial
neighbourhood of the optimal point. On the other hand, (11) re-expresses this
improvement in performance in terms of the rarity parameter $\beta$.

######  Remark 1.

Solving (1) accurately requires accurate approximation of the objective
function $f$ by its IS-weighted sample average $\hat{f}_{is,n}$. To this end,
one seeks to find an alternate measure that resembles the conditional
distribution of the random vector $\bm{X}$ in region
$\\{\bm{x}:\ell(\bm{x},\bm{\theta})>u\\}$. Whenever $\bm{X}$ and the function
$\ell(\bm{x},\bm{\theta})$ satisfy Assumptions 1 and 2 respectively, the
distribution of $\bm{T}_{h}(\bm{X})$ resembles that of $\bm{X}$ in the region
$\\{\bm{x}:\ell(\bm{x};\bm{\theta})>u\\}$, irrespective of the choice of
decision parameter $\bm{\theta}$ (this can be seen through an application of
[Proposition 5.1]deo2021achieving). Proposition 1 showcases that this
translates to an exponential reduction in the variance error incurred while
estimating the objective function $f$ using the importance weighted sample
average $\hat{f}_{is,n}$.

######  Example 1.

To illustrate the difficulty in applying IS to the optimisation context, and
the implications of Proposition 1, consider a simple two dimensional setting,
where $\ell(\theta,x_{1},x_{2})=\theta x_{1}+(1-\theta)x_{2}$. Notice that the
region of concentration of the conditional distribution changes greatly when
$\theta$ is changed from $0.2$ to $0.8$. Therefore, an IS distribution which
provides a large variance reduction in the objective for $\theta=0.2$ may not
for $\theta=0.8$. However, the uniformity of variance reduction in (10)
demonstrates that by use of IS transformations this difficulty may be
mitigated. For the same linear loss but with $\bm{X}\in\mathbb{R}^{5}$, Figure
1(b) shows that the worst case standard error over $S_{r}$ is much smaller for
IS than it is for SAA for $\beta=0.001$. This demonstrates that significant
benefits are obtained through use of IS transformations for realistic values
of $\beta$ (rather than the theorem truly holding only asymptotically).

Proposition 1 can be used to derive variance guarantees on the error in
optimising CVaR with (4).

###### Theorem 1

Under Assumptions 1 and 2, we have
$\sqrt{n}(\hat{c}_{is,n}-c_{\beta})\Rightarrow\sigma_{is}(\beta)N(0,1)$ where

$\sigma_{is}^{2}(\beta)=\text{var}\left((\ell(\bm{Z};\bm{\theta}^{*})-v_{\beta}(\bm{\theta}^{*}))^{+}\mathcal{L}_{h}\right)$
(12)

Further, the limiting variance $\sigma^{2}_{is}(\beta)$ of the IS estimator is
smaller than the limiting variance $\sigma_{c}^{2}(\beta)$ of the SAA
estimator (1) as given by the relationship,

$\frac{\sigma^{2}_{is}(\beta)}{\sigma^{2}_{c}(\beta)}=o(\beta^{1-\varepsilon}),\text{
for any }\varepsilon>0,$

Considering the proposed change of measure for CVaR optimisation, Theorem 1
guarantees a sample complexity of $o(\beta^{-\varepsilon})$ as $\beta\searrow
0,$ where $\varepsilon>0$ can be made arbitrarily small. With the variance
reduction guarantee holding for any choice of hyper-parameter $h>0,$ an
effective $h$ can be chosen via cross-validation without incurring a change of
scaling in sample complexity as demonstrated in Section 4.

(a) Concentration of excess distribution with $\theta=0.2,0.8$

(b) Worst case standard errors over $S_{r}$

Figure 1: Figures 1(a) plots the conditional distribution of $\theta
X_{1}+(1-\theta)X_{2}$ conditional on it exceeding $3.5$. Figure 1(b) shows
the worst case relative errors (ratio of standard deviation to mean) for the
functions $f_{is}(\cdot)$ over $S_{r}$ and compares them with SAA.

##  3 A VANILLA RETROSPECTIVE APPROXIMATION SCHEME FOR CVAR OPTIMIZATION

The procedure described in Section 2 arrives at the optimal CVaR by generating
a single sample path $\\{\bm{X}_{1},\ldots,\bm{X}_{n}\\}$ and then solving the
resulting IS-weighted SAA problem. However, in practical implementation,
solutions are not available in closed form, and sample path problems such as
(4) need to be solved using deterministic optimisation algorithms. Thus,
obtaining a solution to within a small tolerance level may be a
computationally challenging task. With this in mind, sequential procedures are
often deployed to solve SAA (see
bayraksan2011sequential,pasupathy2010choosing). One such method is
Retrospective Approximation (RA), which generates a sequence of SAA sub-
problems, and progressively increases sample sizes, while solving the
resulting SAA sub-problems with progressively reducing error tolerances.
Further, the solution of each stage is fed into as the starting point for the
optimiser in the subsequent stage. This eases the overall computational burden
as follows: the initial stages obtain a rough estimate of the optimal
parameters while not expending too much computation. The later stages are
accurate, since the initial solution to optimisation algorithm used to solve
the SAA sub-problem is close to the true solution.

In our context, given a sequence of sample sizes $\\{m_{k}:k\geq 1\\}$, this
amounts to minimising a sequence of random functions
$\\{f_{is,m_{k}}(u,\bm{\theta}):k\geq 1\\}$, where $f_{is,n}(u,\bm{\theta})\to
f(u,\bm{\theta})$ in probability as $n\to\infty$. Further, at each stage $k$,
we set the initial solution to the $(k+1)th$ sub-problem to be
$(u_{k},\bm{\theta}_{k})$ the minimiser of $f_{m_{k}}$. Algorithm 1 gives a RA
based implementation of the IS-weighted SAA algorithm from Section 2. Denote
the solution of the sample path problem in (16) by
$(u_{k}^{*},\bm{\theta}_{k}^{*})$. Notice that this is the true solution of
the optimisation problem and not the one computed to $\varepsilon_{k}$
precision. Our first result shows that use of IS significantly reduces sample
errors in the solution to the sequence of sub-problems. Define

$\bm{g}(u;\bm{\theta})=\left[1-\beta^{-1}{E}\left(\bm{I}(\ell(\bm{X},\bm{\theta})\geq
u)\right),\beta^{-1}E\left(\nabla_{\bm{\theta}}\ell(\bm{X},\bm{\theta})\bm{I}(\ell(\bm{X},\bm{\theta})\geq
u)\right)\right]\text{ and }$ (13)
$\mathbf{G}(\bm{x};u,\bm{\theta},r)=\left[1-\beta^{-1}\bm{I}(\ell(\bm{x},\bm{\theta})\geq
u)\mathcal{L}_{r}(\bm{x})\ ,\
\beta^{-1}\left(\nabla_{\bm{\theta}}\ell(\bm{x},\bm{\theta})\bm{I}(\ell(\bm{x},\bm{\theta})\geq
u)\right)\mathcal{L}_{r}(\bm{x})\right],$ (14)

where, $\nabla_{\bm{\theta}}$ denotes the derivative with respect to
$\bm{\theta}$.

######  Proposition 2.

Suppose Assumptions 1 and 2 are satisfied. Further, let the problem in (1)
have a unique minimiser. Then, as $k\to\infty$, we have

$\sqrt{m_{k}}[(u_{k}^{*},\bm{\theta}_{k}^{*})-(v_{\beta}^{*}(\bm{\theta}),\bm{\theta}^{*})]\xrightarrow{L}N(\bm{0},\Sigma_{h}),\text{
where
}\Sigma_{h}=[\nabla\bm{g}(u^{*};\bm{\theta}^{*})]^{-\intercal}\text{var}\left(\bm{G}(\bm{Z};u^{*},\bm{\theta}^{*},h)\right)[\nabla\bm{g}(u^{*};\bm{\theta}^{*})]^{-1}$
(15)

Further, the Frobenius norm of the limiting covariance satisfies
$\|\Sigma_{h}\|_{F}=o(\beta^{-\varepsilon})$ for any $\varepsilon>0$.

Proposition 2 showcases the utility of using IS in conjunction with RA to get
a better performance in terms of estimation errors. In particular, the second
part above suggests that the sample complexity of estimation grows as
$o(\beta^{-\varepsilon})$ rather than $\tilde{O}(\beta^{-1})$.

Input: Target tail probability level $\beta$, initial solution
$u_{0},\bm{\theta}_{0}$, sample sizes $m_{1},\ldots,m_{k}$, error tolerances
$\varepsilon_{1},\ldots,\varepsilon_{k}$.

For $k\geq 1$, do

1\. Transform the samples: For each sample $i=1,\ldots,m_{k},$ compute the
transformation,
$\displaystyle\bm{Z}_{i}=\bm{T}_{h}(\bm{X}_{i}):=\bm{X}_{i}[s_{h}]^{\kappa(\bm{X}_{i})},$
2\. Solve the IS based optimisation:
$\hat{c}_{is,m_{k}}:=\inf_{u,\bm{\theta}}\left[u+\frac{1}{m_{k}\beta}\sum_{i=1}^{m_{k}}\big{(}\ell({\bm{Z}}_{i},\bm{\theta})-u\big{)}^{+}\mathcal{L}_{h,i},\right]\text{
to a tolerance of }\varepsilon_{k}$ (16) with an initial solution
$(u_{k-1},\bm{\theta}_{k-1})$. Return also the optimiser of (16),
$(u_{k},\bm{\theta}_{k})$.

Algorithm 1 Retrospective Approximation based CVaR Optimisation –
$\textsc{OPT}(h)$

Notice that Proposition 2 gives no mention of the computation required to
solve the optimisation problem in Algorithm 1, and instead focuses on the
quality of the optimal solutions to the sample path sub-problems. In the next
discussion, we argue the computational benefits obtained by the use of
retrospective approximation continue to hold even when IS is used. Assumption
3 below imposes a mild condition on $m_{k},\varepsilon_{k}$ (see [Assumptions
C.1-C.3]pasupathy2010choosing).

######  Assumption 3.

Suppose that the sequence $\\{(\varepsilon_{k},m_{k}):k\geq 1\\}$ satisfy the
following requirements:

  1. 1.

If the optimisation procedure used to solve the individual sample path
problems exhibits linear convergence,
$\liminf_{k\to\infty}\varepsilon_{k-1}\sqrt{m}>0$. If this procedure exhibits
polynomial convergence, $\liminf_{k\to\infty}{\log
1/\sqrt{m_{k-1}}}({\log\varepsilon_{k}})^{-1}>0$.

  2. 2.

$\limsup_{k\to\infty}(\sum_{j=1}^{k}m_{j})^{2}/\varepsilon_{k}^{2}<\infty$

  3. 3.

$\limsup_{k\to\infty}m_{k}^{-1}\sum_{j=1}^{k}m_{j}<\infty$.

Observe that there are two competing errors in the solution output by
Algorithm 1. First, there is the sample error, which can be quantified by the
limit theorem in Proposition 2. Second, there is the error due to solving the
optimisation problem imperfectly (up to $\varepsilon_{k}$ accuracy).
Assumption 3 imposes conditions so that these errors are balanced out and an
optimal rate of convergence of the solution is obtained. Notice that the total
work done in running Algorithm 1 for $k$ epochs equals $W_{k}=\sum_{j\leq
k}N_{j}m_{j}$, where $N_{j}$ is the (random) number of calls made to the
deterministic algorithm used to solve (16) to $\varepsilon_{j}$ accuracy.
Proposition 3 establishes that the work normalised error for the RA based IS-
weighted CVaR optimisation remains bounded.

######  Proposition 3.

Under Assumptions 1 \- 3, we have
$W_{k}[\sigma_{is}(\beta)]^{-2}\|(u_{k},\bm{\theta}_{k})-(v_{\beta}(\bm{\theta}^{*}),\bm{\theta}^{*})\|_{2}^{2}=O_{p}(1)$.

Notice the appearance of the limiting variance $\sigma_{is}(\beta)$ from (3)
above. This suggests that the IS scheme continues to enjoy a significant
reduction over SAA in work normalised errors. Indeed this is showcased in
experiments in Section 5, where we show that the amount of computing effort
required to solve the CVaR optimisation problem to a fixed accuracy is
substantially less when IS is used.

##  4 ENHANCED RETROSPECTIVE APPROXIMATION WITH ADAPTIVITY IN IS CHOICE

One major shortcoming of Algorithm 1 is the lack of optimisation over the
cross-validation parameter $h$. In this discussion, we develop an enhanced
version of the retrospective approximation based CVaR optimisation which
includes a subroutine to iteratively arrive at a value of $h$ which gives a
lower error. To simplify matters, we assume the existence of a noisy oracle,
which given $n$ samples of data and $(u,\bm{\theta})$ returns a value
$\hat{h}_{n}(u,\bm{\theta})$, such that
$\hat{h}_{n}(u_{n},\bm{\theta}_{n})\xrightarrow{P}h(u,\bm{\theta})$ whenever
$(u_{n},\bm{\theta}_{n})\to(u,\bm{\theta})$ as $n\to\infty$. We further assume
that $h(u^{*},\bm{\theta}^{*})$ is a good hyper-parameter choice for problem
(the specific implementation for $h_{n}$ and $h$ is presented in Section 5).
Under this assumption, Algorithm 2 presents an enhancement to the RA based IS-
weighted CVaR optimisation scheme.

Input: Target tail probability level $\beta$, samples $\bm{X}_{1},\ldots,$
from $f_{\bm{X}}(\cdot)$, initial seeds $u_{0},\bm{\theta}_{0},h_{0}$.

For $k\geq 1$, do

1\. IS-Weighted CVaR optimisation: With a sample size of $m_{k,1}$ and error
tolerance $\varepsilon_{k,1}$, implement steps 1 and 2 of Algorithm 1 starting
from $(u_{k-1},\bm{\theta}_{k-1})$ and with $h=h_{k-1}$.

2 Update the cross validation parameter:  Using a sample size of $m_{k,2}$,
return the updated value $h_{k}=\hat{h}_{m_{k,2}}(u_{k},\bm{\theta}_{k})$

Algorithm 2 Enhanced Retrospective Approximation based CVaR Optimisation

It possesses the following two benefits over prominent adaptive IS schemes:

For every value of $\\{h(u,\bm{\theta}):u>0,\bm{\theta}\in\Theta\\}$, the
distribution for $\bm{Z}$ may be generated easily by suitably transforming the
samples generated from the distribution of $\bm{X}$. The ability to cheaply
obtain samples from an alternate distribution, as assumed with existing IS
techniques, could be an issue when the IS distribution family need to be
expressive enough to approximate the zero variance measures for complex
optimization objectives.

Similarly, given a value of $(u,\bm{\theta})$, Algorithm 2 deals with the
selection of $h(u,\bm{\theta})$ directly from samples of $\bm{X}$, rather than
assuming a black box access to a good choice of $h(u,\bm{\theta}).$ For a
fixed choice of $u,\bm{\theta},$ the selection of $h(u,\bm{\theta})$ is
fortunately a one-dimensional problem which can be tackled effectively with
bisection or grid search and using common random numbers. Using the same
samples of $\bm{X}$ to reduce the variance in IS parameter selection is not
possible with most existing IS approaches. As will be explained shortly, the
computed values of the hyper-parameter will converge in probability to
$h(u^{*},\bm{\theta}^{*})$. Since the entire cross validation procedure only
requires access to samples of $\bm{X}$, it is more easily implementable than
adaptive techniques which require black box access to a good importance
sampler at each stage.

Observe that given the $(k-1)$th stage solution,
$(h_{k-1},u_{k-1},\bm{\theta}_{k-1})$ and $m_{k,1}$ number of samples, as a
consequence of the discussion in [Theorem 2]pasupathy2010choosing,
$(u_{k},\bm{\theta}_{k})=(u^{*},\bm{\theta}^{*})+O_{p}(m_{k,1}^{-1/2}+\varepsilon_{k})$.
Further, notice that with $m_{k,2}$ samples used for cross-validation in the
first stage, one has
$\hat{h}_{m_{k,2}}(u_{k-1},\bm{\theta}_{k-1})=h(u_{k},\bm{\theta}_{k})+o_{{p}}(1)$.
From the continuity of $h(\cdot)$ in the parameter, and the consistency of
$(u_{k-1},\bm{\theta}_{k-1})$ as $k\to\infty$, we further have that
$\hat{h}_{m_{k,2}}(u_{k-1},\bm{\theta}_{k-1})=h(u^{*},\bm{\theta}^{*})+o_{{p}}(1)$.
The above discussion suggests that Algorithm 2 arrives at a statistically
consistent estimate of the desired optimal hyper-parameter. It is of note to
observe that the convergence to the optimal $h(u^{*},\bm{\theta}^{*})$ appears
to be insensitive to the choice of $h_{0}$. Indeed, this is demonstrated
experimentally in Section 5.

A natural question to ask given the above discussion, is whether one may
derive a limit theorem which states that as $k\to\infty$,
$\sqrt{m}_{k}(\hat{c}_{is,m_{k}}-c_{\beta})\xrightarrow{L}N(0,[\sigma_{is}^{*}(\beta)]^{2})$,
where $[\sigma_{is}^{*}(\beta)]^{2}$ is the variance of the objective
evaluated at $(u^{*},\bm{\theta}^{*})$ under the distribution
$\bm{Z}_{h(u^{*},\bm{\theta}^{*})}$. Addressing this is an interesting follow-
up direction. We next attempt to justify the benefits obtained by the use of
Algorithm 2 through numerical experiments.

##  5 NUMERICAL EXPERIMENTS

###  5.1 Minimization of CVaR for Linear Portfolios

In order to demonstrate the efficacy of our algorithms in practice, we
consider solving the constrained minimum CVaR minimisation portfolio
optimisation problem. In this setting,
$\ell(\bm{x},\bm{\theta})=\bm{\theta}^{\intercal}\bm{x}$, and the set
$\Theta=\\{\bm{\theta}:\bm{\theta}^{\intercal}\bm{1}=1\\}$. The marginals of
$\bm{X}$ are taken to have the c.d.f.s $F_{i}(x)=P(X_{i}\leq
x)=1-\mathrm{e}^{-x^{\alpha_{i}}}$ where $\alpha_{i}=0.5\ \forall i$.
Dependence is modelled through a Gaussian copula whose covariance matrix
$\bm{R}$ is designed to capture a realistic degree of correlation among
various asset returns. We compare the effectiveness of the estimators returned
by Algorithms 1 and 2 with plain SAA using three performance measures:

First, we use relative mean square error as a metric to compare the quality of
the optimal value output by each of the algorithms. That is, given the output
of an algorithm, call it $\hat{c}_{\beta}$, we compute the relative root mean
square error as
$\sqrt{\frac{1}{K}\sum_{i=1}^{K}({\hat{c}_{\beta,i}}/{c_{\beta}}-1)^{2}}$
where $\hat{c}_{\beta,i}$ are outputs of independent runs of the algorithm and
$c_{\beta}$ is defined in (1). We use relative errors as metrics of
performance both here as well in the next experiment to facilitate a scale
free (i.e: $\beta$ independent) comparison between different algorithms. In
order to compute the optimal value $c_{\beta}$, we solve the SAA problem (2)
with $N=10^{6}$ samples. In each case, the value of $\beta$ is varied from
$\beta=0.037$ to $\beta=5\times 10^{-4}$. In order to ensure a fair comparison
between each of the three methods, we keep the total sample budget fixed and
equal to $2500$. For the implementation in Algorithm 2, this budget is divided
among two epochs as $m_{1}=500$ and $m_{2}=2000$. Further, for every
$(u,\bm{\theta})$, we let
$\hat{h}_{n}(u,\bm{\theta})=\inf_{h}\frac{1}{n}\sum_{i=1}^{n}\mathbf{I}(\ell(\bm{Z}_{h,i};\bm{\theta})\geq
u)\mathcal{L}_{h,i}^{2}$. Observe that from the law of large numbers this
converges to the second moment of the gradient of the objective function with
respect to $u$. Since the final errors in decision (see Proposition 2 for
example) depend on this quantity having a smaller variance,
$\hat{h}_{n}(u,\bm{\theta})$ as chosen above is reasonable. For Algorithm 1
for all $\beta$, we choose $h=2.5$. Figure 2(a) details the results.

(a) Relative RMSE in CVaR optimisation

(b) Relative regret in CVaR optimisation

Figure 2: Figure 2(a) compares the relative RMSE in CVaR optimisation using
Algorithms 1 and 2 with plain SAA. Figure 2(b) compares the corresponding out-
of-sample regret. In each of the figures, cross-marks respectively denote
estimated quantities.

We next use relative regret as a metric to compare out-of-sample performance
each of the algorithms. Given the output of an algorithm $\hat{\bm{\theta}}$,
we define the relative regret as
$\hat{r}_{\beta}(\hat{\bm{\theta}})=\frac{1}{K}\sum_{i=1}^{K}({c_{\beta}(\hat{\bm{\theta}}_{i})}/{c_{\beta}}-1)\times
100\%$, where $\hat{\bm{\theta}}_{i}$ are outputs of independent runs of the
algorithm. The first observation is that IS-RA significantly outperforms SAA
in terms of relative errors obtained. Secondly, adaptively optimising over $h$
as in Algorithm 2 leads to a further improvement over Algorithm 1. Finally,
figure 3(a) displays the relative regret incurred as a function of the
starting point of the $h_{0}$. Notice that there appears to be a degree of
robustness in the performance of the algorithm to the specific value of
$h_{0}$ selected.

Finally, in order to compare in terms of the effort required to obtain a
desired out of sample accuracy, we compute the number of samples required by
each method to obtain $1\%$ relative regret; refer Figure3(b). Observe that
for $\beta=0.037$, for IS, this is $\approx 600$, while SAA requires $\approx
5500$ samples. This difference is even more pronounced when $\beta=0.003$,
where SAA requires roughly $28000$ samples, while IS only requires $1175$.

(a) Relative regret as a function of $h_{0}$

(b) Samples required to obtain $1\%$ relative regret

Figure 3: Figure 3(a) plots the relative regret incurred in optimisation as a
function of the initial value $h_{0}$. Figure 3(b) plots the number of samples
required by each method to obtain $1\%$ relative regret. In each of the
figures, cross-marks respectively denote estimated quantities.

###  5.2 A Portfolio Credit Risk Example

Consider the problem of selecting a portfolio of loans from among $K$ classes
(industries). Given a that a set of common market factors $\bm{X}$ realise a
value $\bm{x}$, a loan belonging to class $i$ defaults with probability
$p_{i}(\bm{x})$ independently of everything else. We further suppose that the
loss given default for loan $k$ is given by $e_{k}$ (random), independent of
the market variable. Denoting $n_{i}$ as the number of loans belonging to
class $i$, the total loss given a portfolio $\bm{\theta}$ is
$\ell(\bm{\theta},\bm{X})=\sum_{i=1}^{K}\theta_{i}(\sum_{k=1}^{n_{k}}e_{k}\mathbf{I}(\text{loan
k defaults}))$. The objective is to find the portfolio with the minimum CVaR,
subject to the total returns from the investment exceeding a nominal
threshold, that is $\bm{\theta}^{\intercal}\bm{r}\geq q$ where the portfolio
returns are given by $\bm{r}$. In this example, $K=2$ and
$\bm{X}\in\mathbb{R}^{4}$ has Weibull marginals with a Gaussian copula used to
model joint dependence. We further assume that $p_{i}(\bm{x})$ has a logistic
form. Importance sampling is deployed only on the common market factors
$\bm{X}$, and the IS weighted algorithm is solved as before. In this
implementation, $N=2000$ while $n_{1}=n_{2}=5000$ and $e_{i}$ are uniform
random variables. Figure 4 displays the results of the experiment. Notice that
the use of IS on the common market variable significantly improves lowers
errors in optimal VaR estimation. The relative error in CVaR estimation is
roughly $1\%$ when IS is used as opposed to over $50\%$ without the use of IS.

Figure 4: Error in CVaR Computation for the portfolio credit risk problem

##  6 OUTLINES OF KEY PROOF STEPS

Proof of Proposition 1: For notational convenience, let $M_{2,u}$ denote the
second moment of $[\mathcal{L}_{h}(\ell(\bm{Z},\bm{\theta})-u)^{+}]$. Note
that this is an upper bound for the variance which we want to compute.
Further, define $t(u)=\Lambda_{\min}(u^{1/\rho})$, $q_{\infty}(u)$ as the
inverse of $\Lambda_{\min}(u)$ and
$\bm{Y}_{u}=[t(u)]^{-1}\bm{\Lambda}(\bm{X})$. Changing variables from $\bm{X}$
to $\bm{Y}_{u}$ in the expectation below (see (EC.16) onward in the proof of
Lemma EC.6 of deo2021achieving for detailed steps in a similar change of
variables exercise), we obtain
$M_{2,u}=\mathbb{E}\left[(\ell(\bm{X},\bm{\theta})-u)^{2}\frac{f_{\bm{X}}(\bm{X})}{{f}_{\bm{X}}(\bm{T}^{-1}_{h}(\bm{X}))}J(\bm{T}^{-1}_{h}(\bm{X}));\ell(\bm{X},\bm{\theta})\geq
u\right]$. This simplifies to
$\mathbb{E}\left[\exp(-t(u)F_{u}(\bm{Y}_{u}))\right]$ where
$F_{u}(\bm{p})=a_{u}(\bm{p})+b_{u}(\bm{p})+c_{u}(\bm{p})-2d[t(u)]^{-1}\log
t(u)+\chi_{\textnormal{lev}_{1}^{+}(\ell_{u,\bm{\theta}})}(\bm{p})$ where
$a_{u}(\bm{p})=[t(u)]^{-1}\left[\log f_{\bm{Y}}(\bm{\psi}_{u}(t(u)\bm{p})-\log
f_{\bm{Y}}(t(u)\bm{p})\right]$,
$b_{u}(\bm{p})=[t(u)]^{-1}\left[\sum_{i=1}^{d}\left[\log\lambda_{i}(\bm{T}^{-1}_{h,i}(\bm{q}(t(u)\bm{p})))-\log\lambda_{i}(q_{i}(t(u)p_{i}))\right]-\log
J_{h}(\bm{T}^{-1}(\bm{q}(t(u)\bm{p})))\right]$ and
$c_{u}(\bm{p})=-2[t(u)]^{-1}\log(\ell(\bm{q}(t(u)\bm{p},\bm{\theta})-u))$.
Define $p_{u}=P(\ell(\bm{X},\bm{\theta})\geq u)$. Then an application of the
general Varadhan’s integral lemma [Varadhan (1988), Theorem 2.1], for all
$(u,\bm{\theta})$, we have that $M_{2,u}=o(p_{u}^{2-o(1)})$ (see the proof of
[Theorem 1]deo2021efficient). Similarly, we have that
$\text{var}[(\ell(\bm{X},\bm{\theta})-u)\mathbf{I}(\ell(\bm{X},\bm{\theta})\geq
u)]=\Theta(p_{u})$. Combining these two proves the claim for a fixed
$(u,\bm{\theta})$. Notice that $\inf\\{u:u\in
S_{r}\\}\geq\log^{\varepsilon}1/\beta$ for all $\beta$ small enough. Further,
$\inf_{h>h_{\min}}s_{h}=o(-\log^{\varepsilon}\beta)$. The uniformity over
$h,u,\bm{\theta}$ is now obtained using the uniform convergence of $\ell$ as
in Assumption 2 and then following the proof of [Lemma D.2]deo2021achieving.
For part 2, notice that from the above discussion, for any $(u,\bm{\theta})$,
the ratio of the variances of $\hat{f}_{is,1}$ and $\hat{f}_{1}$ is given by
$o(p_{u}^{1-\varepsilon})$. Under Assumptions 1 and 2, it can be established
that for some continuous function $a(\cdot)$,
$C_{\beta}(\bm{\theta})=a(\bm{\theta})(1+o(1))q_{\infty}^{\rho}(-\log\beta)$.
As a consequence of [Theorem 7.33]rockafellar2009variational,
$\bm{\theta}^{*}$ which minimises $C_{\beta}(\bm{\theta})$ satisfies
$\bm{\theta}^{*}=\bm{\theta}_{1}+o(1)$, where $\bm{\theta}_{1}\in\arg\min
a(\bm{\theta})$. Then, the optimal
$u^{*}=v_{\beta}(\bm{\theta}^{*})=I^{*}q_{\infty}^{\rho}(-\log\beta)[1+o(1)]$,
where $I^{*}=\inf_{\bm{\theta}}a(\bm{\theta})$. This implies that $u\in
S_{r}$, $u\geq(1-r)I^{*}q_{\infty}^{\rho}(-\log\beta)$. Finally, it can be
seen that $\log p_{u}=\Lambda_{\min}(u/a(\bm{\theta}))(1+o(1))$ (apply
[Theorem 4.2]deo2021achieving). Plugging everything together, we have that
$\log p_{u}\leq(1-r)^{\alpha_{\min}}(I^{*}/a(\bm{\theta}))\log\beta$. Now, the
continuity of $a(\cdot)$ implies that whenever $(u,\bm{\theta})\in S_{r}$, and
$h>h_{\min}$, $p_{u}\leq\beta^{1/1-\rho(r)}$, for $\rho(r)\to 0$. ∎

Proof of Theorem 1:  Notice first that the limiting function in the
optimisation, $f$ may be written as the expected value of $\hat{f}_{is,1}$.
Further, observe that due to Assumptions 1 and 2 , $\hat{f}_{is,1}$ has a
finite variance. Therefore, by an application [Shapiro (1991), Theorem 3.2],
one obtains the first part of Theorem 1. The quantification of the variance
reduction is obtained following the proof of Proposition 1, with
$p_{v_{\beta(\bm{\theta}^{*})}}=\beta$.∎

Proof of Proposition 2 It is sufficient to verify the conditions from [Theorem
4]pasupathy2010choosing. First, note that the optimal solutions
$(u_{k}^{*},\bm{\theta}_{k}^{*})$ are zeros
$\nabla\hat{f}_{is,m_{k}}(u,\bm{\theta})$. Further, we have that
$E\bm{G}(u,\bm{\theta})=g(u,\bm{\theta})$ for all $(u,\bm{\theta})$. Finally,
owing to Assumptions 1 and 2, the usual conditions required for a CLT hold,
and we therefore have
$\sqrt{m_{k}}(\nabla\hat{f}_{is,m_{k}}(u,\bm{\theta})-\bm{g}(u,\bm{\theta}))\to
N(0,\text{var}(\bm{G}(\bm{Z};u,\bm{\theta};h)))$. By the assumptions of the
theorem, there exists a unique minimiser, and the gradients of $\bm{g}$ is
non-singular. Finally, the uniform convergence assumption holds as a result of
the discussion following the statement of [Theorem 4]pasupathy2010choosing. ∎

Proof of Proposition 3 This follows from Assumption 3, Proposition 2 and the
proof of [Theorem 5]pasupathy2010choosing.∎

##  Appendix A Expression for the Jacobian $J_{h}$

Recall that $J_{h}(\cdot)$ is the Jacobian of the transformation
$\bm{T}_{h}(\cdot)$. This is given by

$J_{h}(\bm{x}):=\left[\prod_{i=1}^{d}\tilde{J}_{i}(\bm{x})\right]\times\frac{s_{h}^{\bm{1}^{\intercal}\bm{\kappa}(\bm{x})}}{\max_{i=1,\ldots,d}\tilde{J}_{i}(\bm{x})}\text{
where
}\tilde{J}_{i}(\bm{x}):=1+\frac{\rho^{-1}\log(s_{h})}{\|\log(1+|\bm{x}|)\|_{\infty}}\frac{|x_{i}|}{1+|x_{i}|},\quad
i=1,\ldots,d.$

## ACKNOWLEDGEMENTS

The authors acknowledge support from Singapore Ministry of Education grant
MOE2019-T2-2-163.

## References

  * Acerbi and Tasche (2002) Acerbi, C., and D. Tasche. 2002.  “On the coherence of expected shortfall”.  Journal of Banking & Finance 26(7):1487–1503.
  * Ahn and Zheng (2021) Ahn, D., and L. Zheng. 2021.  “Efficient simulation for linear programming under uncertainty”.  In 2021 Winter Simulation Conference (WSC), 1–12.  IEEE.
  * Arief et al. (2021) Arief, M., Z. Huang, G. K. S. Kumar, Y. Bai, S. He, W. Ding, H. Lam, and D. Zhao. 2021.  “Deep Probabilistic Accelerated Evaluation: A Robust Certifiable Rare-Event Simulation Methodology for Black-Box Safety-Critical Systems”.  In International Conference on Artificial Intelligence and Statistics, 595–603.  Proceedings of Machine Learning Research.
  * Bai et al. (2022) Bai, Y., Z. Huang, H. Lam, and D. Zhao. 2022, feb.  “Rare-Event Simulation for Neural Network and Random Forest Predictors”.  ACM Trans. Model. Comput. Simul..  Just Accepted.
  * Bardou et al. (2009) Bardou, O., N. Frikha, and G. Pagès. 2009.  “Computing VaR and CVaR Using Stochastic Approximation and Adaptive Unconstrained Importance Sampling”.  Monte Carlo Methods and Applications 15(3):173–210.
  * Bassamboo et al. (2005) Bassamboo, A., S. Juneja, and A. Zeevi. 2005.  “Expected Shortfall in Credit Portfolios with Extremal Dependence”.  In Proceedings of the 2005 Winter Simulation Conference, edited by M. E. Kuhl, N. M. Steiger, F. B. Armstrong, and J. A. Joines, 1849 – 1858\.  Piscataway, New Jersey: Institute of Electrical and Electronics Engineers, Inc.
  * Bayraksan and Morton (2011) Bayraksan, G., and D. P. Morton. 2011.  “A sequential sampling procedure for stochastic programming”.  Operations Research 59(4):898–913.
  * Chen and Schmeiser (2001) Chen, H., and B. W. Schmeiser. 2001, mar.  “Stochastic root finding via retrospective approximation”.  IIE Transactions 33(3):259–275.
  * de Haan and Ferreira (2007) de Haan, L., and A. Ferreira. 2007.  Extreme Value Theory: An Introduction.  Springer Science & Business Media.
  * Deo and Murthy (2021a) Deo, A., and K. Murthy. 2021a.  “Achieving Efficiency in Black Box Simulation of Distribution Tails with Self-structuring Importance Samplers”.  arXiv preprint arXiv:2102.07060.
  * Deo and Murthy (2021b) Deo, A., and K. Murthy. 2021b.  “Efficient Black-Box Importance Sampling for VaR and CVaR Estimation”.  arXiv preprint arXiv:2106.10236.
  * Egloff and Leippold (2010) Egloff, D., and M. Leippold. 2010.  “Quantile Estimation with Adaptive Importance Sampling”.  The Annals of Statistics 38(2):1244–1278.
  * Fairbrother et al. (2019) Fairbrother, J., A. Turner, and S. W. Wallace. 2019, nov.  “Problem-driven scenario generation: an analytical approach for stochastic programs with tail risk measure”.  Mathematical Programming 191(1):141–182.
  * Glasserman et al. (2000) Glasserman, P., P. Heidelberger, and P. Shahabuddin. 2000.  “Variance Reduction Techniques for Estimating Value-at-Risk”.  Management Science 46(10):1349–1364.
  * Glynn (1996) Glynn, P. W. 1996.  “Importance Sampling for Monte Carlo Estimation of Quantiles”.  In Mathematical Methods in Stochastic Simulation and Experimental Design: Proceedings of the 2nd St. Petersburg Workshop on Simulation, 180–185.  Publishing House of St. Petersburg University.
  * He et al. (2021) He, S., G. Jiang, H. Lam, and M. C. Fu. 2021.  “Adaptive Importance Sampling for Efficient Stochastic Root Finding and Quantile Estimation”.  arXiv preprint arXiv:2102.10631.
  * Heidelberger (1995) Heidelberger, P. 1995, January.  “Fast Simulation of Rare Events in Queueing and Reliability Models”.  ACM Trans. Model. Comput. Simul. 5(1):43–85.
  * Juneja and Shahabuddin (2006) Juneja, S., and P. Shahabuddin. 2006.  “Chapter 11 Rare-Event Simulation Techniques: An Introduction and Recent Advances”.  In Simulation, edited by S. G. Henderson and B. L. Nelson, Volume 13 of Handbooks in Operations Research and Management Science, 291 – 350. Elsevier.
  * Lemaire and Pagès (2010) Lemaire, V., and G. Pagès. 2010.  “Unconstrained recursive importance sampling”.  The Annals of Applied Probability 20(3):1029–1067.
  * McNeil et al. (2015) McNeil, A. J., R. Frey, and P. Embrechts. 2015.  Quantitative Risk Management: Concepts, Techniques and Tools.  Princeton university press.
  * Pasupathy (2010) Pasupathy, R. 2010.  “On choosing parameters in retrospective-approximation algorithms for stochastic root finding and simulation optimization”.  Operations Research 58(4-part-1):889–901.
  * Rockafellar and Uryasev (2000) Rockafellar, R. T., and S. Uryasev. 2000.  “Optimization of Conditional Value-at-Risk”.  Journal of Risk 2:21–42.
  * Rockafellar and Wets (2009) Rockafellar, R. T., and R. J.-B. Wets. 2009.  Variational analysis, Volume 317.  Springer Science & Business Media.
  * Rubinstein (1997) Rubinstein, R. Y. 1997.  “Optimization of computer simulation models with rare events”.  European Journal of Operational Research 99(1):89–112.
  * Shapiro (1991) Shapiro, A. 1991.  “Asymptotic analysis of stochastic programs”.  Annals of Operations Research.
  * Varadhan (1988) Varadhan, S. R. S. 1988.  “Large deviations and Applications”.  In École d’Été de Probabilités de Saint-Flour XV–XVII, 1985–87, edited by P.-L. Hennequin, 1–49.  Berlin, Heidelberg: Springer Berlin Heidelberg.

## AUTHOR BIOGRAPHIES

ANAND DEO is a Postdoctoral Researcher at Singapore University of Technology
and Design. His research interests are in quantitative risk management,
operations research and machine learning. His e-mail address is
<EMAIL_ADDRESS>
KARTHYEK MURTHY is an Assistant Professor in Singapore University of
Technology and Design. His research centers around building models and methods
for incorporating competing considerations such as risk, robustness, and
fairness in data-driven optimization problems affected by uncertainty. His
e-mail address is<EMAIL_ADDRESS>
TIRTHO SARKER is a Research Assistant at Singapore University of Technology
and Design. His e-mail address is<EMAIL_ADDRESS>


$\chi_{\mathcal{F}\cup\\{P^{r}\\},G,v}^{(e+1)}=\chi_{\mathcal{F}\cup\\{P^{r}\\},H,w}^{(e+1)}$.
We again use the fact that we can associate with every $v\in V_{G}$ a unique
vertex $w\in V_{H}$ such that
$\chi_{\mathcal{F}\cup\\{P_{1}^{r},P_{2}^{r}\\},G,v}^{(d)}=\chi_{\mathcal{F}\cup\\{P_{1}^{r},P_{2}^{r}\\},H,w}^{(d)}$.
In particular, this implies that
$\chi_{\mathcal{F}\cup\\{P_{1}^{r},P_{2}^{r}\\},G,v}^{(e)}=\chi_{\mathcal{F}\cup\\{P_{1}^{r},P_{2}^{r}\\},H,w}^{(e)}$
and
$\chi_{\mathcal{F}\cup\\{P_{1}^{r},P_{2}^{r}\\},G,v}^{(e+1)}=\chi_{\mathcal{F}\cup\\{P_{1}^{r},P_{2}^{r}\\},H,w}^{(e+1)}$.
From the definition of the $\text{-}\mathsf{WL}$-test, it must also be the
case that the multisets
$\\{\chi_{\mathcal{F}\cup\\{P_{1}^{r},P_{2}^{r}\\},G,v^{\prime}}^{(e)}\mid
v^{\prime}\in N_{G}(v)\\}$ and
$\\{\chi_{\mathcal{F}\cup\\{P_{1}^{r},P_{2}^{r}\\},H,w^{\prime}}^{(e)}\mid
v^{\prime}\in N_{H}(w)\\}$ must be equal as well, i.e., we can find a one-to-
one corresponence between neighbors of $v$ in $G$ and neighbors of $w$ in $H$
that have the same label. From the induction hypothesis we then have that
$\chi_{\mathcal{F}\cup\\{P^{r}\\},G,v}^{(e)}=\chi_{\mathcal{F}\cup\\{P^{r}\\},H,w}^{(e)}$
and also that the multisets
$\\{\chi_{\mathcal{F}\cup\\{P^{r}\\},G,v^{\prime}}^{(e)}\mid v^{\prime}\in
N_{G}(v)\\}$ and $\\{\chi_{\mathcal{F}\cup\\{P^{r}\\},H,w^{\prime}}^{(e)}\mid
v^{\prime}\in N_{H}(w)\\}$ are equal, which implies, by the definition of the
$\mathsf{WL}$-test, that
$\chi_{\mathcal{F}\cup\\{P^{r}\\},G,v}^{(e+1)}=\chi_{\mathcal{F}\cup\\{P^{r}\\},H,w}^{(e+1)}$,
as was to be shown.

###  C.3 Proof of Theorem 3

We show that $\mathcal{F}\cup\\{Q^{r}\\}\text{-}\mathsf{WL}$, where $Q^{r}$ is
pattern whose core has treewidth $k$, is more expressive than
$\mathcal{F}\text{-}\mathsf{WL}$ if every pattern $P^{r}\in\mathcal{F}$
satisfies one of the following conditions: (i) $P^{r}$ has treewidth $<k$; or
(ii) $P^{r}$ does not map homomorphically to $Q^{r}$.

Let $c(Q)^{r}$ to denote the (rooted) core of $Q$, in which the root of
$c(Q)^{r}$ is any vertex which is the image of the root of $Q^{r}$ in a
homomorphism from $Q^{r}$ to $c(Q)^{r}$. By assumption, $c(Q)^{r}$ has
treewidth $k$.

Clearly, $\mathcal{F}\text{-}\mathsf{WL}$ is upper bounded by
$\mathcal{F}\cup\\{Q^{r}\\}\text{-}\mathsf{WL}$. Thus, all we need for the
proof is to find two graphs that are indistinguishable by
$\mathcal{F}\text{-}\mathsf{WL}$ but are in fact distinguished by
$\mathcal{F}\cup\\{Q^{r}\\}\text{-}\mathsf{WL}$.

Those two graphs are, in fact, the graphs $G$ and $H$ constructed for
$c(Q)^{r}$ (of treewidth $k$) in the proof of Theorem 2. From that proof, we
know that:

  1. (a)

$\mathsf{hom}(c(Q),G)=0$ and $\mathsf{hom}(c(Q),H)\neq 0$; and

  2. (b)

$G\equiv_{\mathsf{C}^{k}}H$.

We note that (a) immediately implies that $G$ and $H$ can be distinguished by
$\mathcal{F}\cup\\{Q^{r}\\}\text{-}\mathsf{WL}$. In fact, they are
distinguished in already by the initial labelling in round $0$. We next show
that $G$ and $H$ are indistinguishable by $\mathcal{F}\text{-}\mathsf{WL}$.

Let us first present a small structural result that helps us deal with
patterns in $\mathcal{F}$ satisfying the second condition of the Theorem.

######  Lemma 4.

If a rooted pattern $P^{r}$ does not map homomorphically to $Q^{r}$, then
$\mathsf{hom}(P,G)=\mathsf{hom}(P,H)=0$

###### Proof.

We use the following property of graphs $G$ and $H$, which can be directly
observed from their construction (and was already noted in Atserias et al.
(2007) and Bova & Chen (2019)). Define $G^{r}$ and $H^{r}$ by setting as their
root any vertex $(a_{r},f)$, for $a_{r}$ the root of $c(Q)^{r}$. Then there is
a homomorphism from $G^{r}$ to $c(Q)^{r}$, and there is a homomorphism from
$H^{r}$ to $c(Q)^{r}$.

Now, any homomorphism $h$ from $P^{r}$ to $G$ can be extended to a
homomorphism from $P^{r}$ to $Q^{r}$: we compose $h$ with the homomorphism
mentioned above from $G$ to $c(Q)^{r}$, which by definition again maps
homomorphically to $Q^{r}$. Since by definition we have that $P^{r}$ does not
map to $Q^{r}$, $h$ cannot exist. The proof for $H$ is analogous. ∎

Now, let $\mathcal{F}^{\prime}$ be the set of patterns obtained by removing
from $\mathcal{F}$ all patterns which do not map homomorphically to $Q^{r}$.
By Lemma 4, we have that $G$ and $H$ are distinguished by the
$\mathcal{F}\text{-}\mathsf{WL}$-test if and only if they are distinguished by
$\mathcal{F}^{\prime}\text{-}\mathsf{WL}$.

But all patterns in $\mathcal{F}^{\prime}$ must have treewidth less than $k$,
and by (b) $G$ and $H$ are indistinguishable by $k\text{-}\mathsf{WL}$.
Proposition 2 then implies that $G$ and $H$ are indistinguishable by
$\mathcal{F}\text{-}\mathsf{WL}$, as desired.

##  Appendix D Connections to related work

We here provide more details of how $\mathcal{F}$-$\mathsf{MPNN}\text{s}$
connect to $\mathsf{MPNN}\text{s}$ from the literature which also augment the
initial labelling.

##### Vertex degrees.

We first consider so-called degree-aware $\mathsf{MPNN}\text{s}$ (Geerts et
al., 2020) in which the message functions of the $\mathsf{MPNN}\text{s}$ may
depend on the vertex degrees. The Graph Convolution Networks ($\mathsf{GCN}$s)
(Kipf & Welling, 2017) are an example of such $\mathsf{MPNN}\text{s}$. Degree-
aware $\mathsf{MPNN}\text{s}$ are known to be equivalent, in terms of
expressive power, to standard $\mathsf{MPNN}\text{s}$ in which the initial
labelling is extended with vertex degrees (Geerts et al., 2020). Translated to
our setting, we can simply let
$\mathcal{F}=\\{\raisebox{0.0pt}{\includegraphics[height=12.91663pt]{L1.pdf}}\\}$
since
$\mathsf{hom}(\raisebox{0.0pt}{\includegraphics[height=12.91663pt]{L1.pdf}},G^{v})$
is equal to the degree of vertex $v$ in $G$. When considering graphs without
an initial vertex labelling (or a uniform labelling which assigns every vertex
the same label), our characterisation (Theorem 1) implies
$G\equiv_{\raisebox{0.0pt}{\includegraphics[height=6.02777pt]{L1.pdf}}\text{-}\mathsf{WL}}^{(d)}H$
if and only if $\mathsf{hom}(T,G)=\mathsf{hom}(T,H)$ for every
$\\{\raisebox{0.0pt}{\includegraphics[height=12.91663pt]{L1.pdf}}\\}$-pattern
tree of depth at most $d$. This in turn is equivalent to
$\mathsf{hom}(T,G)=\mathsf{hom}(T,H)$ for every (standard) tree of depth at
most $d+1$. Indeed,
$\\{\raisebox{0.0pt}{\includegraphics[height=12.91663pt]{L1.pdf}}\\}$-pattern
trees of depth at most $d$ are simply trees of depth $d+1$. Combining this
with the characterisation of $\mathsf{WL}$ by Dvorak (2010) and Dell et al.
(2018), we thus have for unlabelled graphs that
$G\equiv_{\raisebox{0.0pt}{\includegraphics[height=6.02777pt]{L1.pdf}}\text{-}\mathsf{WL}}^{(d)}H$
if and only if $G\equiv_{\mathsf{WL}}^{(d+1)}H$. So, by considering
$\mathcal{F}=\\{\raisebox{0.0pt}{\includegraphics[height=12.91663pt]{L1.pdf}}\\}$-$\mathsf{MPNN}\text{s}$
one gains one round of computation compared to considering standard
$\mathsf{MPNN}\text{s}$. To lift this to labeled graphs, instead of
$\mathcal{F}=\\{\raisebox{0.0pt}{\includegraphics[height=12.91663pt]{L1.pdf}}\\}$
one has to include labeled versions of the single edge pattern, in order to
count the number of neighbours of a specific label for each vertex. This is
done, e.g., by Ishiguro et al. (2020), who use the $\mathsf{WL}$ labelling
obtained after the first round to augment the initial vertex labelling. This
corresponds indeed by adding $\mathsf{hom}(T^{r},G^{v})$ as feature for every
labeled tree of depth one. This results in that
$G\equiv_{\raisebox{0.0pt}{\includegraphics[height=6.02777pt]{L1.pdf}}\text{-}\mathsf{WL}}^{(d)}H$
if and only if $G\equiv_{\mathsf{WL}}^{(d+1)}H$ for labelled graphs.

##### Walk counts.

The Graph Feature Networks by Chen et al. (2019a) can be regarded as a
generalisation of the previous approach. Instead of simply adding vertex
degrees, the number of walks of certain lengths emanating from vertices are
added. Translated to our setting, this corresponds to considering
$\\{L_{2},L_{3},\ldots,L_{\ell}\\}$-$\mathsf{MPNN}\text{s}$, where $L_{\ell}$
denotes a rooted path of length $\ell$. For unlabelled graphs, our
characterisation (Theorem 1) implies that
$G\equiv_{L_{1},\ldots,L_{\ell}\text{-}\mathsf{WL}}^{(d)}H$ is upper bounded
by $G\equiv_{\mathsf{WL}}^{(d+\ell)}H$, simply because every
$\\{L_{2},L_{3},\ldots,L_{\ell}\\}$-pattern tree of depth $d$ is a standard
tree of depth at most $d+\ell$.

##### Cycles.

Li et al. (2019) extend $\mathsf{MPNN}\text{s}$ by varying the notion of
neighbourhood over which is aggregated. One particular instance corresponds to
an aggregation of features, weighted by the number of cycles of a certain
length in each vertex (see discussion at the end of Section 4 in Li et al.
(2019)). Translated to our setting, this corresponds to considering
$\\{C_{\ell}\\}$-$\mathsf{MPNN}\text{s}$ where $C_{\ell}$ denotes the cycle of
length $\ell$. As mentioned in the main body of the paper, these extend
$\mathsf{MPNN}\text{s}$ and result in architectures bounded by
$2\text{-}\mathsf{WL}$ (Proposition 4). This is in line with Theorem 3 from Li
et al. (2019) stating that their framework strictly extends
$\mathsf{MPNN}\text{s}$ and thus $1\text{-}\mathsf{WL}$.

##### Isomorphism counts.

Another, albeit similar, approach to add structural information to the initial
labelling is taken in the paper Graph Structure Networks by Bouritsas et al.
(2020). The idea there is to extend the initial features with information
about how often a vertex $v$ appears in a subgraph of $G$ which is isomorphic
to $P$. More precisely, Bouritsas et al. (2020) consider a connected
unlabelled graph $P$ as pattern and partition its vertex set $V_{P}$ orbit-
wise. That is, $V_{P}=\biguplus_{i=1}^{o_{P}}V_{P}^{i}$ where $o_{P}$ denotes
the number of orbits of $P$. Here, $v,v^{\prime}\in V_{P}^{i}$ whenever there
is an automorphism $h$ in $\mathsf{Aut}(P)$ mapping $v$ to $v^{\prime}$. Next,
they consider all distinct subgraphs $G_{1},\ldots,G_{k}$ in $G$ which are
isomorphic to $P$, denoted by $P\cong G_{j}$ for $j\in[k]$. We write
$P\cong_{f}G_{j}$ when $P\cong G_{j}$ using a specific isomorphism $f$. Then
for each orbit partition $i\in[o_{P}]$ and vertex $v\in V$, they define:

$\mathsf{iso}(P,G,v,i)=|\\{G_{j}\cong P\mid v\in V_{G_{j}},\text{and there
exists an $f$ s.t. }G_{j}\cong_{f}P\text{ and }f(v)\in V_{P}^{i},j\in[k]\\}|.$

That is, the number of subgraphs $G_{j}$ in $G$ that can be isomorphically
mapped to $P$ are counted, provided that this can be done by an isomorphism
which maps vertex $v$ in $G_{j}$ (and thus $G$) to one of the vertices in the
$i$th orbit partition $V_{P}^{i}$ of the pattern. A similar notion is proposed
for edges, which we will not consider here. Similar to our extended features,
the initial features of each vertex $v$ is then augmented with
$\big{(}\mathsf{iso}(P,G^{v},i)\mid P\in\mathcal{F},i\in[o_{P}]\big{)}$ for
some set $\mathcal{F}$ of patterns. Standard $\mathsf{MPNN}\text{s}$ are
executed on these augmented initial features. We refer to Bouritsas et al.
(2020) for more details.

We can view the above approach as an instance of our framework. Indeed, given
a pattern $P$ in $\mathcal{F}$, for each orbit partition, we replace $P$ by a
different rooted version $P^{r_{i}}$, where $r_{i}$ is a vertex in
$V_{P}^{i}$. Which vertex in the orbit under consideration is selected as root
is not important (because they are equivalent by definition of orbit). We then
see that the standard notion of subgraph isomorphism counting directly
translates to the quantity used in Bouritsas et al. (2020):

$\mathsf{sub}(P^{r_{i}},G^{v}):=\text{number of subgraphs in $G$ containing
$v$, isomorphic to $P^{r_{i}}$}=\mathsf{iso}(P,G,v,i).$

It thus remains to express $\mathsf{sub}(P^{r_{i}},G^{v})$ in terms of
homomorphism counts. This, however, follows from Curticapean et al. (2017) in
which it is shown that $\mathsf{iso}(P^{r_{i}},G^{v})$ can be computed by a
linear combination of $\mathsf{hom}(Q^{r_{i}},G^{v})$ where $Q^{r_{i}}$ ranges
over all graphs on which $P^{r_{i}}$ can be mapped by means of a surjective
homomorphism. For a given $P^{r_{i}}$, the finite set of such patterns is
called the spasm of $P^{r_{i}}$ in Curticapean et al. (2017) and can be easily
computed.

In summary, given the set $\mathcal{F}$ of patterns in Bouritsas et al.
(2020), we first replace every $P\in\mathcal{F}$ by its rooted versions
$P^{r_{i}}$, for $i\in[o_{P}]$, and then expand the resulting set of rooted
patterns, by the spasms of each of these patterns. Let us denote by
$\mathcal{F}^{*}$ the final set of rooted patterns. It now follows that
$\mathsf{hom}(Q^{r},G^{v})$ for $Q^{r}\in\mathcal{F}^{*}$ provides sufficient
information to extract $\mathsf{sub}(P^{r_{i}},G^{v})$ and thus also
$\mathsf{iso}(P,G,v,i)$ for every $P\in\mathcal{F}$ and orbit part
$i\in[o_{P}]$. As a consequence, the $\mathsf{MPNN}\text{s}$ from Bouritsas et
al. (2020) are bounded by $\mathcal{F}^{*}$-$\mathsf{MPNN}\text{s}$ and thus
$\mathcal{F}^{*}\text{-}\mathsf{WL}$. Conversely, given an
$\mathcal{F}$-$\mathsf{MPNN}$ one can, again using results by Curticapean et
al. (2017), define a set $\mathcal{F}^{+}$ of patterns, such that the subgraph
isomorphism counts of patterns in $\mathcal{F}^{+}$ can be used to compute the
homomorphism counts of patterns in $\mathcal{F}$. Hence,
$\mathcal{F}$-$\mathsf{MPNN}\text{s}$ are upper bounded by the
$\mathsf{MPNN}\text{s}$ considered in Bouritsas et al. (2020) using patterns
in $\mathcal{F}^{+}$. This is in agreement with Curticapean et al. (2017) in
which it is shown that homomorphism counts, subgraph isomorphism counts and
other notions of pattern counts are all interchangeable. Nevertheless, by
using homomorphism counts one can gracefully extend known results about
$\mathsf{WL}$ and $\mathsf{MPNN}\text{s}$, as we have shown in the paper, and
add little overhead.

##  Appendix E Additional experimental information

###  E.1 Experimental setup

One of the crucial questions when studying the effect of adding structural
information to the initial vertex labels is whether these additional labels
enhance the performance of graph neural networks. In order to reduce the
effect of specific implementation details of $\mathsf{GNN}\text{s}$ and choice
of hyper-parameters, we start from the $\mathsf{GNN}$ implementations and
choices made in the benchmark by Dwivedi et al. (2020)777The original
implementations can be found on
https://github.com/graphdeeplearning/benchmarking-gnns. and only change the
initial vertex labels, while leaving the $\mathsf{GNN}\text{s}$ themselves
unchanged. This ensures that we only measure the effect of augmenting initial
features with homomorphism counts. We will use the $\mathsf{GNN}\text{s}$ from
the benchmark, without extended features, as our baselines. For the same
reasons, we use datasets proposed in the benchmark for their ability to
statistically separate the performance of $\mathsf{GNN}\text{s}$. All other
parameters are taken as in Dwivedi et al. (2020) and we refer to that paper
for more details.

##### Selected $\mathsf{GNN}\text{s}$

Dwivedi et al. (2020) divide the benchmarked $\mathsf{GNN}\text{s}$ into two
classes: the $\mathsf{MPNN}\text{s}$ and the “theoretically designed”
$\mathsf{WL}$-$\mathsf{GNN}\text{s}$. The first class is found to perform
stronger and train faster. Hence, we chose to include the five following
$\mathsf{MPNN}$ models from the benchmark:

  * •

Graph Attention Network ($\mathsf{GAT}$) as described in Velickovic et al.
(2018)

  * •

Graph Convolutional Network ($\mathsf{GCN}$) as described in Kipf & Welling
(2017)

  * •

$\mathsf{GraphSage}$ as described in Hamilton et al. (2017)

  * •

Mixed Model Convolutional Networks ($\mathsf{MoNet}$) as described in Monti et
al. (2017)

  * •

$\mathsf{GatedGCN}$ as described in Bresson & Laurent (2017).

For $\mathsf{GatedGCN}$ we used the version in which positional encoding
(Belkin & Niyogi, 2003) is added to the vertex features, as it is empirically
shown to be the strongest performing version of this model by for the selected
datasets (Dwivedi et al., 2020). We denote this version by
$\mathsf{GatedGCN}_{E,PE}$, referring to the presence of edge features and
this positional encoding. Details, background and a mathematical formalisation
of the message passing layers of these models can be found in the
supplementary material of Dwivedi et al. (2020).

As explained in the experimental section of the main paper, we enhance the
vertex features with the log-normalised counts of the chosen patterns in every
vertex of every graph of every dataset. The first layers of some models of
(Dwivedi et al., 2020) are adapted to take in this variation in input size.
All other layers where left identical to their original implementation as
provided by Dwivedi et al. (2020).

##### Hardware, compute and resources

All models for ZINC, PATTERN and COLLAB were trained on a GeForce GTX 1080 337
Ti GPU, for CLUSTER a Tesla V100-SXM3-32GB GPU was used. Tables 7, 10, 13 and
16 report the training times for all combination of models and additional
feature set. A rough estimate of the CO2 emissions based on the total
computing times of reported experiments ($2\,074$ hours GeForce GTX 1080,
$372$ hours Tesla V100-SXM3-32GB), the computing times of not-included
experiments ($1\,037$ hours GeForce GTX 1080, $181$ hours Tesla
V100-SXM3-32GB), the GPU types (GeForce GTX 1080, Tesla V100-SXM3-32GB) and
the geographical location of our cluster results in a carbon emission of $135$
kg CO2 equivalent. This estimation was conducted using the MachineLearning
Impact calculator presented in Lacoste et al. (2019).

###  E.2 Graph learning tasks

We here report the full results of our experimental evaluation for graph
regression (Section E.2.1), link prediction (Section E.2.2) and vertex
classification (Section E.2.3) as considered in Dwivedi et al. (2020). More
precisely, a full listing of the patterns and combinations used and the
obtained results for the test sets can be found in Tables 5, 8, 11 and 14.
Average training time (in hours) and the number of epochs are reported in
Tables 7, 10, 13 and 16. Finally, the total number of model parameters are
reported in Tables 6, 9, 12 and 15. All averages and standard deviations are
over 4 runs with different random seeds. The main take-aways from these
results are included in the main paper.

####  E.2.1 Graph regression with the ZINC dataset

Just as in Dwivedi et al. (2020) we use a subset ($12$K) of ZINC molecular
graphs ($250$K) dataset (Irwin et al., 2012b) to regress a molecular property
known as the constrained solubility. For each molecular graph, the vertex
features are the types of heavy atoms and the edge features are the types of
bonds between them. The following are taken from Dwivedi et al. (2020):
Splitting. ZINC has $10\,000$ train, $1\,000$ validation and $1\,000$ test
graphs.
Training.888Here and in the next tasks we are using the parameters used in the
code accompanying Dwivedi et al. (2020). In the paper, slightly different
parameters are used. For the learning rate strategy, an initial learning rate
is set to $5\times 10^{-5}$ , the reduce factor is $0.5$, and the stopping
learning rate is $1\times 10^{-6}$, the patience value is 25 and the maximal
training time is set to $12$ hours.
Performance Measure The performance measure is the mean absolute error (MAE)
between the predicted and the ground truth constrained solubility for each
molecular graph.
Number of layers 16 MPNN layers are used for every model.

Table 5: Full results of the mean absolute error (predicted constrained solubility vs. the ground truth) for selected cycle combinations and $\mathsf{GNN}\text{s}$ on the ZINC data set. In the last two rows we compare between homomorphism counts (hom) and subgraph isomorphism counts (iso). Pattern set $\mathcal{F}$ | $\mathsf{GAT}$ | $\mathsf{GCN}$ | $\mathsf{GraphSage}$ | $\mathsf{MoNet}$ | $\mathsf{GatedGCN}_{E,PE}$
---|---|---|---|---|---
None | 0,47$\pm$0,02 | 0,35$\pm$0,01 | 0,25$\pm$0,01 | 0,44$\pm$0,01 | 0,34$\pm$0,05
$\\{C_{3}\\}$ | 0,45$\pm$0,01 | 0,36$\pm$0,01 | 0,25$\pm$0,00 | 0,44$\pm$0,00 | 0,30$\pm$0,01
$\\{C_{4}\\}$ | 0,34$\pm$0,02 | 0,29$\pm$0,02 | 0,26$\pm$0,01 | 0,30$\pm$0,01 | 0,27$\pm$0,06
$\\{C_{5}\\}$ | 0,44$\pm$0,02 | 0,34$\pm$0,02 | 0,23$\pm$0,01 | 0,42$\pm$0,01 | 0,27$\pm$0,03
$\\{C_{6}\\}$ | 0,31$\pm$0,00 | 0,27$\pm$0,02 | 0,25$\pm$0,01 | 0,30$\pm$0,01 | 0,26$\pm$0,09
$\\{C_{3},C_{4}\\}$ | 0,33$\pm$0,01 | 0,27$\pm$0,01 | 0,24$\pm$0,02 | 0,32$\pm$0,01 | 0,23$\pm$0,03
$\\{C_{5},C_{6}\\}$ | 0,28$\pm$0,01 | 0,26$\pm$0,01 | 0,23$\pm$0,01 | 0,28$\pm$0,01 | 0,20$\pm$0,03
$\\{C_{4},C_{5},C_{6}\\}$ | 0,24$\pm$0,00 | 0,21$\pm$0,00 | 0,20$\pm$0,00 | 0,25$\pm$0,01 | 0,16$\pm$0,02
$\\{C_{3},C_{4},C_{5},C_{6}\\}$ | 0,23$\pm$0,00 | 0,21$\pm$0,00 | 0,20$\pm$0,01 | 0,26$\pm$0,02 | 0,18$\pm$0,02
$\\{C_{3},\ldots,C_{10}\\}$ (hom) | 0,22$\pm$0,01 | 0,20$\pm$0,00 | 0,19$\pm$0,00 | 0,2376$\pm$0,01 | 0,1352$\pm$0,01
$\\{C_{3},\ldots,C_{10}\\}$ (iso) | 0,24$\pm$0,01 | 0,22$\pm$0,01 | 0,16$\pm$0,01 | 0,2408$\pm$0,01 | 0,1357 $\pm$ 0,01
Table 6: Total model parameters for selected cycle combinations and $\mathsf{GNN}\text{s}$ on the ZINC data set. In the last two rows we compare between homomorphism counts (hom) and subgraph isomorphism counts (iso). Pattern set $\mathcal{F}$ | $\mathsf{GAT}$ | $\mathsf{GCN}$ | $\mathsf{GraphSage}$ | $\mathsf{MoNet}$ | $\mathsf{GatedGCN}_{E,PE}$
---|---|---|---|---|---
None | 358 273 | 360 742 | 388 963 | 401 148 | 408 135
$\\{C_{3}\\}$ | 358 417 | 360 887 | 389 071 | 401 238 | 408 205
$\\{C_{4}\\}$ | 358 417 | 360 887 | 389 071 | 401 238 | 408 205
$\\{C_{5}\\}$ | 358 417 | 360 887 | 389 071 | 401 238 | 408 205
$\\{C_{6}\\}$ | 358 417 | 360 887 | 389 071 | 401 238 | 408 205
$\\{C_{3},C_{4}\\}$ | 358 561 | 361 032 | 389 179 | 401 328 | 408 275
$\\{C_{5},C_{6}\\}$ | 358 561 | 361 032 | 389 179 | 401 328 | 408 275
$\\{C_{4},C_{5},C_{6}\\}$ | 358 705 | 361 177 | 389 287 | 401 418 | 408 345
$\\{C_{3},C_{4},C_{5},C_{6}\\}$ | 358 849 | 361 322 | 389 395 | 401 508 | 408 415
$\\{C_{3},\ldots,C_{10}\\}$ (hom) | 359 425 | 361 902 | 389 827 | 401 868 | 408 695
$\\{C_{3},\ldots,C_{10}\\}$ (iso) | 359 425 | 361 902 | 389 827 | 401 868 | 408 695
Table 7: Average training time in hours and number of epochs for selected
cycle combinations and $\mathsf{GNN}\text{s}$ on the ZINC data set. In the
last two rows we compare between homomorphism counts (hom) and subgraph
isomorphism counts (iso).

Model: $\mathsf{GAT}$ $\mathsf{GCN}$ $\mathsf{GraphSage}$ $\mathsf{MoNet}$
$\mathsf{GatedGCN}_{E,PE}$ Pattern set $\mathcal{F}$ Time Epochs Time Epochs
Time Epochs Time Epochs Time Epochs None 2,40 377 10,99 463 2,46 420 1,53 345
12,08 136 $\\{C_{3}\\}$ 2,88 444 12,03 363 2,03 500 0,91 298 12,07 148
$\\{C_{4}\\}$ 2,30 351 11,36 324 2,31 396 1,70 382 12,06 139 $\\{C_{5}\\}$
2,42 375 12,03 333 1,70 444 1,06 370 12,06 202 $\\{C_{6}\\}$ 2,40 369 9,98 421
2,58 446 1,25 288 12,08 136 $\\{C_{3},C_{4}\\}$ 2,98 461 12,03 332 2,56 458
1,41 321 12,09 132 $\\{C_{5},C_{6}\\}$ 2,76 422 12,04 319 2,67 464 1,53 356
12,06 137 $\\{C_{4},C_{5},C_{6}\\}$ 2,45 381 10,13 419 1,67 463 1,04 382 12,04
229 $\\{C_{3},C_{4},C_{5},C_{6}\\}$ 2,65 408 10,38 420 2,09 503 1,26 364 12,08
135 $\\{C_{3},\ldots,C_{10}\\}$ (hom) 2,65 428 12,03 350 2,76 478 1,48 363
12,06 175 $\\{C_{3},\ldots,C_{10}\\}$ (iso) 2,78 497 11,72 419 2,63 547 1,58
440 11,62 148

####  E.2.2 Link Prediction with the Collab dataset

Another set used in Dwivedi et al. (2020) is COLLAB, a link prediction dataset
proposed by the Open Graph Benchmark (OGB) (Hu et al., 2020) corresponding to
a collaboration network between approximately $235$K scientists, indexed by
Microsoft Academic Graph. Vertices represent scientists and edges denote
collaborations between them. For vertex features, OGB provides
$128$-dimensional vectors, obtained by averaging the word embeddings of a
scientist’s papers. The year and number of co-authored papers in a given year
are concatenated to form edge features. The graph can also be viewed as a
dynamic multi-graph, since two vertices may have multiple temporal edges
between if they collaborate over multiple years. The following are taken from
Dwivedi et al. (2020):

Splitting. We use the real-life training, validation and test edge splits
provided by OGB. Specifically, they use collaborations until 2017 as training
edges, those in 2018 as validation edges, and those in 2019 as test edges.
Training. All GNNs use the same learning rate strategy: an initial learning
rate is set to $1\times 10^{-3}$, the reduce factor is $0.5$, the patience
value is 10, and the stopping learning rate is $1\times 10^{-5}$.
Performance Measure. We use the evaluator provided by OGB (Hu et al., 2020),
which aims to measure a model’s ability to predict future collaboration
relationships given past collaborations. Specifically, they rank each true
collaboration among a set of $100\,000$ randomly-sampled negative
collaborations, and count the ratio of positive edges that are ranked at
$K$-place or above (Hits@K). The value $K=50$ as this gives the best value for
statistically separating the performance of $\mathsf{GNN}\text{s}$. Number of
layers 3 MPNN layers are used for every model.

Table 8: Full Results (Hits @50) for all selected pattern combinations and $\mathsf{GNN}\text{s}$ on the COLLAB data set. Pattern set $\mathcal{F}$ | $\mathsf{GAT}$ | $\mathsf{GCN}$ | $\mathsf{GraphSage}$ | $\mathsf{MoNet}$ | $\mathsf{GatedGCN}_{E,PE}$
---|---|---|---|---|---
None | 50,32$\pm$0,55 | 51,36$\pm$1,30 | 49,81$\pm$1,56 | 50,33$\pm$0,68 | 51,00$\pm$2,54
$\\{K_{3}\\}$ | 52,87$\pm$0,87 | 53,57$\pm$0,89 | 50,18$\pm$1,38 | 51,10$\pm$0,38 | 51,57$\pm$0,68
$\\{K_{4}\\}$ | 51,33$\pm$1,42 | 52,84$\pm$1,32 | 51,76$\pm$1,38 | 51,13$\pm$1,60 | 49,43$\pm$1,85
$\\{K_{5}\\}$ | 52,41$\pm$0,89 | 54,60$\pm$1,01 | 50,94$\pm$1,30 | 51,39$\pm$1,23 | 50,31$\pm$1,59
$\\{K_{3},K_{4}\\}$ | 52,68$\pm$1,82 | 53,49$\pm$1,35 | 50,88$\pm$1,73 | 50,97$\pm$0,68 | 51,36$\pm$0,92
$\\{K_{3},K_{4},K_{5}\\}$ | 51,81$\pm$1,17 | 54,32$\pm$1,02 | 49,94$\pm$0,23 | 51,01$\pm$1,00 | 51,11$\pm$1,06
Table 9: Total number of model parameters for all selected pattern combinations and $\mathsf{GNN}\text{s}$ on the COLLAB data set. Pattern set $\mathcal{F}$ | $\mathsf{GAT}$ | $\mathsf{GCN}$ | $\mathsf{GraphSage}$ | $\mathsf{MoNet}$ | $\mathsf{GatedGCN}_{E,PE}$
---|---|---|---|---|---
None | 25 992 | 40 479 | 39 751 | 26 487 | 27 440
$\\{K_{3}\\}$ | 26 049 | 40 553 | 39 804 | 26 525 | 27 475
$\\{K_{4}\\}$ | 26 049 | 40 553 | 39 804 | 26 525 | 27 475
$\\{K_{5}\\}$ | 26 049 | 40 553 | 39 804 | 26 525 | 27 475
$\\{K_{3},K_{4}\\}$ | 26 106 | 40 627 | 39 857 | 26 563 | 27 510
$\\{K_{3},K_{4},K_{5}\\}$ | 26 163 | 40 701 | 39 910 | 26 601 | 27 545
Table 10: Average training times and number of epochs for all selected pattern
combinations and $\mathsf{GNN}\text{s}$ on the COLLAB data set.

Model: $\mathsf{GAT}$ $\mathsf{GCN}$ $\mathsf{MoNet}$ $\mathsf{GraphSage}$
$\mathsf{GatedGCN}_{E,PE}$ Pattern set $\mathcal{F}$ Time #Epochs Time #Epochs
Time #Epochs Time #Epochs Time #Epochs None 0,81 167 0,85 141 1,62 190 12,05
115,67 2,22 167 $\\{K_{3}\\}$ 0,67 165 0,90 153 1,70 184 12,10 67,00 2,48 186
$\\{K_{4}\\}$ 1,06 188 0,95 160 2,16 188 12,04 113,50 1,26 188 $\\{K_{5}\\}$
0,50 167 1,13 165 1,04 193 12,05 124,00 1,82 174 $\\{K_{3},K_{4}\\}$ 1,20 189
0,86 128 2,15 189 12,05 113,25 1,51 183 $\\{K_{3},K_{4},K_{5}\\}$ 0,44 149
0,90 134 0,98 186 12,05 124,00 1,84 177

####  E.2.3 Vertex classification with PATTERN and CLUSTER

Finally, also used in Dwivedi et al. (2020) are the PATTERN and CLUSTER graph
data sets, generated with the Stochastic Block Model (SBM) (Abbe, 2018), which
is widely used to model communities in social networks by modulating the
intra- and extra-communities connections, thereby controlling the difficulty
of the task. A SBM is a random graph which assigns communities to each vertex
as follows: any two vertices are connected with probability $p$ if they belong
to the same community, or they are connected with probability $q$ if they
belong to different communities (the value of $q$ acts as the noise level).

For the PATTERN dataset, the goal of the vertex classification problem is the
detection of a certain pattern $P$ embedded in a larger graph $G$. The graphs
in $G$ consist of $5$ communities with sizes randomly selected between
$[5,35]$. The parameters of the SBM for each community is $p=0.5$, $q=0.35$,
and the vertex features in $G$ are generated using a uniform random
distribution with a vocabulary of size $3$, i.e., $\\{0,1,2\\}$. Randomly,
$100$ patterns $P$ composed of $20$ vertices with intra-probability
$p_{P}=0.5$ and extra-probability $q_{P}=0.5$ are generated (i.e., $50\%$ of
vertices in $P$ are connected to $G$). The vertex features for $P$ are also
generated randomly using values in $\\{0,1,2\\}$. The graphs consist of
$44$-$188$ vertices. The output vertex labels have value $1$ if the vertex
belongs to $P$ and value $0$ belongs to $G$.

For the CLUSTER dataset, the goal of the vertex classification is the
detection of which cluster a vertex belongs. Here, six SBM clusters are
generated with sizes randomly selected between $[5,35]$ and probabilities
$p=0.55$ and $q=0.25$. The graphs consist of $40$-$190$ vertices. Each vertex
can take an initial feature value in range $\\{0,1,2,\ldots,6\\}$. If the
value is $i$ then the vertex belongs to class $i-1$. If the value is $0$, then
the class of the vertex is unknown and need to be inferred. There is only one
labelled vertex that is randomly assigned to each community and most vertex
features are set to $0$. The output vertex labels are defined as the
community/cluster class labels.

The following are taken from Dwivedi et al. (2020):

Splitting The PATTERN dataset has $10\,000$ train, $2\,000$ validation and
$2\,000$ test graphs. The CLUSTER dataset has $10\,000$ train, $1\,000$
validation and $1\,000$ test graphs. We save the generated splits and use the
same sets in all models for fair comparison.
Training For all $\mathsf{GNN}\text{s}$, an initial learning rate is set to
$1\times 10^{-3}$, the reduce factor is $0.5$, the patience value is $10$, and
the stopping learning rate is $1\times 10^{-5}$ .
Performance measure The performance measure is the average vertex-level
accuracy weighted with respect to the class sizes. Number of layers 16 MPNN
layers are used for every model.

Table 11: Full results of the weighted accuracy for selected pattern combinations and $\mathsf{GNN}\text{s}$ on the CLUSTER data set. Pattern set $\mathcal{F}$ | $\mathsf{GAT}$ | $\mathsf{GCN}$ | $\mathsf{MoNet}$ | $\mathsf{GraphSage}$ | $\mathsf{GatedGCN}_{E,PE}$
---|---|---|---|---|---
None | 70,86$\pm$0,06 | 70,64$\pm$0,39 | 71,15$\pm$0,33 | 72,25$\pm$0,52 | 74,28$\pm$0,15
$\\{K_{3}\\}$ | 71,60$\pm$0,15 | 64,88$\pm$4,16 | 72,21$\pm$0,19 | 72,97$\pm$0,23 | 74,14$\pm$0,12
$\\{K_{4}\\}$ | 71,40$\pm$0,24 | 60,64$\pm$2,93 | 72,14$\pm$0,19 | 72,57$\pm$0,19 | 74,16$\pm$0,24
$\\{K_{5}\\}$ | 71,26$\pm$0,39 | 66,60$\pm$1,47 | 72,34$\pm$0,09 | 72,60$\pm$0,24 | 74,23$\pm$0,07
$\\{K_{3},K_{4}\\}$ | 71,80$\pm$0,28 | 50,94$\pm$22,98 | 72,32$\pm$0,27 | 73,03$\pm$0,25 | 74,17$\pm$0,13
$\\{K_{3},K_{4},K_{5}\\}$ | 71,63$\pm$0,26 | 63,03$\pm$3,72 | 72,32$\pm$0,36 | 72,65$\pm$0,13 | 74,03$\pm$0,19
Table 12: Total number of model parameters for all selected pattern combinations and $\mathsf{GNN}\text{s}$ on the CLUSTER data set. Pattern set $\mathcal{F}$ | $\mathsf{GAT}$ | $\mathsf{GCN}$ | $\mathsf{MoNet}$ | $\mathsf{GraphSage}$ | $\mathsf{GatedGCN}_{E,PE}$
---|---|---|---|---|---
None | 395 396 | 362 849 | 399 373 | 386 835 | 406 755
None | 395 396 | 362 849 | 399 373 | 386 835 | 406 755
$\\{K_{3}\\}$ | 395 396 | 362 849 | 399 373 | 386 835 | 406 755
$\\{K_{4}\\}$ | 395 548 | 362 995 | 399 463 | 386 943 | 406 825
$\\{K_{5}\\}$ | 395 700 | 363 141 | 399 553 | 387 051 | 406 895
$\\{K_{3},K_{4}\\}$ | 395 700 | 363 141 | 399 553 | 387 051 | 406 895
$\\{K_{3},K_{4},K_{5}\\}$ | 396 004 | 363 433 | 399 733 | 387 267 | 407 035
Table 13: Training times (in hours) and number of epochs for all selected
pattern combinations and $\mathsf{GNN}\text{s}$ on the CLUSTER data set.

Model: GAT $\mathsf{GCN}$ $\mathsf{MoNet}$ $\mathsf{GraphSage}$
$\mathsf{GatedGCN}_{E,PE}$ Pattern set $\mathcal{F}$ Time #Epochs Time #Epochs
Time #Epochs Time #Epochs Time #Epochs None 1,62 109 2,83 117 1,54 125 0,95
101 10,40 92 $\\{K_{3}\\}$ 1,52 107 2,67 85 1,72 145 1,08 102 11,01 89
$\\{K_{4}\\}$ 1,18 107 1,94 80 1,62 149 0,90 102 10,23 90 $\\{K_{5}\\}$ 1,23
106 2,30 84 1,68 143 0,92 99 10,68 91 $\\{K_{3},K_{4}\\}$ 1,53 102 1,97 82
1,89 153 0,94 99 10,80 90 $\\{K_{3},K_{4},K_{5}\\}$ 1,62 105 1,96 82 1,95 157
0,97 100 10,25 91

Table 14: Full results of the weighted accuracy for selected pattern combinations and $\mathsf{GNN}\text{s}$ on the PATTERN data set. Pattern set $\mathcal{F}$ | $\mathsf{GAT}$ | $\mathsf{GCN}$ | $\mathsf{MoNet}$ | $\mathsf{GraphSage}$ | $\mathsf{GatedGCN}_{E,PE}$
---|---|---|---|---|---
None | 78,83$\pm$0,60 | 71,42$\pm$1,38 | 85,90$\pm$0,03 | 70,78$\pm$0,19 | 86,15$\pm$0,08
$\\{K_{3}\\}$ | 84,34$\pm$0,09 | 61,54$\pm$2,20 | 86,59$\pm$0,02 | 84,75$\pm$0,11 | 85,02$\pm$0,20
$\\{K_{4}\\}$ | 84,43$\pm$0,40 | 63,40$\pm$1,55 | 86,60$\pm$0,02 | 84,51$\pm$0,06 | 85,40$\pm$0,28
$\\{K_{5}\\}$ | 83,47$\pm$0,11 | 64,18$\pm$3,88 | 86,57$\pm$0,02 | 83,73$\pm$0,10 | 85,63$\pm$0,22
$\\{K_{3},K_{4}\\}$ | 85,44$\pm$0,24 | 81,29$\pm$2,82 | 86,58$\pm$0,02 | 85,85$\pm$0,13 | 85,80$\pm$0,20
$\\{K_{3},K_{4},K_{5}\\}$ | 85,50$\pm$0,23 | 82,49$\pm$0,48 | 86,63$\pm$0,03 | 85,88$\pm$0,15 | 85,56$\pm$0,33
Table 15: Total number of model parameters for selected pattern combinations and $\mathsf{GNN}\text{s}$ on the PATTERN data set. Pattern set $\mathcal{F}$ | $\mathsf{GAT}$ | $\mathsf{GCN}$ | $\mathsf{MoNet}$ | $\mathsf{GraphSage}$ | $\mathsf{GatedGCN}_{E,PE}$
---|---|---|---|---|---
None | 394 632 | 362 117 | 398 921 | 386 291 | 406 403
$\\{K_{3}\\}$ | 394 784 | 362 263 | 399 011 | 386 399 | 406 473
$\\{K_{4}\\}$ | 394 784 | 362 263 | 399 011 | 386 399 | 406 473
$\\{K_{5}\\}$ | 394 784 | 362 263 | 399 011 | 386 399 | 406 473
$\\{K_{3},K_{4}\\}$ | 394 936 | 362 409 | 399 101 | 386 507 | 406 543
$\\{K_{3},K_{4},K_{5}\\}$ | 395 088 | 362 555 | 399 191 | 386 615 | 406 613
Table 16: Training times (in hours) and number of epochs for selected pattern
combinations and $\mathsf{GNN}\text{s}$ on the PATTERN data set.

Model: GAT $\mathsf{GCN}$ $\mathsf{MoNet}$ $\mathsf{GraphSage}$
$\mathsf{GatedGCN}_{E,PE}$ Pattern set $\mathcal{F}$ Time Epochs Time Epochs
Time Epochs Time Epochs Time Epochs None 1,96 87 3,41 102 1,68 116 0,77 103
10,32 101 $\\{K_{3}\\}$ 0,97 97 2,58 80 1,42 107 0,69 105 9,12 95
$\\{K_{4}\\}$ 0,90 90 2,68 80 1,46 106 0,67 95 9,47 94 $\\{K_{5}\\}$ 0,89 95
2,36 80 1,26 100 0,58 98 9,14 99 $\\{K_{3},K_{4}\\}$ 2,11 91 3,62 98 1,68 108
0,86 97 9,50 87 $\\{K_{3},K_{4},K_{5}\\}$ 1,02 91 3,26 94 1,48 109 0,76 102
8,84 88