# Graph sampling for node embedding

_Li-Chun Zhang_

###### Abstract

Node embedding is a central topic in graph representation learning.
Computational efficiency and scalability can be challenging to any method that
requires full-graph operations. We propose sampling approaches to node
embedding, with or without explicit modelling of the feature vector, which aim
to extract useful information from both the eigenvectors related to the graph
Laplacien and the given values associated with the graph.

#### Keywords:

graph representation learning, graph sampling, eigen neighbour function,
normalised Laplacian, sample estimation equation

##  1 Introduction

Valued graph is a ubiquitous form of data. Graph representation learning can
be used for node classification, edge prediction, community detection and
various other purposes; see e.g. Hamilton (2020). A key to all the most
successful approaches is to utilise the links between nodes, rather than
simply treating the nodes as statistically independent units. In this context,
node embedding either refers to the generation of a low-dimension feature
vector for each node in the given graph or the generated vectors as inputs to
downstream tasks, which aims to capture the information about a node’s
position and neighbourhood structure in the graph as well as the associated
values.

Node embedding can be obtained by so-called traditional methods, which do not
require an explicit model of the features being targeted, such as based on
factorisation or random walks; see e.g. Khoshraftar and An (2022), Hamilton
(2020). More recently, the modelling approach has become very popular, where
embeddings are typically obtained using graph neural networks (GNN). See e.g.
Zhou et al. (2020) for a review of GNN methods and applications.

Figure 1 shows the well-known valued graph of Zachary’s karate club (ZKC),
where edges represent social relationships between the 34 members (i.e. nodes)
which are indicative of the members split (Zazhary, 1977). As an example of
node embedding discussed in the literature, every node can be assigned one of
four classes by modularity-based clustering (Brandes et al., 2008). Kipf and
Welling (2017) show that this community structure can be closely captured by
2-vector node embeddings, using a 3-layer graph convolutional networks (GCN),
and these embeddings are comparable to those obtained by DeepWalk (Perozzi et
al., 2014) which does not require explicit modelling.

Figure 1: Members split (encircled or not) of Zachary’s karate club

Meanwhile, the $x$-coordinates of the nodes in Figure 1 are actually given by
the eigenvector corresponding to the smallest non-zero eigenvalue of the
normalised Laplacian of this graph, which can be attributed to the traditional
factorisation-based Laplacian eigenmaps technique of Belkin and Niyogi (2001).
This illustrates clearly how embedding can facilitate node classification,
since perfect classification of the members split can be obtained using these
scalar $x_{i}$ as the node embedding, $i=1,...,34$, given a classification-
threshold parameter value between $(x_{9},x_{10})=(0.0528,0.0556)$ for the
nodes 9 and 10.

However, computational efficiency and scalability can be challenging to any of
the aforementioned approaches, if the implementation requires full-graph
operations. For instance, full-graph training of GCN can be hampered by slow
convergence and scales poorly to very large graphs (Liu et al., 2021).

Sampling can be beneficial to both computational efficiency and scalability. A
general statistical framework for graph sampling has been developed in the
recent years (Zhang, 2021a; Zhang and Patone, 2017), including the theories
for probability-based breadth-first or depth-first sampling methods. In this
paper, we propose two sampling-based approaches to node embedding, which aims
to extract useful features from both the eigenvectors related to the graph
Laplacien and the given values associated with the graph.

In the first approach, node embeddings are given by a parametric neighbour
function that holds exactly (or approximately) for the eigenvectors of a
chosen Laplacian matrix. The function value of any given node depends directly
on those of its adjacent nodes and, thus, recursively on all the nodes
connected to it in the graph. The parameters of the neighbour function can be
defined given the whole graph and the associated values. Then, for situations
where full-graph training needs to be avoided, supervised learning from sample
graphs can yield valid estimates of the full-graph parameters.

In the second approach, nonparametric supervised embeddings are first
generated for the sampled nodes. These are then updated by repeated sampling
from the given graph, i.e. for the nodes that are sampled on each occasion.
The approach facilitates parallel computing, where many small sample graphs
can be processed separately from each other following the same procedure. In
the extreme case, the required graph computation can be minimised by letting
each node sample consist of only one selected node and its adjacent nodes. As
will be illustrated, although the expected node embedding over all possible
samples would differ somewhat to that based on the full-graph training, they
can be just as useful for node classification tasks.

##  2 Sampling-based estimation

Let $G=(U,A)$ be an undirected simple graph with $U$ as the set of nodes,
$N=|U|$, and $A$ as the adjacency matrix, where $a_{ij}=1$ if an edge exists
from $i$ to $j$ and 0 otherwise, and $a_{ij}\equiv a_{ji}$ for any $i,j\in U$.
Provided they are available, let $v_{G}=(y_{U},\omega_{A})$ contain all the
given values (or vectors) associated with the nodes and edges, where
$y_{U}=\\{y_{i}:i\in U\\}$ and $\omega_{A}=\\{\omega_{ij}:i,j\in
U,a_{ij}=1\\}$.

Let node embedding yield $x_{U}=\\{x_{i}:i\in U\\}$, where $x_{i}$ is the
feature value (or vector) of node $i$. Let $\theta$ be the unknown constants
required for node embedding, such as the parameters of a model of $x_{U}$ or
simply $\theta=x_{U}$. Given $(G,v_{G})$, let the graph-fit of $\theta$,
denoted by $\theta_{0}$, be a solution to the estimating equation

$u(\theta_{0})=\frac{\partial\Delta(\theta)}{\partial\theta}|_{\theta=\theta_{0}}=0$

where we require the loss function $\Delta(\theta)$ to be a sum over all the
nodes, i.e.

$\Delta(\theta)=\sum_{i\in U}\Delta_{i}(\theta)$

Note that $\Delta_{i}$ may depend on the node $i$ as well as its adjacent
nodes, denoted by $\nu_{i}=\\{j\in U:a_{ij}a_{ji}~=1\\}$ and referred to as
its neighbourhood in the graph, but we write simply $\Delta_{i}$ (instead of
$\Delta_{i,\nu_{i}}$) insofar as $\nu_{i}$ is fixed given $i$.

Denote by $G_{s}=(U_{s},A_{s})$ a sample graph that is selected from $G$
according to a given graph sampling method, where $U_{s}$ contains the
observed nodes of $U$ and $A_{s}$ the observed elements of the adjacency
matrix $A$. In particular, the way by which the edges drive graph sampling is
called the _observation procedure (OP)_ ; see Zhang (2021a), Zhang and Patone
(2017). For undirected simple graphs considered in this paper, we assume the
OP is incident such that, given an initial sample of nodes $s_{0}$, both
$a_{ij}$ and $a_{ji}$ are observed for any $i\in s_{0}$. Moreover, the values
in $v_{G}$ are assumed to be observed together with the associated nodes and
edges. Finally, an OP may be applied repeatedly, i.e. to the observed nodes
outside the initial $s_{0}$, and the set of nodes to which the OP is applied
is called the _seed sample_ , denoted by $s$. By the general definition of
sample graph in Zhang (2021a), $A_{s}$ contains then all the elements of $A$
which are observed from the seed sample, and $U_{s}=s\cup\\{i,j:a_{ij}\in
A_{s},a_{ij}=1\\}$.

Given $G_{s}$, we shall obtain an estimator of $\theta_{0}$, denoted by
$\hat{\theta}_{0}$, by solving a _sample estimating equation (SEE)_ , such as

$u_{s}(\hat{\theta}_{0})=\sum_{i\in U_{s}}w_{i}u_{i}(\hat{\theta}_{0})=0$

where
$u_{i}(\hat{\theta}_{0})=\partial\Delta_{i}/\partial\theta|_{\theta=\hat{\theta}_{0}}$
and $w_{i}$ is some chosen weight. In particular, the SEE is _unbiased_ over
repeated sampling if, for any $\theta_{0}$, we have

$E\\{u_{s}(\theta_{0})\\}=0$ (1)

where the expectation is over the sampling distribution of $G_{s}$ given fixed
$(G,v_{G})$.

###  2.1 Targeted random walk (TRW)

There are many kinds of random walk in graphs; see e.g. Masuda et al. (2017).
We consider the TRW (Zhang, 2021a). Starting from any node, $O_{0}\in U$, the
transition probabilities from $O_{t}=i$ to $O_{t+1}=j$ are given by

$p_{ij}=\begin{cases}\frac{1}{d_{i}+r}\big{(}1+\frac{r}{N}\big{)}&\text{if
}a_{ij}=1\\\ \frac{r}{d_{i}+r}\big{(}\frac{1}{N}\big{)}&\text{if
}a_{ij}=0\end{cases}$

where $d_{i}$ is the node degree and $r$ is a chosen tuning constant
(Avrachenkov, et al., 2010). That is, one either moves randomly over one of
the edges incident to the current node or jumps randomly to any node
(including the current one). The Markov process $\\{O_{t}:t\geq 0\\}$ is
irreducible by virtue of random jump. At equilibrium, the stationary
probability is given by

$\pi_{i}=\Pr(O_{t}=i)\propto d_{i}+r$

The same $\pi_{i}$ can as well be obtained by faster-mixing lagged random walk
(Zhang, 2021b), or one can design any individual stationary probabilities by
lagged Metropolis-Hastings walk (Zhang, 2022).

Applying the OP to the current $O_{t}=i$ enables one to observe $\nu_{i}$ and
the associated values $y_{\nu_{i}}=\\{y_{j}:j\in\nu_{i}\\}$ and
$\omega_{\nu_{i}}=\\{(\omega_{ij},\omega_{ji}):j\in\nu_{i}\\}$, if they exist.
Given any extraction of $n$ states at equilibrium, as the seed sample of TRW
and denoted by $s=\\{t_{1},...,t_{n}\\}$, an unbiased SEE (1) can be given as

$u_{s}(\theta)=\sum_{k=1}^{n}\sum_{i\in
U}\mathbb{I}(O_{t_{k}}=i)\frac{u_{i}(\theta)}{n(d_{i}+r)}$

since

$E\\{u_{s}(\theta_{0})\\}=\frac{1}{n}~\sum_{k=1}^{n}\sum_{i\in
U}\frac{\pi_{i}}{d_{i}+r}u_{i}(\theta_{0})\propto u(\theta_{0})=0$

Under TRW, the states $\\{O_{t_{1}},...,O_{t_{n}}\\}$ are correlated. The
simplest approach to unbiased variance estimation is to run multiple TRWs
independently and obtain, say, $\hat{\theta}_{0,l}$ for $l=1,...,L$. The
combined estimator and variance estimator are then given by

$\hat{\theta}_{0}=\frac{1}{L}~\sum_{l=1}^{L}\hat{\theta}_{0,l}\qquad\text{and}\qquad\hat{V}(\hat{\theta}_{0})=\frac{1}{L(L-1)}\sum_{l=1}^{L}(\hat{\theta}_{0,l}-\hat{\theta}_{0})^{2}$

###  2.2 Snowball sampling (SBS)

Let $s_{0}$ be an initial node sample, $s_{0}\subset U$. For $t=1,...,T$, let

$s_{t}=\nu(s_{t-1})\setminus\bigcup_{r=0}^{t-1}s_{r}$

be the $t$-th wave sample of nodes, which are outside of $s_{t-1}$ but
adjacent to some nodes in $s_{t-1}$. The sampling is terminated if
$s_{t}=\emptyset$ for some $t<T$, in which case $s_{r}=\emptyset$ for
$r=t,...,T$. The seed sample of $T$-wave SBS is

$s=\bigcup_{t=0}^{T-1}s_{t}$

by repeated incident OP, and the elements of $A_{s}$ are $s\times U\cup
U\times s$. See Zhang (2021a), Frank and Snijders (1994), and Goodman (1961).

Different strategies (Zhang, 2021a, Ch. 5) can be used to define the SEE. For
instance, given the sample graph $G_{s}=(U_{s},A_{s})$, let $F_{i,s}$ contain
all the nodes in $G_{s}$ that can lead to $i\in s$ by $T$-wave SBS, in
contrast to $F_{i}$ that contains all such nodes in $G$. We can let
$\delta_{i}=1$ only if $i\in s$ results from $s_{0}\cap F_{i,s}\neq\emptyset$,
and 0 otherwise — including when $i\in s$ only because of
$s_{0}\cap(F_{i}\setminus F_{i,s})\neq\emptyset$ — although the distinction
vanishes in the special case of 1-wave SBS, where

$s=s_{0}\qquad\text{and}\qquad F_{i,s}=F_{i}=\nu_{i}\cup\\{i\\}$

Next, let

$w_{i}^{-1}=\Pr(s_{0}\cap F_{i,s}\neq\emptyset)=E(\delta_{i})$

according to the design of the initial sample $s_{0}$. We have then

$E\\{u_{s}(\theta_{0})\\}=\sum_{i\in
U}E(\delta_{i})w_{i}u_{i}(\theta_{0})=u(\theta_{0})=0$

Moreover, by expansion around $\theta_{0}$, we have

$0=\sum_{i\in s}w_{i}u_{i}(\hat{\theta}_{0})\approx\sum_{i\in
s}w_{i}u_{i}(\theta_{0})+\mathcal{H}_{s}(\hat{\theta}-\theta_{0})$

where $\mathcal{H}_{s}=E\big{(}\sum_{i\in s}w_{i}\partial
u_{i}(\theta)/\partial\theta|_{\theta=\theta_{0}}\big{)}$. The sampling
variance of $\hat{\theta}_{0}$ follows as

$V(\hat{\theta}_{0})\approx\mathcal{H}_{s}^{-1}\left\\{\sum_{i,j\in
U}\Big{(}w_{i}w_{j}\Pr(\delta_{i}\delta_{j}=1)-1\Big{)}u_{i}(\theta_{0})u_{j}(\theta_{0})^{\prime}\right\\}\mathcal{H}_{s}^{-1}$

We refer to Zhang (2021a) for relevant details of T-wave SBS.

##  3 Eigen neighbour function

Let $L=D-A$ be the Laplacian matrix of $G$, where $D$ is the diagonal matrix
of degrees $d_{i}$. The normalised graph Laplacian is

$\dot{L}=D^{-\frac{1}{2}}LD^{-\frac{1}{2}}=I-D^{-\frac{1}{2}}AD^{-\frac{1}{2}}=I-M$

Denote by $\lambda_{0}$ the smallest non-zero eigenvalue of $\dot{L}$ and by
$z_{0}$ the corresponding eigenvector. The components of $z_{0}/\sqrt{d}$ are
more similar for adjacent nodes than any other eigenvector $z$ with
corresponding non-zero eigenvalue $\lambda_{z}$. Thus, $z_{0}$ may be a useful
embedding for any node values that are more similar among adjacent nodes than
otherwise, an example of which is the ZKC members split. However, one
generally cannot estimate $z_{0}$ based on a subgraph of $G$.

Regarding any binary node classification $y$, for which the eigenvector
$z_{0}$ is a good predictor, consider generalised linear node classification
in combination with eigen neighbour function (ENF) for embedding. For any
$i\in U$, let

$\begin{cases}\eta_{i}=g(p_{i})=(1,x_{i})\psi\\\
x_{i}=\xi\sum_{j\in\nu_{i}}m_{ij}y_{j}\end{cases}$ (2)

where $g(\cdot)$ is a link function mapping $p_{i}=p(x_{i};\psi)\in(0,1)$ to
the linear scale, $x_{i}$ is a scalar feature (or embedding) of node $i$ which
depends on its adjacent nodes, $m_{ij}$ is an element of the normalised
adjacency matrix $M$, $\xi$ is a scalar parameter for embedding and $\psi$ is
a (column) vector parameter for classification given $x_{i}$. Both $\xi$ and
$\psi$ are parameters that need to be estimated.

#### Classification

One may consider $p_{i}$ as a probability for binary $y_{i}$, according to
which a classifier can simply be $\mathbb{I}(p_{i}>0.5)$. Given $(G,x_{U})$,
the graph-fit $\psi_{0}$ can be obtained by the loss and estimation equation
below

$\displaystyle\Delta(\psi)\propto\sum_{i\in
U}y_{i}\log(p_{i})+(1-y_{i})\log(1-p_{i})$ $\displaystyle
u(\psi)=\frac{\partial\Delta(\psi)}{\partial\psi}~=\sum_{i\in
U}\big{(}\frac{\partial
p_{i}}{\partial\psi}\big{)}\frac{y_{i}-p_{i}}{p_{i}(1-p_{i})}=0$

For example, the logistic link yields $p_{i}=1/\big{(}1+e^{-\eta_{i}}\big{)}$,
such that

$u(\psi)=\sum_{i\in
U}(y_{i}-p_{i})(1,x_{i})^{\prime}\quad\text{and}\quad\frac{\partial
u(\psi)}{\partial\psi}=-\sum_{i\in
U}p_{i}(1-p_{i})(1,x_{i})^{\prime}(1,x_{i})$

The tanh link yields $p_{i}=\frac{1}{2}(1+t_{i})$ where
$t_{i}=\big{(}e^{\eta_{i}}-e^{-\eta_{i}}\big{)}/\big{(}e^{\eta_{i}}+e^{-\eta_{i}}\big{)}$,
such that

$u(\psi)=\sum_{i\in
U}(y_{i}-p_{i})(1,x_{i})^{\prime}\quad\text{and}\quad\frac{\partial
u(\psi)}{\partial\psi}=-\frac{1}{2}~\sum_{i\in
U}(1+t_{i})(1-t_{i})(1,x_{i})^{\prime}(1,x_{i})$

#### Embedding

The underlying neighbour function for $x_{i}$ in (2) can be written as

$\mu=\xi\dot{\mu}~\qquad\text{and}\qquad\dot{\mu}~=M\mu$

in terms of an $N$-vector $\mu^{\prime}=(\mu_{1},...,\mu_{N})$. Since any
eigenvector $z$ of $\dot{L}$ satisfies $z=(1-\lambda_{z})^{-1}Mz$, the feature
vector $x=\xi My$ would hopefully combine useful information in the
eigenvectors of $\dot{L}$ and the given labels $y$; hence the term ENF. Let
$e=y-\mu$ be the discrepancy of $\mu$ satisfying $\mu=\xi_{0}\dot{\mu}$, such
that

$y-\xi_{0}\dot{y}=(\mu+e)-\xi_{0}(\dot{\mu}+\dot{e})=e-\xi_{0}\dot{e}$

Let the total loss be

$\Delta(\xi)=\sum_{i\in U}\Delta_{i}(\xi)=\sum_{i\in
U}(y_{i}-\xi\dot{y}_{i})^{2}$ (3)

Given $(G,y_{U})$, the graph-fit $\xi_{0}$ minimising $\Delta(\xi)$ is given
by

$\xi_{0}=\big{(}\dot{y}^{\prime}~\dot{y}~\big{)}^{-1}\big{(}\dot{y}^{\prime}y\big{)}$

#### Supervised learning

Estimates of the graph-fit $(\xi_{0},\psi_{0})$ can be obtained from the
respective SEEs given any sample graph $G_{s}$ and associated $y$-values. In
particular, under 1-wave SBS with seed sample $s=s_{0}$, yielding
$U_{s}=\bigcup_{i\in s}\nu_{i}\cup s$ and the associated $y_{U_{s}}$, let the
SEE for $\xi$ be

$u_{s}(\xi)=\frac{\partial}{\partial\xi}\sum_{i\in
s}\pi_{i}^{-1}(y_{i}-\xi\dot{y}_{i})^{2}\propto\sum_{i\in
s}\pi_{i}^{-1}\dot{y}_{i}(y_{i}-\xi\dot{y}_{i})=0$

such that

$\hat{\xi}_{0}=\Big{(}\sum_{i\in
s}\frac{\dot{y}_{i}^{2}}{\pi_{i}}\Big{)}^{-1}\Big{(}\sum_{i\in
s}\frac{\dot{y}_{i}y_{i}}{\pi_{i}}\Big{)}$

Moreover, let $\hat{x}_{i}=\hat{\xi}_{0}\dot{y}_{i}$ and obtain
$\hat{\psi}_{0}$ from the plug-in SEE for $\psi$ given by

$u_{s}(\psi)=\sum_{i\in s}\pi_{i}^{-1}(y_{i}-p_{i})(1,\hat{x}_{i})^{\prime}=0$

###  3.1 Example: ZKC

Denote the members split shown in Figure 1 by a binary label $y_{i}$, for all
$i\in U$. As mentioned before, given the eigenvector $z_{0}$, a perfect
classifier is $y_{i}=\mathbb{I}(z_{0,i}>\tau)$ for any $i\in U$, given any
parameter value $\tau\in(0.0528,0.0556)$.

Figure 2: Normalised $x=\xi_{0}\dot{y}$ against eigenvector $z_{0}$ (top-left)
or $E(\hat{x})$ (top-right), given $\hat{x}=\hat{\xi}_{0}\dot{y}$ and
$\hat{\xi}_{0}$ by snowball sampling with $|s|=5$. Logistic-ENF
$p(x;\psi_{0})$ against $z_{0}$ (bottom-left) or tanh-ENF $p$ (bottom-right).
Members split marked.

The top-left plot of Figure 2 shows the normalised ENF embedding
$x=\xi_{0}\dot{y}$ given the graph-fit $\xi_{0}=0.955$, where $x^{\prime}x=1$.
The ENF yields a feature vector $x$ that is similar to $z_{0}$, but less
discriminant than $z_{0}$ for the $y$-classification, since there clearly
exists an interval of $x$ in which $y$ can be both 0 and 1. The top-right plot
of Figure 2 shows $E(\hat{x})$ against $x$, where
$\hat{x}=\hat{\xi}_{0}\dot{y}$ and $\hat{\xi}_{0}$ is obtained by 1-wave SBS
with initial simple random sampling (SRS) of sample size $|s|=5$. The sample-
graph estimator $\hat{\xi}_{0}$ is unbiased for $\xi_{0}$ and embedding $x$.

The bottom-left plot of Figure 2 shows the logistic-ENF
$p(\xi_{0}\dot{y};\psi_{0})$ against the eigenvector $z_{0}$ given the graph-
fit $\psi_{0}=(-4.631,15.747)$. Applying $\mathbb{I}(p_{i}>0.5)$ to all the 34
nodes would incorrectly classify two 1-nodes and two 0-nodes. The bottom-right
plot of Figure 2 shows the logistic-ENF against the tanh-ENF, given the graph-
fit $\psi_{0}=(-2.315,7.874)$ for the latter, which shows that the choice of
link does not matter for this valued graph.

###  3.2 Discussion

In the statistical literature of spatial and network analysis, e.g. Ord
(1975), Friedkin (1990) and Leenders (2002), the linear autoregressive
predictor

$\mu=\rho H\mu+c\qquad\Leftrightarrow\qquad\mu=(I-\rho H)^{-1}c$

builds on node features $c_{i}=y_{i}^{\prime}\beta$ for $i\in U$, where
$\beta$ is a parameter vector, $H=[h_{ij}]$ and $h_{ij}=0$ for any
$j\notin\nu_{i}$. Provided $\|\rho H\|<1$, we have

$\mu=(I+\rho H+\rho^{2}H^{2}+\cdots)c$

showing how $\mu_{i}$ depends on $c_{i}$, those of its 1-hop neighbours, 2-hop
neighbours, etc. For $l=1,...,L$, the $l$-th order approximation is given by
iterating

$\mu^{(l)}=\rho H\mu^{(l-1)}+c\quad\text{and}\quad\mu^{(0)}=c$

Such network influencing has certain similarity to messaging-and-aggregating
by an $L$-layer GNN (e.g. Hamilton, 2020), e.g.

$\mu^{(l)}=\sigma\Big{(}\rho_{l}H\mu^{(l-1)}+\phi_{l}\mu^{(l-1)}\Big{)}$

where $\sigma(\cdot)$ is an activation function such as $\max(0,\cdot)$, and
$\\{(\rho_{l},\phi_{l}):1\leq l\leq L\\}$ are the parameters to be learned.
The initial $\mu^{(0)}=c$ can be directly given by known node features or
functions of them.

The ENF $\dot{x}=My$ is a one-layer GNN with identity activation function
$\sigma(\cdot)$, where $x^{(0)}=y$, $H=M$, $\xi=\rho_{1}$ and $\phi_{1}=0$.
The sampling-based training approach to ENF above can be adapted to other
model-based embeddings. For instance, let $x^{(0)}=y$, and let $x$ be given by
a 2-layer GNN as

$x=\sigma\Big{(}\rho_{2}H\sigma\big{(}\rho_{1}Hy+\phi_{1}y\big{)}+\phi_{2}\sigma\big{(}\rho_{1}Hy+\phi_{1}y\big{)}\Big{)}=\mathring{y}$

with parameters $\xi=(\rho_{2},\rho_{2},\phi_{2},\phi_{1})$. Similarly to (3),
one can use the loss

$\Delta_{i}(\xi)=(y_{i}-\mathring{y}_{i})^{2}$

for the underlying 2-hop neighbour function satisfying
$\mu_{i}=\mathring{\mu}_{i}(\xi_{0})$. For sampling-based estimation of
$\xi_{0}$, one now needs to observe $y_{j}$ of any node $j$ that belongs to
the 2-hop neighbourhood of a seed sample node $i$. In the case of $s=s_{0}$,
one can implement the 2-hop incident OP, which yields $A_{s}=\tilde{s}\times
U\cup U\times\tilde{s}$ where $\tilde{s}~=\bigcup_{i\in s}\nu_{i}\cup s$, and
$U_{s}=\bigcup_{i\in s}\nu_{i}^{2}\cup s$ where
$\nu_{i}^{2}=\bigcup_{j\in\nu_{i}}\nu_{j}$, and the associated $y_{U_{s}}$.

##  4 Normalised Laplacian embedding

Consider a different approach to eigenvector-like embeddings, which operates
without an explicit model but through an appropriately chosen linear system
$Pz=0$, which is satisfied exactly (or nearly so) by any eigenvector $z$ of
$\dot{L}$. In particular, for any pair of eigenvalue and eigenvector
$(\lambda_{z},z)$, we can let

$P_{\lambda_{z}}(M)=(1-\lambda_{z})I-M$

or

$P_{\lambda_{z}}(\tilde{M})=\mbox{Diag}\Big{(}1-\frac{\lambda_{z}d_{i}}{1+d_{i}}\Big{)}-\tilde{M}$

where

$\tilde{M}=(I+D)^{-\frac{1}{2}}(I+A)(I+D)^{-\frac{1}{2}}$

It is convenient to define $\tilde{\nu}_{i}=\\{j\in U:p_{ij}\neq 0\\}$ as the
looped neighbourhood of node $i$, to be distinguished from $\nu_{i}$, since
$P$ has non-zero diagonal elements. Note that both $P_{\lambda}(M)$ and
$P_{\lambda}(\tilde{M})$ are symmetric matrices.

Given any tuning constant $\lambda\in(0,2)$, let the corresponding ($y$-)
supervised normalised Laplacian embedding (SNLE) be given by minimising

$~\Delta(x)=x^{\prime}P_{\lambda}^{\prime}P_{\lambda}x+\gamma(y-x)^{\prime}(y-x)$
(4)

over $x$, where $P_{\lambda}$ is given by replacing $\lambda_{z}$ with
$\lambda$ in $P_{\lambda_{z}}$, and $\gamma$ is a regularisation constant.
Setting $\partial\Delta(x)/\partial x$ to $0$ yields the graph-fit

$x_{0}=\big{(}I+\gamma^{-1}P_{\lambda}^{\prime}P_{\lambda}\big{)}^{-1}y$

Any $x$ satisfying $x^{\prime}P_{\lambda}^{\prime}P_{\lambda}x\approx 0$ would
entail $P_{\lambda}x\approx\bm{0}$, analogously to $P_{\lambda_{z}}z=0$. As we
illustrate below, although the loss(4) is defined whether or not $\lambda$ is
an eigenvalue of $\dot{L}$, it can generate embeddings that are correlated
with any eigenvector of the normalised Laplacian of the graph — hence the term
SNLE.

###  4.1 Example: ZKC

We have $\lambda_{0}=0.132$ for the ZKC graph (Figure 1). As illustrated in
Figure 3, full-graph embedding $x_{0}$ can closely resemble $z_{0}$ for a
range of $\lambda$, given $\gamma=0$.

Figure 3: Eigenvector $z_{0}$ vs. normalised $x_{0}$ given $\lambda=0.1$
(left) or $0.15$ (right) Figure 4: Ranks of best-correlated eigenvector with
SNLE $x_{0}$ as $\lambda$ varies, with fixed $\gamma\equiv 0.001$. Left,
$P_{\lambda}(M)$; right, $P_{\lambda}(\tilde{M})$.

Figure 4 shows the rank of the eigenvector that is best correlated with
$x_{0}$ given by (4), as $\lambda$ varies from 0.1 to 1.9 while $\gamma$ is
fixed at 0.0001. The eigenvectors are ranked from 1 to 34 in the decreasing
order of their corresponding eigenvalues, where the rank of $z_{0}$ is 33
(marked top horizontal dashed line). It is seen that SNLE (4) can yield
different eigenvector-like embeddings as $\lambda$ varies. Note that the non-
monotonicity seen here for $P_{\lambda}(\tilde{M})$ can be removed by reducing
$\gamma$, e.g. to 0.0001, the details of which are omitted here.

###  4.2 Sample-graph embedding

Rewrite the loss $\Delta(x)$ in (4) as a sum over all the nodes as

$\Delta(x)=\sum_{i\in U}\dot{x}_{i}^{2}+\gamma\sum_{i\in U}(y_{i}-x_{i})^{2}$

where $\dot{x}_{i}=\sum_{j\in\tilde{\nu}_{i}}p_{ij}x_{j}$. By 1-wave SBS with
seed sample $s=s_{0}$, one observes $U_{s}=\bigcup_{i\in s}\nu_{i}\cup
s=\bigcup_{k\in s}\tilde{\nu}_{k}$ and the associated $y_{U_{s}}$, as well as
the sub-matrices $A_{s}=A_{sU_{s}}\cup A_{U_{s}s}$, where
$A_{sU_{s}}=[a_{kj}]_{n\times n_{U_{s}}}$ for $k\in s$ and $j\in U_{s}$ and
$n=|s|$ and $n_{U_{s}}=|U_{s}|$. Given any fixed $N$-vector $x$, an unbiased
estimator of $\Delta(x)$ over repeated sampling of $G_{s}=(U_{s},A_{s})$ is
given by

$\Delta_{s}=\sum_{k\in s}\pi_{k}^{-1}\dot{x}_{k}^{2}+\gamma\sum_{j\in
U_{s}}\dot{\pi}_{j}^{-1}(y_{j}-x_{j})^{2}$

where $\pi_{k}=\Pr(i\in s)$, and $\dot{\pi}_{j}=\Pr(j\in U_{s})$. For
instance, under SRS of $s$, we have $\pi_{k}=n/N$ and
$\dot{\pi}_{j}=1-\binom{N-d_{j}-1}{n}/\binom{N}{n}$.

Given any $(G_{s},y_{U_{s}})$, supervised learning is feasible for $j\in
U_{s}$. We have

$\frac{1}{2}~\frac{\partial\Delta_{s}}{\partial x_{j}}=\sum_{k\in
s\cap\tilde{\nu}_{j}}\frac{p_{jk}}{\pi_{k}}\sum_{l\in\tilde{\nu}_{k}}p_{kl}x_{l}-\frac{\gamma}{\dot{\pi}}(y_{j}-x_{j})$

where $\partial\dot{x}_{k}/\partial x_{j}=p_{jk}$ if $k\in
s\cap\tilde{\nu}_{j}$ given $k\in s$, and 0 otherwise. Let $y_{U_{s}}$ and
$x_{U_{s}}$ be vectors associated with $U_{s}$. Let
$P_{sU_{s}}=[p_{kj}]_{n\times n_{U_{s}}}$ for $k\in s$ and $j\in U_{s}$.
Define the $n\times n$-matrix $W_{s}=\mbox{Diag}(\pi_{i}^{-1})$ associated
with $s$, and the $n_{U_{s}}\times n_{U_{s}}$-matrix
$W_{U_{s}}=\mbox{Diag}(\dot{\pi}_{j}^{-1})$ associated with $U_{s}$. Let the
SEE be given by

$u_{s}=\frac{\partial\Delta_{s}}{\partial
x_{U_{s}}}=P_{sU_{s}}^{\prime}W_{s}P_{sU_{s}}x_{U_{s}}-\gamma
W_{U_{s}}(y_{U_{s}}-x_{U_{s}})=0$

Let the solution to the SEE be

$\hat{x}_{U_{s}}=\big{(}\gamma^{-1}W_{U_{s}}^{-1}P_{sU_{s}}^{\prime}W_{s}P_{sU_{s}}+I\big{)}^{-1}y_{U_{s}}$

For any $i\in U$, let the sample-graph SNLE $x_{i}$ over repeated sampling to
be

$x_{i}=E(\hat{x}_{i}\mid i\in U_{s})$ (5)

Whereas the sample-graph embedding (5) surely converges to the full-graph
embedding (4) as $|s|\rightarrow|U|$, the two are unequal generally, unless
$\tilde{\nu}_{j}\equiv s\cap\tilde{\nu}_{j}$ for any $j\in U_{s}$ under graph
sampling. The reason is that $\partial\Delta/\partial x_{j}$ when $\nu_{j}$ is
fully observed may differ to $\partial\Delta_{s}/\partial x_{j}$, for any
$j\in U_{s}\setminus s$, i.e.

$\frac{\partial\Delta}{\partial
x_{j}}=2\sum_{k\in\tilde{\nu}_{j}}p_{jk}\dot{x}_{k}\neq\frac{\partial\Delta_{s}}{\partial
x_{j}}=2\sum_{k\in s\cap\tilde{\nu}_{j}}p_{jk}\dot{x}_{k}$

Thus, derivation and expectation are not exchangeable,

$E(\partial\Delta_{s}/\partial x_{j})\neq\partial E(\Delta_{s})/\partial
x_{j}$

The matter is the same for any graph sampling method, as long as there is a
distinction between a seed sample node (whose neighbourhood is observed by
definition of incident OP) and another sample node in $U_{s}$ (whose
neighbourhood is not always observed when sampling is terminated).

Figure 5: SNLE $x_{0}(G,y_{U})$ vs. sample embedding $E(\hat{x})$ given
$\gamma=0.1$ (left) or $\gamma=1$ (right), using $\lambda=0.1$,
$P_{\lambda}(\tilde{M})$, seed sample size $|s|=1$. Members split marked.

Nevertheless, Figure 5 compares full-graph embedding (4) for the ZKC data to
sample-graph embedding by (5), which is obtained by 1-wave SBS where the
initial sample consists of a single randomly selected node. Both embeddings
are tuned towards $z_{0}$ here given $P_{0.1}(\tilde{M})$. Using $|s|=1$ that
minimises the required graph computation, we obtain the sample-graph embedding
$E(\hat{x})$ with $\gamma=0.1$ (left) or $\gamma=1$ (right), which readily
yields $p(E(\hat{x});\psi)$ that can be more discriminating than
$p(x_{0};\psi)$ or $p(z_{0};\psi)$ for $y$-like labels, now that $E(\hat{x})$
separates the 1-nodes and 0-nodes by a wider margin in both plots.

## References

  * [1] Avrachenkov, K., Ribeiro, B. and Towsley, D. (2010). Improving Random Walk Estimation Accuracy with Uniform Restarts. _Research report_ , RR-7394, INRIA. inria-00520350
  * [2] Belkin, M. and Niyogi, P. (2001). Laplacian eigenmaps and spectral techniques for embedding and clustering. Advances in Neural Information Processing Systems, 14:586-691.
  * [3] Brandes, U., Delling, D., Gaertler, M., Gorke, R., Hoefer, M., Nikoloski, Z. and Wagner, D. (2008). On modularity clustering. IEEE Transactions on Knowledge and Data Engineering, 20(2):172-188.
  * [4] Frank O. and Snijders T. (1994). Estimating the size of hidden populations using snowball sampling. Journal of Official Statistics, 10:53-53.
  * [5] Friedkin, N.E. (1990). Social Networks in Structural Equation Models. Social Psychology Quarterly, 53:316-328.
  * [6] Goodman, L.A. (1961). Snowball sampling. Annals of Mathematical Statistics, 32:148-170.
  * [7] Hamilton, W.L. (2020). Graph Representation Learning. Synthesis Lectures on Artificial Intelligence and Machine Learning, Vol. 14, No. 3 , pp. 1-159. Morgan and Claypool.
  * [8] Khoshraftar, S. and An, A. (2022). A Survey on graph representation learning methods. https://doi.org/10.48550/arXiv.2204.01855
  * [9] Kipf, T. and Welling, M. (2017). Semi-supervised classification with graph convolutional networks. https://doi.org/10.48550/arXiv.1609.02907
  * [10] Leenders, R. (2002). Modeling social influence through network autocorrelation: constructing the weight matrix. Social Networks, 24:21-47.
  * [11] Masuda, N., Porter, M.A. and Lambiotte, R. (2017) Random walks and diffusion on networks. _Physics Reports_ , 716-717:1-58.
  * [12] Ord, K. (1975). Estimation Methods for Models of Spatial Interaction. Journal of the American Statistical Association, 70:120-126.
  * [13] Perozzi, B., Al-Rfou, R. and Skiena, S. (2014). Deepwalk: Online learning of social representations. https://doi.org/10.48550/arXiv.1403.6652
  * [14] Zachary, W. W. (1977). An information flow model for conflict and fission in small groups. Journal of Anthropological Research, 33:452-473.
  * [15] Zhang, L.-C. (2022). Graph spatial sampling. https://doi.org/10.48550/arXiv.2203.11734
  * [16] Zhang, L.-C. (2021a). _Graph sampling._ CRC Press.
  * [17] Zhang, L.-C. (2021b). Graph sampling by lagged random walks. Stat, https://onlinelibrary.wiley.com/doi/abs/10.1002/sta4.444
  * [18] Zhang, L.-C. and Patone, M. (2017). Graph sampling. _Metron_ , 75:277-299. https://doi.org/10.1007/s40300-017-0126-y
  * [19] Zhou, J., Cui, G., Hu, S., Zhang, Z., Yang, C., Liu, Z., Wang, L., Li, C. and Sun, M. (2020). Graph neural networks: A review of methods and applications. AI Open, 1: 57-81.


# Minimisation of Polyak-Łojasewicz Functions Using Random Zeroth-Order
Oracles

Amir Ali Farzin and Iman Shames The authors are with CIICADA Lab, School of
Engineering, ANU<EMAIL_ADDRESS>
###### Abstract

The application of a zeroth-order scheme for minimising Polyak-Łojasewicz (PL)
functions is considered. The framework is based on exploiting a random oracle
to estimate the function gradient. The convergence of the algorithm to a
global minimum in the unconstrained case and to a neighbourhood of the global
minimum in the constrained case along with their corresponding complexity
bounds are presented. The theoretical results are demonstrated via numerical
examples.

##  I Introduction

Zeroth-order (ZO) or derivative-free optimisation schemes are of interest when
the gradient (or subgradient in case of non-differentiable cost functions)
information is not readily available. A common scenario is the case where the
value of the cost function, and not its higher order derivatives, is the only
information available to the solver [1] and [2]. Sometimes, even if the
gradient value is theoretically available, evaluating the gradient might incur
high computational costs [3]. ZO methods provide a way forward for solving
such optimisation problems as well. The majority of existing ZO methods aim to
construct an estimate of the gradient of the function and use this estimate as
a surrogate for the gradient. The method analysed in this paper is no
exception to this general trend.

Various ZO optimisation methods have been designed and analysed for different
problem classes. In [4], a constrained stochastic composite optimisation
problem was studied where the function is possibly non-convex. The proposed
algorithm in [4] relied on the existence of an unbiased variance-bounded
estimator of the gradient. In [5], the authors considered the unconstrained
zeroth-order optimisation problem and defined a random oracle which was an
unbiased variance-bounded estimator of the gradient of a smoothed version of
the original function.

The Polyak-Łojasiewicz (PL) inequality was first introduced in [6] and the
convergence of gradient descent method under PL assumption was first proved
there. It is observed that a range of different cost functions satisfy the PL
condition [7]. Karimi et al. used the PL inequality to provide a new proof
technique for analysing various first-order gradient descent methods and used
proximal PL condition to analyse non-smooth cases [8]. In [9], a variant of
the direct search method was employed to solve stochastic minimisation and
saddle point problems. The direct search algorithm which is a derivative-free
scheme and obtained the complexity bounds for the convergence of PL functions.

In this paper, we aim to fill a gap in the existing literature on random
zeroth order oracles. We specifically leverage the class of random method
proposed in [5] for optimisation problems, and establish its performance for
the case where the cost functions satisfy the (proximal) PL condition.
Specifically, we establish the convergence properties and the complexity
bounds of these methods for both constrained and unconstrained cases. By doing
so, we hope to shed some light on the behaviour of such algorithms when
applied to a class of benignly nonconvex problems.

The outline of this paper is as follows. In Section II we introduce the
necessary background information. The problems of interest are outlined in
III. The convergence properties and the complexity of the algorithms for
solving the problems of interest are presented in Section IV. Numerical
examples are presented in V and conclusions and possible future research
directions come in the end.

##  II Preliminaries

In this section we provide the necessary definitions and background material
required for presenting the results of this paper.

Consider a function $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$. The gaussian
smoothed version of $f$, termed $f_{\mu}(x)$, is defined below:

$\displaystyle\begin{split}&f_{\mu}(x)\overset{\mathrm{def}}{=}\frac{1}{\kappa}\underset{E}{\int}f(x+\mu
u)e^{-\frac{1}{2}||u||^{2}}\mathrm{d}u,\\\
&\kappa\overset{\mathrm{def}}{=}\underset{E}{\int}e^{-\frac{1}{2}||u||^{2}}\mathrm{d}u=\frac{(2\pi)^{n/2}}{[detB]^{\frac{1}{2}}},\end{split}$
(1)

where vector $u\in R^{n}$ is sampled from zero mean Gaussian distribution with
positive definite correlation operator $B^{-1}$ and the positive scalar
smoothing parameter $\mu$. Define the random oracle $g_{\mu}$ as

$g_{\mu}(x)\overset{\mathrm{def}}{=}\frac{f(x+\mu u)-f(x)}{\mu}Bu,$ (2)

where $u$ and $B$ are defined above. The projection operator on a convex set
$\mathcal{Z}$ is defined as

$\mathrm{Proj}_{\mathcal{Z}}(x)\overset{\text{def}}{=}\arg\underset{z\in\mathcal{Z}}{\min}\|z-x\|^{2},$

where $\|\cdot\|$ is the Euclidean norm of its argument if it is a vector, and
the corresponding induced norm if the argument is a matrix.

######  Definition 1 ($C^{1,1}$ Functions).

The differentiable function $f(x):D\rightarrow\mathbb{R}$ with
$D\subseteq\mathbb{R}^{n}$ as its domain is in $C^{1,1}$ if its gradient is
Lipschitz, i.e.,

$||\nabla f(x)-\nabla f(y)||\leq L_{1}(f)||x-y||,\forall x,y\in D,$ (3)

where $L_{1}(f)>0$ is the _gradient Lipschitz constant_.

######  Remark 1.

Any $f(x)\in C^{1,1}$ satisfies

$f(y)\leq f(x)+\langle\nabla f(x),y-x\rangle+\frac{L_{1}(f)}{2}||y-x||^{2}.$
(4)

######  Definition 2 (PL Functions [6]).

The differentiable function $f(x):D\rightarrow\mathbb{R}$ with
$D\subseteq\mathbb{R}^{n}$ as its domain is termed a _PL function_ if it
satisfies the Polyak-Łojasewicz (PL) condition, i.e.,

$\dfrac{1}{2}||\nabla f(x)||^{2}\geq l(f(x)-f^{\ast}),\;\forall x\in D,$ (5)

where $l>0$ is the _PL constant_ and $f^{\ast}=f(x^{\ast}).$

We have the following result for PL functions.

######  Lemma 1.

PL inequality implies that every stationary point of the function is a global
minimum.

###### Proof.

Assume that $x$ is an arbitrary stationary point, i.e., $||\nabla f(x)||=0$.
Substituting $x$ in (5) yields

$f(x)-f^{\ast}\leq 0\implies f(x)=f^{\ast}.$

Hence, any stationary point corresponds to the global minimum. ∎

######  Definition 3 (Proximal PL Functions).

Consider the function $F(x)=f(x)+h(x)$ where $f(x):D\rightarrow\mathbb{R}$
with $D\subseteq\mathbb{R}^{n}$ is a differentiable function over as its
domain and $h(x)$ is possibly nondiffrentiable. The function $F(x)$ is a
proximal PL (PPL) function if it satisfies the PPL condition on a set
$\mathcal{X}$ with positive constant $l$, for all $x\in\mathcal{X}$, i.e.,

$\displaystyle\frac{1}{2}Q(x,L_{1}(f))\geq l(F(x)-F(x^{\ast})),$ (6)

where

$\displaystyle\begin{split}Q(x,a)\overset{\mathrm{def}}{=}&-2a\underset{z\in\mathcal{X}}{\min}\\{\frac{a}{2}||z-x||^{2}+\\\
&\langle\nabla f(x),z-x\rangle+h(z)-h(x)\\}\end{split}$ (7)

with $a$ being a positive scalar.

The Proximal PL condition is a generalisation of the PL condition (5). To see
more examples, refer to [8].

To state complexity results we use the big-O notation in the sense defined
below.

######  Definition 4 (The big O-notation).

Suppose $f(x)$ and $g(x)$ are two positive scalar functions defined on some
subset of the real numbers. We write $f(x)=\mathcal{O}(g(x))$, and say $f(x)$
is in the order of $g(x)$, if and only if there exist constants $\bar{K}$ and
$M$ such that $f(x)\leq Mg(x)$ for all $x>K$.

##  III Problems of Interest

In this paper we study the performance of a random zeroth-order method for
optimising PL functions for both unconstrained and constrained cases.
Specifically, first, we study the following unconstrained problem

$\min_{x\in R^{n}}\quad f(x)$ (8)

where $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is in $C^{1,1}$, satisfies the
PL condition (5), and is bounded below.

Later, we will focus on the following constrained problem:

$\min_{x\in\mathcal{X}}\quad f(x),$ (9)

where $\mathcal{X}$ is a compact convex set in $\mathbb{R}^{n}$ with diameter
$d_{x}$ and $f(x)$ is in $C^{1,1}$ and
$F(x)=f(x)+\mathcal{I}_{\mathcal{X}}(x)$ satisfies (6) where
$\mathcal{I}_{\mathcal{X}}(x)$ is the indicator function of set $\mathcal{X}$,
i.e.,

$\mathcal{I}_{\mathcal{X}}(x)=\begin{cases}0&x\in\mathcal{X}\\\
\infty&x\not\in\mathcal{X}.\end{cases}$

The boundedness of $\mathcal{X}$ and continuity of $f(x)$ guarantees that the
problem has a solution.

##  IV Main Results

In this section, we establish the convergence properties and the complexity
bounds of a well-known class of random zeroth-order algorithms proposed by
Nesterov and Spokoiny in [5] for minimising (proximal) PL functions. We
consider both unconstrained and constrained cases.

###  IV-A Unconstrained Problem

In this subsection, we study the unconstrained problem (8). The framework
introduced in [5] is recalled in Algorithm 1, where $x_{0}$ is the initial
guess, $h_{k}$ is the step size and $N$ is the number of iterations.

Algorithm 1 $\mathrm{RS}_{\mu}$

1:Input: $x_{0},\mu,h_{k},N$

2:for $k=0,\dots,N$ do

3: Generate $u$

4: Calculate $g_{\mu}(x_{k})$ using (2)

5: $x_{k+1}=x_{k}-h_{k}g_{\mu}(x_{k})$

6:end for

7:return $x_{N}$

The following theorem characterises the convergence properties of Algorithm 1
applied to problem (8).

######  Theorem 1.

Let the sequence $\\{x_{k}\\}_{k\geq 0}$ be generated by Algorithm 1 (RSμ),
where $f(x)\in C^{1,1}$ and satisfies PL condition. Then, for any $N\geq 0$,
with $h_{k}=\frac{1}{4(n+4)L_{1}(f)}$, we have

$\displaystyle\begin{split}\frac{1}{N+1}&\sum_{k=0}^{N}(\Phi_{k}-f^{\ast})\leq\frac{8(n+4)L_{1}(f)}{l}\left[\frac{f(x_{0})-f^{\ast}}{N+1}\right.\\\
&+\left.\frac{3\mu^{2}(n+4)}{32}L_{1}(f)\right]+\frac{\mu^{2}}{4l}L_{1}^{2}(f)(n+6)^{3},\end{split}$
(10)

where $\Phi_{k}\overset{def}{=}E_{\mathcal{U}_{k}}[f(x_{k})]$ for all $k\geq
1$ and $\Phi_{0}=f(x_{0})$. Also,
$\mathcal{U}_{k}=\\{u_{0},u_{1},\dots,u_{k-1}\\}$ and $n$ is the dimension of
$x$.

###### Proof.

From [5, Equation (12)], we know that if $f\in C^{1,1}$, then $f_{\mu}\in
C^{1,1}$ and $L_{1}(f_{\mu})\leq L_{1}(f)$. Thus writing (4) for $f_{\mu}(x)$
at points $x_{k}$ and $x_{k+1}$

$\displaystyle\begin{split}f_{\mu}(x_{k+1})\leq&f_{\mu}(x_{k})-h_{k}\langle\nabla
f_{\mu}(x_{k}),g_{\mu}(x_{k})\rangle\\\
&+\frac{1}{2}h_{k}^{2}L_{1}(f_{\mu})||g_{\mu}(x_{k})||^{2}.\end{split}$ (11)

From [5, Equation (21)], we know for a function $f(x)$ we have

$\nabla f_{\mu}(x)=\frac{1}{\kappa}\int\frac{f(x+\mu
u)-f(x)}{\mu}e^{-\frac{1}{2}||u||^{2}}Bu\mathrm{d}u.$

From the definition of $g_{\mu}(x)$ and the term obtained for $\nabla
f_{\mu}(x)$, we have $E_{u}[g_{\mu}(x)]=\nabla f_{\mu}(x)$ (it means
$g_{\mu}(x)$ is an unbiased estimator of $\nabla f_{\mu}(x)$). Taking the
expectation in $u_{k}$ yields

$\displaystyle\begin{split}E_{u_{k}}(f_{\mu}(x_{k+1}))\leq&f_{\mu}(x_{k})-h_{k}||\nabla
f_{\mu}(x_{k})||^{2}\\\
&+\frac{1}{2}h_{k}^{2}L_{1}(f_{\mu})E_{u_{k}}(||g_{\mu}(x_{k})||^{2}).\end{split}$
(12)

For a function $f\in C^{1,1}$ from [5, Lemma 5], we have

$\displaystyle\begin{split}E_{u}[||g_{\mu}(x)||^{2}]\leq&4(n+4)||\nabla
f_{\mu}(x)||^{2}\\\ &+3\mu^{2}L_{1}^{2}(f)(n+4)^{3}\end{split}$ (13)

Substituting (13) in (12) and noting $L_{1}(f_{\mu})\leq L_{1}(f)$, we obtain

$\displaystyle\begin{split}E_{u_{k}}(f_{\mu}(x_{k+1}))&\leq
f_{\mu}(x_{k})-h_{k}||\nabla f_{\mu}(x_{k})||^{2}+\\\
&\frac{1}{2}h_{k}^{2}L_{1}(f)\left(4(n+4)||\nabla
f_{\mu}(x_{k})||^{2}\right.\\\
&\left.+3\mu^{2}L_{1}^{2}(f)(n+4)^{3}\right).\end{split}$ (14)

Fixing $h_{k}=\hat{h}\overset{def}{=}\frac{1}{4(n+4)L_{1}(f)}$:

$\displaystyle\begin{split}E_{u_{k}}(f_{\mu}(x_{k+1}))&\leq
f_{\mu}(x_{k})-\frac{1}{2}\hat{h}||\nabla f_{\mu}(x_{k})||^{2}\\\
&+\frac{3\mu^{2}}{32}L_{1}(f)(n+4).\end{split}$ (15)

Taking the expectation of this inequality with respect to
$\mathcal{U}_{k-1}=\\{u_{0},u_{1},\dots,u_{k-1}\\}$, we obtain

$\Phi_{k+1}\leq\Phi_{k}-\frac{1}{2}\hat{h}\Xi_{k}^{2}+\frac{3\mu^{2}(n+4)}{32}L_{1}(f),$
(16)

where $\Xi_{k}^{2}=E_{\mathcal{U}_{k}}(||\nabla f_{\mu}(x_{k})||^{2})$.
Considering $f^{\ast}\leq f(x_{N+1})$ and summing (16) from $k=0$ to $k=N$ and
divide it by $N+1$, we get

$\displaystyle\begin{split}\frac{1}{N+1}\sum_{k=0}^{N}\Xi_{k}^{2}&\leq
8(n+4)L_{1}(f)[\frac{f(x_{0})-f^{\ast}}{N+1}\\\
&+\frac{3\mu^{2}(n+4)}{32}L_{1}(f)],\end{split}$ (17)

From [5, Lemma 4], it is known that for a function $f\in C^{1,1}$

$||\nabla f(x)||^{2}\leq 2||\nabla
f_{\mu}(x)||^{2}+\frac{\mu^{2}}{2}L_{1}^{2}(f)(n+6)^{3}.$ (18)

From (5) and (18) one concludes

$\displaystyle\begin{split}2l(f(x_{k})-&f^{\ast})\leq||\nabla
f(x_{k})||^{2}\\\ &\leq 2||\nabla
f_{\mu}(x_{k})||^{2}+\frac{\mu^{2}}{2}L_{1}^{2}(f)(n+6)^{3}.\end{split}$ (19)

Taking the expectation of this inequality with respect to $\mathcal{U}_{k}$

$E_{\mathcal{U}_{k}}(2l(f(x_{k})-f^{\ast}))\leq\theta_{k}^{2}\leq
2\Xi_{k}^{2}+\frac{\mu^{2}}{2}L_{1}^{2}(f)(n+6)^{3},$ (20)

where $\theta_{k}^{2}=E_{\mathcal{U}_{k}}(||\nabla f(x_{k})||^{2})$. Summing
the inequality from $k=0$ to $k=N$ and dividing it by $N+1$, yields

$\displaystyle\begin{split}\frac{1}{N+1}\sum_{k=0}^{N}&(\Phi_{k}-f^{\ast})\leq\frac{1}{2l(N+1)}\sum_{k=0}^{N}\theta_{k}^{2}\\\
&\leq\frac{1}{l(N+1)}\sum_{k=0}^{N}\Xi_{k}^{2}+\frac{\mu^{2}}{4l}L_{1}^{2}(f)(n+6)^{3}.\end{split}$
(21)

∎

In a practical implementation, we might be interested in identifying the
“best” solution guess at any step $N$. To this aim, define
$\hat{x}_{N}\overset{\mathrm{def}}{=}\arg\underset{x}{\min}[f(x):x\in\\{x_{0},\dots,x_{N}\\}]$.
Hence,

$E_{\mathcal{U}_{N-1}}(f(\hat{x}_{N})-f^{\ast})\leq\frac{1}{N+1}\sum_{k=0}^{N}(\Phi_{k}-f^{\ast}).$
(22)

######  Remark 2 ($\mathrm{RS}_{\mu}$ Complexity, Parameter Selection, and
Solution Error Bound).

If $\mu$ is in the order of
$\mathcal{O}(\frac{\sqrt{l\epsilon}}{n^{3/2}L_{1}(f)})$ and $N$ is in the
order of $\mathcal{O}(\frac{nL_{1}(f)}{l\epsilon})$, it is guaranteed that
$E_{\mathcal{U}_{N-1}}(f(\hat{x}_{N}))-f^{\ast}\leq\epsilon$ for some positive
scalar $\epsilon$.

In comparison with the non-convex smooth case in [5], from the properties of
PL functions the convergence is to a global minimum (see Lemma 1). Also,
comparing the bounds in these two cases, for the case where $l>1$ one can
obtain the same error bound, i.e. $\epsilon$, with fewer iterations. The
number of iterations is inversely proportional with $l$.

###  IV-B Constrained Problem

Now we will focus on the problem (9). This problem can be reformulated as

$\underset{x\in R^{n}}{\min}\;f(x)+\mathcal{I}_{\mathcal{X}}(x),$

where $\mathcal{I}_{\mathcal{X}}(x)$ is the indicator function of
$\mathcal{X}$ . The new scheme for constrained problem is defined in Algorithm
2.

Algorithm 2 $\mathrm{RSc}_{\mu}$

1:Input: $x_{0},h_{k},\mu,N$

2:for $k=0,\dots,N$ do

3: Generate $u$

4: Calculate $g_{\mu}(x_{k})$ using (2)

5: $\bar{x}_{k+1}=x_{k}-h_{k}g_{\mu}(x_{k})$

6: $x_{k+1}=\mathrm{Proj}_{\mathcal{X}}(\bar{x}_{k+1})$

7:end for

8:return $x_{N}$

In Algorithm 2, the projection is used to ensure that the output sequence is
completely in the feasible set.

Before proceeding further, we define the following auxiliary variables:

$P_{\mathcal{X}}(x,g(x),h)\overset{\text{def}}{=}\frac{1}{h}[x-\mathrm{Proj}_{\mathcal{X}}(x-hg(x))],$
(23)
$s_{k}\overset{\text{def}}{=}P_{\mathcal{X}}(x_{k},g_{\mu}(x_{k}),h_{k}),$
(24) $v_{k}\overset{\text{def}}{=}P_{\mathcal{X}}(x_{k},\nabla
f_{\mu}(x_{k}),h_{k}).$ (25)

In the unconstrained case we had $x_{k+1}=x_{k}-h_{k}g_{\mu}(x_{k})$. In the
constrained case, from Algorithm 2, (23), and (24), we can see that
$x_{k+1}=x_{k}-h_{k}s_{k}.$

Before stating the main result, in what follows, we propose a lower bound for
the value of an operator that plays a crucial role in proving the main result
of this section.

######  Lemma 2.

Consider problem (9) where $f(x)\in C^{1,1}$ is a proximal PL function in the
sense of Definition 3, we have

$\displaystyle\begin{split}T(x_{k},L_{1}&(f))\geq 2l(f(x_{k})-f(x^{\ast}))\\\
&-\mu L_{1}(f)^{2}(n+3)^{3/2}d_{x}-2L_{1}(f)d_{x}||\xi_{k}||,\end{split}$ (26)

where

$\displaystyle
T(x,a)\overset{\text{def}}{=}-2a\underset{z\in\mathcal{X}}{\min}$
$\displaystyle\\{\frac{a}{2}||z-x||^{2}$ (27) $\displaystyle+\langle
g_{\mu}(x),z-x\rangle+\mathcal{I}_{\mathcal{X}}(z)-\mathcal{I}_{\mathcal{X}}(x)\\},$

for a positive scalar $a$.

###### Proof.

From the definition of PPL functions, we have

$\displaystyle\begin{split}&2l(f(x_{k})-f(x^{\ast}))\leq Q(x_{k},L_{1}(f))\\\
&\leq-2L_{1}(f)\underset{z\in\mathcal{X}}{\min}(\frac{L_{1}(f)}{2}||z-x_{k}||^{2}\\\
&+\langle\nabla f(x_{k}),z-x_{k}\rangle)\\\
&\leq-2L_{1}(f)\underset{z\in\mathcal{X}}{\min}(\frac{L_{1}(f)}{2}||z-x_{k}||^{2}+\langle\nabla
f_{\mu}(x_{k})\\\ &,z-x_{k}\rangle+\langle\nabla f(x_{k})-\nabla
f_{\mu}(x_{k}),z-x_{k}\rangle)\\\
&\leq-2L_{1}(f)\underset{z\in\mathcal{X}}{\min}(\frac{L_{1}(f)}{2}||z-x_{k}||^{2}+\langle\nabla
f_{\mu}(x_{k})\\\ &,z-x_{k}\rangle-||\nabla f(x_{k})-\nabla
f_{\mu}(x_{k})||||z-x_{k}||),\end{split}$

where the second inequality is a consequence of evaluating (7) for
$x_{k}\in\mathcal{X}$ and noting that $\mathcal{I}_{\mathcal{X}}(x_{k})=0$.
From [5, Lemma 3], for a function $f\in C^{1,1}$ we have

$||\nabla f(x)-\nabla f_{\mu}(x)||\leq\frac{\mu}{2}L_{1}(f)(n+3)^{3/2}.$ (28)

From (28) and $||z-x_{k}||\leq d_{x}$ for all $z\in\mathcal{X}$, we have

$\displaystyle\begin{split}&2l(f(x_{k})-f(x^{\ast}))\leq Q(x_{k},L_{1}(f))\\\
&\leq-2L_{1}(f)\underset{z\in\mathcal{X}}{\min}(\frac{L_{1}(f)}{2}||z-x_{k}||^{2}\\\
&+\langle\nabla f_{\mu}(x_{k}),z-x_{k}\rangle-\frac{\mu
L_{1}(f)(n+3)^{3/2}d_{x}}{2})\\\
&\leq-2L_{1}(f)\underset{z\in\mathcal{X}}{\min}(\frac{L_{1}(f)}{2}||z-x_{k}||^{2}+\langle
g_{\mu}(x_{k})\\\ &,z-x_{k}\rangle-||\xi_{k}||d_{x}-\frac{\mu
L_{1}(f)(n+3)^{3/2}d_{x}}{2}).\end{split}$

where $\xi_{k}\overset{\mathrm{def}}{=}g_{\mu}(x_{k})-\nabla f_{\mu}(x_{k})$.
Rearranging above terms completes the proof. ∎

Before stating the main result regarding the performance of Algorithm 2, we
present the following lemma on the variance of the random oracle $g_{\mu}(x)$.

######  Lemma 3.

Random oracle $g_{\mu}(x)$ is a variance bounded unbiased estimator of $\nabla
f_{\mu}(x)$. That is

$E_{u}[||g_{\mu}(x_{k})-\nabla
f_{\mu}(x_{k})||^{2}]\leq\sigma_{k}^{2},\;\;\sigma_{k}\geq 0.$

###### Proof.

We know that $E_{u}[g_{\mu}(x_{k})]=\nabla f_{\mu}(x_{k})$, so we have
$E_{u}[||g_{\mu}(x_{k})-\nabla f_{\mu}(x_{k})||^{2}]\leq
E_{u}[||g_{\mu}(x_{k})||^{2}]\leq\sigma_{k}^{2}.$ ∎

######  Remark 3.

An upper bound for $E_{u}[||g_{\mu}(x_{k})||^{2}]$ can be obtained. For
example, from [5, Theorem 4] we know for a function $f\in C^{0,0}$ we have
$E_{u}[||g_{\mu}(x)||^{2}]\leq L_{0}(f)^{2}(n+4)^{2}$ and for a function $f\in
C^{1,1}$ we have
$E_{u}[||g_{\mu}(x)||^{2}]\leq\frac{\mu^{2}}{2}L_{1}^{2}(f)(n+6)^{3}+2(n+4)||\nabla
f(x)||^{2}.$ These upper bounds can be used as candidates for
$\sigma_{k}^{2}$.

######  Theorem 2.

Consider problem (9). Let the sequence $\\{x_{k}\\}_{k\geq 0}$ be generated by
RScμ, when $f(x)\in C^{1,1}$ satisfies the proximal PL condition in the sense
of Definition 3. Then, for any $N\geq 0$, with $h_{k}=\frac{1}{L_{1}(f)}$, we
have

$\displaystyle\frac{1}{N+1}\sum_{k=0}^{N}$
$\displaystyle\Phi_{k}-f(x^{\ast})\leq\frac{L_{1}(f)}{l}\frac{f(x_{0})-f(x^{\ast})}{N+1}$
$\displaystyle+\frac{\mu
d_{x}L_{1}(f)^{2}(n+3)^{3/2}}{2l}+\frac{L_{1}(f)d_{x}}{l(N+1)}\sum_{k=0}^{N}\sigma_{k}$
$\displaystyle+\frac{1}{l(N+1)}\sum_{k=0}^{N}\sigma_{k}^{2},$ (29)

where $\Phi_{k}\overset{\mathrm{def}}{=}E_{\mathcal{U}_{k}}[f(x_{k})]$ for all
$k\geq 1$, $\Phi_{0}=f(x_{0})$, $\sigma_{k}$ is given in Lemma 3,
$\mathcal{U}_{k}=\\{u_{0},u_{1},\dots,u_{k-1}\\}$, $n$ is the dimension of
$x$, and $d_{x}$ is the diameter of the feasible set.

###### Proof.

From [5, Equation (12)] we know $f_{\mu}(x)\in C^{1,1}$. Writing (4) for
$f_{\mu}(x)$ at points $x_{k}$ and $x_{k+1}$, yields

$\displaystyle\begin{split}&f_{\mu}(x_{k+1})\leq
f_{\mu}(x_{k})+\frac{L_{1}(f_{\mu})}{2}||x_{k+1}-x_{k}||^{2}\\\
&+\langle\nabla f_{\mu}(x_{k}),x_{k+1}-x_{k}\rangle\\\ &\leq
f_{\mu}(x_{k})+\frac{L_{1}(f_{\mu})}{2}||x_{k+1}-x_{k}||^{2}\\\ &+\langle
g_{\mu}(x_{k}),x_{k+1}-x_{k}\rangle-\langle\xi_{k},x_{k+1}-x_{k}\rangle\\\
&\leq
f_{\mu}(x_{k})-\frac{1}{2L_{1}(f)}T(x_{k},L_{1}(f))+h_{k}\langle\xi_{k},s_{k}\rangle\\\
&\leq f_{\mu}(x_{k})-\frac{1}{2L_{1}(f)}T(x_{k},L_{1}(f))\\\
&+h_{k}\langle\xi_{k},s_{k}-v_{k}\rangle+h_{k}\langle\xi_{k},v_{k}\rangle,\end{split}$

where the third inequality follows from (27), $x_{k}\in\mathcal{X}$, and the
fact that

$\displaystyle T(x_{k},L_{1}(f))$
$\displaystyle=-2L_{1}(f)\underset{z\in\mathcal{X}}{\min}\big{\\{}\frac{L_{1}(f)}{2}||z-x_{k}||^{2}$
$\displaystyle\qquad+\langle g_{\mu}(x_{k}),z-x_{k}\rangle\big{\\}}$
$\displaystyle=-2L_{1}(f)\big{(}\frac{L_{1}(f)}{2}||x_{k+1}-x_{k}||^{2}$
$\displaystyle\qquad+\langle g_{\mu}(x_{k}),x_{k+1}-x_{k}\rangle\big{)}.$

The last equality above follows from the definition of the projection
operator. Taking expected value with respect to $u_{k}$ and considering Lemmas
3 and 4 in the appendix, and $h_{k}=\frac{1}{L_{1}(f)}$, leads to

$\displaystyle\begin{split}E_{u_{k}}[f_{\mu}&(x_{k+1})]\leq f_{\mu}(x_{k})\\\
&-\frac{1}{2L_{1}(f)}E_{u_{k}}[T(x_{k},L_{1}(f))]+\frac{\sigma_{k}^{2}}{L_{1}(f)}.\end{split}$

Taking the expectation with respect to $\mathcal{U}_{k-1}$, we have

$E_{\mathcal{U}_{k}}[T(x_{k},L_{1}(f))]\leq
2L_{1}(f)(\Phi_{k+1}-\Phi_{k})+2\sigma_{k}^{2}.$

Summing over $k=0$ to $k=N$ and dividing it by $N+1$, results in

$\displaystyle\frac{1}{N+1}\sum_{k=0}^{N}E_{\mathcal{U}_{k}}$
$\displaystyle[T(x_{k},L_{1}(f))]\leq\frac{2}{N+1}\sum_{k=0}^{N}\sigma_{k}^{2}+$
$\displaystyle 2L_{1}(f)\frac{f(x_{0})-f(x^{\ast})}{N+1}.$

Also, taking the expectation of (26) with respect to $u_{k}$ and then
$\mathcal{U}_{k-1}$, summing over $k=0$ to $k=N$, dividing it by $N+1$ and
using Lemma 5 from the appendix, yield

$\displaystyle\begin{split}&\frac{2l}{N+1}\sum_{k=0}^{N}\Phi_{k}-f(x^{\ast})\leq\mu
d_{x}L_{1}(f)^{2}(n+3)^{3/2}\\\
&+\frac{1}{N+1}\sum_{k=0}^{N}E_{\mathcal{U}_{k}}[T(x_{k},L_{1}(f))]+\frac{2L_{1}(f)d_{x}}{N+1}\sum_{k=0}^{N}\sigma_{k}.\end{split}$

Thus, we have

$\displaystyle\begin{split}\frac{1}{N+1}\sum_{k=0}^{N}&\Phi_{k}-f(x^{\ast})\leq\frac{L_{1}(f)}{l}\frac{f(x_{0})-f(x^{\ast})}{N+1}\\\
&+\frac{\mu
d_{x}L_{1}(f)^{2}(n+3)^{3/2}}{2l}+\frac{L_{1}(f)d_{x}}{l(N+1)}\sum_{k=0}^{N}\sigma_{k}\\\
&+\frac{1}{l(N+1)}\sum_{k=0}^{N}\sigma_{k}^{2}.\end{split}$

∎

In a practical implementation of the algorithm we might be interested in
keeping track of the best guess for the optimum solution at any given step
$N$. As in the unconstrained case, let this best guess be denoted by
$\hat{x}_{N}=\arg\underset{x}{\min}[f(x):x\in\\{x_{0},\dots,x_{N}\\}]$. Thus,

$E_{\mathcal{U}_{N-1}}(f(\hat{x}_{N})-f^{\ast})\leq\frac{1}{N+1}\sum_{k=0}^{N}(\Phi_{k}-f^{\ast}).$

######  Remark 4 ($\mathrm{RSc}_{\mu}$ Complexity, Parameter Selection, and
Solution Error Bound).

If $\mu\leq\frac{l\epsilon}{d_{x}L_{1}(f)^{2}(n+3)^{3/2}}$, $N$ is in the
order of $\mathcal{O}\left(\frac{L_{1}(f)}{l\epsilon}\right)$, and denoting
$\sigma=\underset{k}{\max}[\sigma_{k}:k\in\\{0,\dots,N\\}]$, then

$\displaystyle
E_{\mathcal{U}_{N-1}}(f(\hat{x}_{N})-f^{\ast})\leq\epsilon+\frac{L_{1}(f)d_{x}\sigma}{l}+\frac{\sigma^{2}}{l},$
(30)

for some positive scalar $\epsilon$. Thus, we can guarantee that there exists
an integer $\bar{N}$ such that for all $N\geq\bar{N}$,
$E_{\mathcal{U}_{N-1}}f(x_{k})$ is in a neighbourhood of $f^{\ast}$ for an
appropriate choice of $\mu$.

Similar phenomena are observed consistently in the literature on constrained
non-convex problems, where this effect has been reported for different
algorithms, see e.g. [4], [10]. Note that similarly to the unconstrained case
$l$ appears in the denominator of the iteration number order term.
Additionally, the two terms on the right-hand side of (30) are inversely
proportional with $l$. Also, it can be seen in this case due to the choice of
$h_{k}$, the number of iterations is not dependent on the dimension of $x$,
but we know that $h_{k}\in(0,\frac{1}{L_{1}(f)}]$ and in fact $N$ is in the
order of $\mathcal{O}(\frac{1}{lh_{k}\epsilon})$ for the more general case.

##  V Numerical Examples

In this section, we consider two scenarios. In both scenarios we study the
performance of Algorithms 1 and 2 applied to unconstrained and constrained
least square problems. Specifically, the objective function is assumed to be
$||Ax-b||^{2}$, where $A\in R^{m\times n}$ ($n\geq m$), $x\in R^{n}$, and
$b\in R^{m}$. The function is not strongly convex but it satisfies the PL
condition with $l=2||A^{T}||^{2}$ and is in $C^{1,1}$ with
$L_{1}(f)=2||A^{T}A||$.

#### Scenario 1

In the first scenario, we consider an unconstrained least-squares problem of
the form introduced above with $m=100$ and $n=1000$. In light of Remark 2, we
choose $\epsilon=0.01$ and consequently we set $\mu=10^{-7}$ and $N=200000$.
Matrix $A$ rows are sampled from $\mathcal{N}(0,I_{m})$ and
$b=A\bar{x}+\omega$, where $\bar{x}$ sampled from $\mathcal{N}(0,1)$ and
$\omega$ from $\mathcal{N}(0,0.01)$. Moreover, the initial condition vector is
sampled from $\mathcal{N}(0,1)$. We explore the performance of the algorithm
for two step sizes $\frac{1}{4(n+4)L_{1}(f)}\approx 10^{-7}$ and $10^{-6}$. We
repeat the example for 25 times. The empirical mean of best guess for the
optimum solution ($f(\hat{x}_{N})$) over 25 runs and upper bound calculated in
(22) is presented in Fig. 1. As it was shown in Theorem 1, by increasing the
step size the convergence to a neighbourhood of the solution would be faster,
but this comes at the expense of larger error bounds.

#### Scenario 2

In this scenario we consider a constrained least squares problem with the same
problem parameter choices as the previous scenario except for $\mu$. From
Remark 4, given the same value for $\epsilon$ as the previous case, we choose
$\mu=10^{-10}$ and $N=200000$. The constraint set $\mathcal{X}$ is assumed to
be $\mathcal{X}=\\{x\in\mathbb{R}^{n}|x_{i}\in[-0.5,0.5],\forall
i\in\\{1,\dots,n\\}\\}$. We explore the performance of the algorithm for two
step sizes $\frac{1}{nL_{1}(f)}\approx 10^{-7}$ and $10^{-6}$ (both are less
than $\frac{1}{L_{1}(f)}$). The empirical mean of objective function value in
iterations generated by Algorithm 2 over $25$ runs are depicted in Fig. 2. The
effect of additional error terms in the constrained case can be seen by
comparing the figures at point $10^{4}$ iterations for $h_{k}=10^{-6}$ cases
and at point $10^{5}$ iterations for $h_{k}=10^{-7}$cases.

Figure 1: The evolution of the empirical mean of $f(\hat{x}_{k})$ and the
calculated upper bound versus the number of iterations in Scenario 1. Figure
2: The evolution of $f(x_{k})$ versus the number of iterations. Note that in
this case $f(x^{\ast})=0$ in Scenario 2.

##  VI Conclusion and Future Research Directions

The application of a zeroth-order scheme using random oracles for minimising
PL functions with or without constraints was discussed. For the unconstrained
problem, the convergence properties of the proposed algorithm were studied.
Additionally, the complexity bounds for the algorithm along with guidelines
for selecting algorithm parameters were introduced. Next, a generalisation of
the PL inequality was exploited to establish the convergence of the algorithm
for solving constrained problems. Similar to the unconstrained case,
complexity bounds were derived. Numerical examples were presented to
illustrate the theoretical results. An immediate future step is extending the
analysis of the constrained case to the case where the constraint set is
unbounded. Another possible future direction is applying similar techniques
for solving minimax problems using zeroth-order oracles of the type studied in
this paper.

###  -A Auxiliary Lemmas

######  Lemma 4.

Let $\xi_{k}\overset{\mathrm{def}}{=}g_{\mu}(x_{k})-\nabla f_{\mu}(x_{k}).$
From [4, Proposition 1], it implies that
$\langle\xi_{k},s_{k}-v_{k}\rangle\leq||\xi_{k}||^{2}$.

######  Lemma 5.

For a function $f(x)$ with smoothed version $f_{\mu}(x)$ and random oracle
$g_{\mu}(x)$, we have

$E_{u_{k}}[||\xi_{k}||]\leq\sigma_{k},\;\;\sigma_{k}\geq 0,$

where $\xi_{k}=g_{\mu}(x_{k})-\nabla f_{\mu}(x_{k})$.

###### Proof.

Due to Lemma 3, $E_{u_{k}}[||\xi_{k}||^{2}]\leq\sigma_{k}^{2}.$ From Jensen’s
inequality, we have

$E_{u_{k}}[||\xi_{k}||]^{2}\leq E_{u_{k}}[||\xi_{k}||^{2}].$

Thus, $E_{u_{k}}[||\xi_{k}||]\leq\sigma_{k}.$ ∎

## References

  * [1] Alejandro I Maass, Chris Manzie, Iman Shames and Hayato Nakada  “Zeroth-Order Optimization on Subsets of Symmetric Matrices With Application to MPC Tuning”  In _IEEE Transactions on Control Systems Technology_ 30.4 IEEE, 2021, pp. 1654–1667
  * [2] Alonso Marco et al.  “Automatic LQR tuning based on Gaussian process global optimization”  In _2016 IEEE international conference on robotics and automation (ICRA)_ , 2016, pp. 270–277  IEEE
  * [3] Léon Bottou, Frank E Curtis and Jorge Nocedal  “Optimization methods for large-scale machine learning”  In _SIAM review_ 60.2 SIAM, 2018, pp. 223–311
  * [4] Saeed Ghadimi, Guanghui Lan and Hongchao Zhang  “Mini-batch stochastic approximation methods for nonconvex stochastic composite optimization”  In _Mathematical Programming_ 155.1-2 Springer, 2016, pp. 267–305
  * [5] Yurii Nesterov and Vladimir Spokoiny  “Random Gradient-Free Minimization of Convex Functions”  In _Foundations of Computational Mathematics_ 17.2, 2017, pp. 527–566  DOI: 10.1007/s10208-015-9296-2
  * [6] Boris T Polyak  “Gradient methods for solving equations and inequalities”  In _USSR Computational Mathematics and Mathematical Physics_ 4.6 Elsevier, 1964, pp. 17–32
  * [7] Chaoyue Liu, Libin Zhu and Mikhail Belkin  “Loss landscapes and optimization in over-parameterized non-linear systems and neural networks”  In _Applied and Computational Harmonic Analysis_ 59 Elsevier, 2022, pp. 85–116
  * [8] Hamed Karimi, Julie Nutini and Mark Schmidt  “Linear convergence of gradient and proximal-gradient methods under the polyak-łojasiewicz condition”  In _Machine Learning and Knowledge Discovery in Databases: European Conference, ECML PKDD 2016, Riva del Garda, Italy, September 19-23, 2016, Proceedings, Part I 16_ , 2016, pp. 795–811  Springer
  * [9] Sotirios-Konstantinos Anagnostidis, Aurelien Lucchi and Youssef Diouane  “Direct-search for a class of stochastic min-max problems”  In _International Conference on Artificial Intelligence and Statistics_ , 2021, pp. 3772–3780  PMLR
  * [10] Sijia Liu et al.  “Zeroth-order stochastic variance reduction for nonconvex optimization”  In _Advances in Neural Information Processing Systems_ 31, 2018


# Magnetic point groups and space groups

Ron Lifshitz
School of Physics and Astronomy
Raymond and Beverly Sackler Faculty of Exact Sciences
Tel Aviv University, Tel Aviv 69978, Israel
E-mail<EMAIL_ADDRESS>
(10 May 2004)

##  1 Introduction

Magnetic groups—also known as antisymmetry groups, Shubnikov groups, Heesch
groups, Opechowski-Guccione groups, as well as dichromatic, 2-color, or simply
black-and-white groups—are the simplest extension to standard point group and
space group theory. They allow directly to describe, classify, and study the
consequences of the symmetry of crystals, characterized by having a certain
property, associated with each crystal site, that can take one of two possible
values. This is done by introducing a single external operation of order two
that interchanges the two possible values everywhere throughout the crystal.
This operation can be applied to the crystal along with any of the standard
point group or space group operations, and is denoted by adding a prime to the
original operation. Thus, any rotation $g$ followed by this external operation
is denoted by $g^{\prime}$.

We begin here with a few typical examples of this two-valued property, some of
which are illustrated in Figure 1. In section 2 we discuss the notion of a
magnetic point group, followed by a discussion on magnetic space groups in
section 3. In section 4 we describe one of the most direct consequences of
having magnetic symmetry in crystals which is the extinction of magnetic Bragg
peaks in neutron diffraction patterns. We finish in section 5 by mentioning
the generalization of magnetic point groups to cases where the property
associated with each crystal site is not limited to having only one of two
possible values.

Let us begin by considering the structure of crystals of the cesium chloride
type. In these crystals the atoms are located on the sites of a body-centered
cubic (bcc) lattice. Atoms of one type are located on the simple cubic lattice
formed by the corners of the cubic cells, and atoms of the other type are
located on the simple cubic lattice formed by the body centers of the cubic
cells. The parent bcc lattice is partitioned in this way into two identical
sublattices, related by a body-diagonal translation. One may describe the
symmetry of the cesium chloride structure as having a simple cubic lattice of
ordinary translations, with a basis containing one cesium and one chlorine
atom. Alternatively, one may describe the symmetry of the crystal as having
twice as many translations, forming a bcc lattice, half of which are primed to
indicate that they are followed by the external operation that exchanges the
chemical identities of the two types of atoms. A similar situation occurs in
crystals whose structure is of the sodium chloride type in which a simple
cubic lattice is partitioned into two identical face-centered cubic
sublattices.

Another typical example is the orientational ordering of magnetic moments
(spins), electric dipole moments, or any other tensorial property associated
with each site in a crystal. Adopting the language of spins, if these can take
only one of two possible orientations—“up” and “down”—as in simple
antiferromagnets, we have the same situation as above, with the two spin
orientations replacing the two chemical species of atoms. In the case of spins
the external prime operation is a so-called antisymmetry operation that
reverses the signs of the spins. Physically, one may think of using the
operation of time inversion to reverse the directions of all the spins without
affecting anything else in the crystal.

Finally, consider a periodic or quasiperiodic scalar function of space $f({\bf
r})$, as in Figure 1(e), whose average value is zero. It might be possible to
extend the point group or the space group describing the symmetry of the
function $f({\bf r})$ by introducing an external prime operation that changes
the sign of $f({\bf r})$.

For the sake of clarity we shall adopt the picture of up-and-down spins, and
the use of time inversion to exchange them, for the remaining of our
discussion. It should be emphasized, though, that most of what we say here
(except for the discussion of extinctions in section 4) applies equally to all
the other two-valued properties. We describe the magnetic crystal using a
scalar spin density field $S({\bf r})$ whose magnitude gives the magnitude of
an atomic magnetic moment, or some coarse-grained value of the magnetic
moment, at position ${\bf r}$. The sign of $S({\bf r})$ gives the direction of
the spin—a positive sign for up spins and a negative sign for down spins.
Thus, the function $S({\bf r})$ can be a discrete set of delta functions
defined on the crystal sites as in Figure 1(b), or a continuous spin density
field as in Figure 1(e).

Figure 1: Several realizations of the simple one-dimensional magnetic space
group $p_{p}\bar{1}$. Note that the number of symmetry translations is doubled
if one introduces an operation $e^{\prime}$ that interchanges the two possible
values of the property, associated with each crystal site. Also note that
spatial inversion $\bar{1}$ can be performed without a prime on each crystal
site, and with a prime ($\bar{1}^{\prime}$) between every two sites. (a) An
abstract representation of the two possible values as the two colors—black and
white. The operation $e^{\prime}$ is the nontrivial permutation of the two
colors. (b) A simple antiferromagnetic arrangement of spins. The operation
$e^{\prime}$ is time inversion which reverses the signs of the spins. (c) A
ferromagnetic arrangement of two types of spins, where $e^{\prime}$ exchanges
them as in the case of two colors. (d) A 1-dimensional version of salt.
$e^{\prime}$ exchanges the chemical identities of the two atomic species. (e)
A function $f(x)$ whose overall sign changes by application of the operation
$e^{\prime}$.

##  2 Magnetic point groups

A $d$-dimensional magnetic point group $G_{M}$ is a subgroup of
$O(d)\\!\times\\!1^{\prime}$, where $O(d)$ is the group of $d$-dimensional
rotations, and $1^{\prime}$ is the time inversion group containing the
identity $e$ and the time inversion operation $e^{\prime}$. Note that we use
the term “rotation” to refer to proper as well as improper rotations such as
mirrors.

Three cases exist: (1) All rotations in $G_{M}$ appear with and without time
inversion; (2) Half of the rotations in $G_{M}$ are followed by time inversion
and the other half are not; and (3) No rotation in $G_{M}$ is followed by time
inversion. More formally, if $G$ is any subgroup of $O(d)$, it can be used to
form at most three types of magnetic point groups, as follows:

  1. 1.

$G_{M}=G\\!\times\\!1^{\prime}$. Thus, each rotation in $G$ appears in $G_{M}$
once by itself and once followed by time inversion. Note that in this case
$e^{\prime}\in G_{M}$.

  2. 2.

$G_{M}=H+g^{\prime}H$, where $G=H+gH$, and $g\not\in H$. Thus, exactly half of
the rotation in $G$, which belong to its subgroup $H$, appear in $G_{M}$ by
themselves, and the other half, belonging to the coset $gH$, are followed by
time inversion. Note that in this case $e^{\prime}\not\in G_{M}$.

  3. 3.

$G_{M}=G$. Thus, $G_{M}$ contains rotations none of which are followed by time
inversion.

Enumeration of magnetic point groups is thus straightforward. Any ordinary
point group $G$ [i.e. any subgroup, usually finite, of $O(d)$] is trivially a
magnetic point group of type 3, and along with the time inversion group gives
another magnetic point group of type 1, denoted as $G1^{\prime}$. One then
lists all distinct subgroups $H$ of index 2 in $G$, if there are any, to
construct additional magnetic point groups of type 2. These are denoted either
by the group-subgroup pair $G(H)$, or by using the International (Hermann-
Mauguin) symbol for $G$ and adding a prime to those elements in the symbol
that do not belong to the subgroup $H$ (and are therefore followed by time
inversion). The set of magnetic groups—here magnetic point groups
$G_{M}$—formed in this way from a single ordinary group $G$—here an ordinary
point group—is sometimes called the magnetic superfamily of the group $G$. For
example, the orthorhombic point group $G=mm2=\\{e,m_{x},m_{y},2_{z}\\}$ has
three subgroups of index 2, $H_{1}=\\{e,m_{x}\\}$, $H_{2}=\\{e,m_{y}\\}$, and
$H_{3}=\\{e,2_{z}\\}$, of which the first two are equivalent through a
relabeling of the $x$ and $y$ axes. This yields a total of four magnetic point
groups: $mm2$, $mm21^{\prime}$, $m^{\prime}m2^{\prime}$ (equivalently
expressed as $mm^{\prime}2^{\prime}$), and $m^{\prime}m^{\prime}2$.

Magnetic point groups of type 3 are equivalent to ordinary non-magnetic point
groups and are listed here only for the sake of completeness. Strictly
speaking, they can be said to describe the symmetry of ferromagnetic
structures with no spin-orbit coupling. We shall not consider groups of this
type any further. Magnetic point groups of type 2 can be used to describe the
point symmetry of finite objects, as described in 2.1, as well as that of
infinite crystals, as described in 2.2. Magnetic point groups of type 1 can be
used to describe the point symmetry of infinite crystals, but because they
include $e^{\prime}$ as a symmetry element they cannot be used to describe the
symmetry of finite objects, as we now define.

###  2.1 Magnetic point groups of finite objects

The magnetic symmetry of a finite object, such as a molecule or a cluster
containing an equal number of two types of spins, can be described by a
magnetic point group. We say that a primed or unprimed rotation from
$O(d)\\!\times\\!1^{\prime}$ is in the magnetic point group $G_{M}$ of a
finite object in $d$ dimensions, if it leaves the object invariant. Clearly,
only magnetic point groups of type 2, listed above, can describe the symmetry
of finite magnetically-ordered structures. This is because time inversion
$e^{\prime}$ changes the orientations of all the spins and therefore cannot
leave the object invariant unless it is accompanied by a nontrivial rotation.
It should be mentioned that in this context, point groups of type 1 are
sometimes called “gray” groups, as they describe the symmetry of “gray”
objects which stay invariant under the exchange of black and white.

###  2.2 Magnetic point groups of crystals

The magnetic point group of a $d$-dimensional magnetically-ordered periodic
crystal, containing an equal number of up and down spins, is defined as the
set of primed or unprimed rotations from $O(d)\\!\times\\!1^{\prime}$ that
leave the crystal invariant to within a translation. The magnetic point group
of a crystal can be either of the first or of the second type listed above. It
is of type 1 if time inversion by itself leaves the crystal invariant to
within a translation. Recall that in this case any rotation in the magnetic
point group can be performed either with or without time inversion. If time
inversion cannot be combined with a translation to leave the crystal
invariant, the magnetic point group is of type 2, in which case half the
rotations are performed without time inversion and the other half with time
inversion.

(a) (b)

Figure 2: Periodic antiferromagnets. (a) A square crystal with magnetic point
group $G_{M}=4mm1^{\prime}$ and magnetic space group $p_{p}4mm$ (also denoted
$p_{c}4mm$). (b) A hexagonal crystal with magnetic point group
$G_{M}=6^{\prime}mm^{\prime}$ and magnetic space group
$p6^{\prime}mm^{\prime}$.

Figure 2(a) shows an example of a magnetically-ordered crystal whose magnetic
point group is of type 1. This is a square crystal with magnetic point group
$G_{M}=4mm1^{\prime}$ where time inversion can be followed by a translation to
leave the crystal invariant. Figure 2(b) shows an example of a crystal whose
magnetic point group is of type 2. This is a hexagonal crystal with point
group $G_{M}=6^{\prime}mm^{\prime}$. Note that all right-side-up triangles
contain a dark circle (spin up) and all up-side-down triangles contain a light
circle (spin down). Time inversion, exchanging the two types of spins, cannot
be combined with a translation to recover the original crystal. Time inversion
must be combined with an operation, such as the 6-fold rotation or the
horizontal mirror, that interchanges the two types of triangles to recover the
original crystal. Note that the vertical mirror (the first $m$ in the
International symbol) leaves the crystal invariant without requiring time
inversion and is therefore unprimed in the symbol.

(a) (b)

Figure 3: Quasiperiodic antiferromagnets. (a) An octagonal quasicrystal with
magnetic point group $G_{M}=8mm1^{\prime}$ and magnetic space group
$p_{p}8mm$, first appeared in an article by Niizeki [J. Phys. A: Math. Gen. 23
(1990) 5011]. (b) A decagonal quasicrystal with magnetic point group
$G_{M}=10^{\prime}m^{\prime}m$ and magnetic space group
$p10^{\prime}m^{\prime}m$, based on the tiling of Li, Dubois, and Kuo [Philos.
Mag. Lett. 69 (1994) 93].

More generally, the magnetic point group of a $d$-dimensional magnetically-
ordered quasiperiodic crystal (quasicrystal), is defined as the set of primed
or unprimed rotations from $O(d)\\!\times\\!1^{\prime}$ that leave the crystal
indistinguishable. This means that the rotated and unrotated crystals contain
the same spatial distributions of finite clusters of spins of arbitrary size.
The two are statistically the same though not necessarily identical. For the
special case of periodic crystals, the requirement of indistinguishability
reduces to the requirement of invariance to within a translation. Figure 3
shows two quasiperiodic examples, analogous to the two periodic examples of
Figure 2. Figure 3(a) shows an octagonal crystal with magnetic point group
$G_{M}=8mm1^{\prime}$. One can see that time inversion rearranges the spin
clusters in the crystal but they all still appear in the crystal with the same
spatial distribution. This is because any finite spin cluster and its time-
reversed image appear in the crystal with the same spatial distribution.
Figure 3(b) shows a decagonal crystal with magnetic point group
$G_{M}=10^{\prime}m^{\prime}m$. In this case, time inversion does not leave
the crystal indistinguishable. It must be combined either with odd powers of
the 10-fold rotation, or with mirrors of the vertical type.

##  3 Magnetic space groups

The full symmetry of a magnetically-ordered crystal, described by a scalar
spin density field $S({\bf r})$, is given by its magnetic space group ${\cal
G}_{M}$. We said that the magnetic point group $G_{M}$ is the set of primed or
unprimed rotations that leave a periodic crystal invariant to within a
translation, or more generally, leave a quasiperiodic crystal
indistinguishable. We still need to specify the distinct sets of translations
${\bf t}_{g}$ or ${\bf t}_{g^{\prime}}$ that can accompany the rotations $g$
or $g^{\prime}$ in a periodic crystal to leave it invariant, or provide the
more general characterization of the nature of the indistinguishability in the
case of quasicrystals.

###  3.1 Periodic crystals

As in the case of point groups, one can take any space group $\cal G$ of a
periodic crystal and form its magnetic superfamily, consisting of one trivial
group, ${\cal G}_{M}=\cal G$, one gray group ${\cal G}_{M}={\cal
G}1^{\prime}$, and possibly nontrivial magnetic groups of the form ${\cal
G}_{M}={\cal H}+(g^{\prime}|{\bf t}_{g^{\prime}}){\cal H}$, where ${\cal H}$
is a subgroup of index 2 of $\cal G$. Only the latter, nontrivial, groups are
relevant as magnetic space groups consisting of operations that leave the
magnetically-ordered crystal invariant. These nontrivial magnetic space groups
are divided into two types depending on their magnetic point groups $G_{M}$:

  1. 1.

Class-equivalent magnetic space groups: The magnetic point group $G_{M}$ is of
type 1 and therefore all rotations in $G$ appear with and without time
inversion. In this case the point groups of $\cal G$ and $\cal H$ are the
same, but the Bravais lattice of $\cal H$ contains exactly half of the
translations that are in the Bravais lattice of $\cal G$.

  2. 2.

Translation-equivalent magnetic space groups: The magnetic point group $G_{M}$
is of type 2. Here the Bravais lattices of $\cal G$ and $\cal H$ are the same,
but the point group $H$ of $\cal H$ is a subgroup of index 2 of the point
group $G$ of $\cal G$.

Enumeration of translation-equivalent magnetic space group types is very
simple. Given a space group $\cal G$, all that one needs to do is to consider
all subgroups $H$ of index 2 of the point group $G$ of $\cal G$, as explained
in section 2. The only additional consideration is that with an underlying
Bravais lattice there may be more than one way to orient the subgroup $H$
relative to the lattice. Take, for example, the magnetic point group
$m^{\prime}m2^{\prime}$ which we earlier said was equivalent to
$mm^{\prime}2^{\prime}$. If the lattice is an $A$-centered orthorhombic
lattice, then the $x$\- and $y$-axes are no longer equivalent and one needs to
consider both options as distinct groups. We say that $Am^{\prime}m2^{\prime}$
and $Amm^{\prime}2^{\prime}$ belong to the same magnetic geometric crystal
class (i.e. have the same magnetic point group), but to two distinct magnetic
arithmetic crystal classes. Translation-equivalent magnetic space groups are
denoted by taking the International symbol for $\cal G$ and adding a prime to
those elements in the symbol that do not belong to $\cal H$ and are therefore
followed by time inversion.

Enumeration of class-equivalent magnetic space group types can proceed in two
alternative ways. Given a space group $\cal G$, one may consider all distinct
sublattices of index 2 of the Bravais lattice of $\cal G$. Alternatively,
given a space group $\cal H$, one may consider all distinct superlattices of
index 2 of the Bravais lattice of $\cal H$. Also note that if one prefers to
look at the reciprocal lattices in Fourier space, instead of the lattices of
real-space translations, the roles of superlattice and sublattice are
exchanged. Although the lattice of $\cal G$ contains twice as many points as
the lattice of $\cal H$, the reciprocal lattice of $\cal G$ contains exactly
half of the points of the reciprocal lattice of $\cal H$.

Because of the different enumeration methods there are conflicting notations
in the literature for class-equivalent magnetic space groups. The Opechowski-
Guccione (OG) notation follows the first approach (sublattices of $\cal G$),
taking the International symbol for the space group $\cal G$ with a subscript
on the Bravais class symbol denoting the Bravais lattice of $\cal H$. The
Belov notation follows the second approach (superlattices of $\cal H$), taking
the International symbol for the space group $\cal H$ with a subscript on the
Bravais class symbol denoting one of the primed translations ${\bf
t}_{e^{\prime}}$ in the coset $(e^{\prime}|{\bf t}_{e^{\prime}}){\cal H}$. The
Lifshitz notation, which follows a third approach (sublattices of $\cal H$ in
Fourier space), resembles that of Belov but generalizes more easily to
quasiperiodic crystals.

There is an additional discrepancy in the different notations which has been
the cause of some confusion over the years. It is due to the fact that when
$e^{\prime}\in G_{M}$ ordinary mirrors or rotations when unprimed may become
glide planes or screw axes when primed, and vice versa. The only consistent
way to avoid this confusion is to use only unprimed elements for the symbols
of class-equivalent magnetic space group types. This is the approach adopted
by both the Belov and the Lifshitz notations, but unfortunately not by the OG
notation, having introduced errors into their list of magnetic space group
types. In any case, there is no need to leave the $1^{\prime}$ at the end of
the symbol, as it is clear from the existence of a subscript on the lattice
symbol that the magnetic point group $G_{M}$ contains $e^{\prime}$.

Let us list all the class-equivalent magnetic space group types with point
group $432$. In the cubic system there are two lattice-sublattice pairs: (1)
The simple cubic lattice $P$ has a face-centered sublattice $F$ of index 2;
and (2) The body-centered cubic lattice $I$ has a simple cubic sublattice $P$
of index 2. Also recall that the reciprocal of $F$ is a bcc lattice in Fourier
space, denoted by $I^{*}$, and the reciprocal of $I$ is an fcc lattice,
denoted $F^{*}$. In the Belov notation we have $F_{s}432$, $F_{s}4_{j}32$,
$P_{I}432$, and $P_{I}4_{j}32$ with $j=1,2,3$. The subscript $s$ is for
“simple”. In the Lifshitz notation, using Fourier space lattices, these groups
become $I^{*}_{P}432$, $I^{*}_{P}4_{j}32$, $P_{F^{*}}432$, and
$P_{F^{*}}4_{j}32$. In the OG notation these groups should become $P_{F}432$,
$P_{F}4_{j}32$, $I_{P}432$, and $I_{P}4_{j}32$. Instead, they list $P_{F}432$,
$P_{F}4_{2}32$, $I_{P}432$, $I_{P}4_{1}32$ $I_{P}4^{\prime}32^{\prime}$, and
$I_{P}4^{\prime}_{1}32^{\prime}$, using primes inconsistently and clearly
missing two groups. We therefore recommend to refrain from using primes when
denoting class-equivalent magnetic space groups.

We refer the reader to references in the “Further reading” section for
complete lists of magnetic space group types. Let us only summarize some
statistics. Starting from the 17 space group types in 2 dimensions, also known
as the 17 plane groups, one can form a total of 80 magnetic space group types
as follows: 17 trivial groups, 17 gray groups, and 46 nontrivial groups of
which 18 are class-equivalent and 28 are translation-equivalent. Starting from
the 230 space group types in 3 dimensions one can form a total of 1651
magnetic space group types as follows: 230 trivial groups, 230 gray groups,
and 1191 nontrivial groups of which 517 are class-equivalent and 674 are
translation-equivalent. These numbers have no particular significance other
than the fact that they are surprisingly large.

###  3.2 Quasiperiodic crystals

####  3.2.1 Lattices and Bravais classes

Magnetically-ordered quasiperiodic crystals, in general, possess no lattices
of real space translations that leave them invariant, but they do have
reciprocal lattices (or Fourier modules) $L$ that can be inferred directly
from their neutron diffraction diagrams. To be a bit more specific, consider
spin density fields with well defined Fourier transforms

$S({\bf r})=\sum_{{\bf k}\in L}S({\bf k})e^{i{\bf k}\cdot{\bf r}},$ (1)

in which the set $L$ contains at most a countable infinity of plane waves. In
fact, the reciprocal lattice $L$ of a magnetic crystal can be defined as the
set of all integral linear combinations of wave vectors ${\bf k}$ determined
from its neutron diffraction diagram, as we expect to see a magnetic Bragg
peak at every ${\bf k}$ in the Fourier transform of $S({\bf r})$, unless it is
forbidden by symmetry, as explained in section 4 below. The rank $D$ of $L$ is
the smallest number of vectors needed to generate $L$ by integral linear
combinations. If $D$ is finite the crystal is quasiperiodic. If, in
particular, $D$ is equal to the number $d$ of spatial dimensions, and $L$
extends in all $d$ directions, the crystal is periodic. In this special case
the magnetic lattice $L$ is reciprocal to a lattice of translations in real
space that leave the magnetic crystal invariant. The reciprocal lattices $L$
of magnetic crystals are classified into Bravais classes, just like those of
ordinary non-magnetic crystals.

####  3.2.2 Phase functions and space group types

The precise mathematical statement of indistinguishability, used earlier to
define the magnetic point group, is the requirement that any symmetry
operation of the magnetic crystal leave invariant all spatially-averaged
$n^{th}$-order autocorrelation functions of its spin-density field,

$C^{(n)}({\bf r}_{1},\ldots,{\bf r}_{n})=\lim_{V\to\infty}{1\over
V}\int_{V}d{\bf r}S({\bf r}_{1}-{\bf r})\cdots S({\bf r}_{n}-{\bf r}).$ (2)

One can prove that two quasiperiodic spin density fields $S({\bf r})$ and
$\hat{S}({\bf r})$ are indistinguishable in this way, if their Fourier
coefficients, defined in Eq. (1), are related by

$\hat{S}({\bf k})=e^{2\pi i\chi({\bf k})}S({\bf k}),$ (3)

where $\chi$ is a real-valued linear function (modulo integers) on $L$, called
a gauge function. By this we mean that for every pair of wave vectors ${\bf
k}_{1}$ and ${\bf k}_{2}$ in the magnetic lattice $L$,

$\chi({\bf k}_{1}+{\bf k}_{2})\equiv\chi({\bf k}_{1})+\chi({\bf k}_{2}),$ (4)

where “$\equiv$” means equal to within adding an integer. Thus, for each
element $g$ or $g^{\prime}$ in the magnetic point group $G_{M}$ of the
crystal, that by definition leaves the crystal indistinguishable, there exists
such a gauge function $\Phi_{g}({\bf k})$ or $\Phi_{g^{\prime}}({\bf k})$,
called a phase function, satisfying

$S(g{\bf k})=\cases{e^{2\pi i\Phi_{g}({\bf k})}S({\bf k}),&$g\in
G_{M}$,\cr\cr-e^{2\pi i\Phi_{g^{\prime}}({\bf k})}S({\bf k}),&$g^{\prime}\in
G_{M}$.}$ (5)

Since for any $g,h\in G$, $S([gh]{\bf k})=S(g[h{\bf k}])$, the corresponding
phase functions for elements in $G_{M}$, whether primed or not, must satisfy
the group compatibility condition,

$\Phi_{g^{*}h^{\dagger}}({\bf k})\equiv\Phi_{g^{*}}(h{\bf
k})+\Phi_{h^{\dagger}}({\bf k}),$ (6)

where the asterisk and the dagger denote optional primes. A spin space group,
describing the symmetry of a magnetic crystal, whether periodic or aperiodic,
is thus given by a lattice $L$, a magnetic point group $G_{M}$, and a set of
phase functions $\Phi_{g}({\bf k})$ and $\Phi_{g^{\prime}}({\bf k})$,
satisfying the group compatibility condition (6).

Magnetic space groups are classified in this formalism into magnetic space
group types by organizing sets of phase functions satisfying the group
compatibility condition (6) into equivalence classes. Two such sets $\Phi$ and
$\hat{\Phi}$ are equivalent if: (I) they describe indistinguishable spin
density fields related as in (3) by a gauge function $\chi$; or (II) they
correspond to alternative descriptions of the same crystal that differ by
their choices of absolute length scales and spatial orientations. In case (I)
$\Phi$ and $\hat{\Phi}$ are related by a gauge transformation,

${\hat{\Phi}}_{g^{*}}({\bf k})\equiv\Phi_{g^{*}}({\bf k})+\chi(g{\bf k}-{\bf
k}),$ (7)

where, again, the asterisk denotes an optional prime.

In the special case that the crystal is periodic ($D=d$) it is possible to
replace each gauge function by a corresponding $d$-dimensional translation
${\bf t}$, satisfying $2\pi\chi({\bf k})={\bf k}\cdot{\bf t}$, thereby
reducing the requirement of indistinguishability to that of invariance to
within a translation. Only then does the point group condition (5) become

$S(g{\bf r})=\cases{S({\bf r}-{\bf t}_{g}),&$g\in G_{M}$,\cr\cr-S({\bf r}-{\bf
t}_{g}),&$g^{\prime}\in G_{M}$,}$ (8)

the gauge transformation (7) becomes a mere shift of the origin, and the whole
description reduces to that given in the section 3.1 for periodic crystals.

We refer the reader to references in the “Further reading” section for details
regarding the enumeration of magnetic space groups of quasiperiodic crystals
and for complete lists of such groups.

##  4 Extinctions in neutron diffraction of antiferromagnetic crystals

We said earlier that every wave vector ${\bf k}$ in the lattice $L$ of a
magnetic crystal is a candidate for a diffraction peak unless symmetry forbids
it. We are now in a position to understand exactly how this happens. Given a
wave vector ${\bf k}\in L$ we examine all magnetic point-group operations $g$
or $g^{\prime}$ for which $g{\bf k}={\bf k}$. For such elements, the point-
group condition (5) can be rewritten as

$S({\bf k})=\cases{e^{2\pi i\Phi_{g}({\bf k})}S({\bf k}),&$g\in
G_{M}$,\cr\cr-e^{2\pi i\Phi_{g^{\prime}}({\bf k})}S({\bf k}),&$g^{\prime}\in
G_{M}$,}$ (9)

requiring $S({\bf k})$ to vanish unless $\Phi_{g}({\bf k})\equiv 0$, or
$\Phi_{g^{\prime}}({\bf k})\equiv 1/2$, or unless both conditions are
satisfied when both $g$ and $g^{\prime}$ are in $G_{M}$. It should be noted
that the phase values in Eq. (9), determining the extinction of $S({\bf k})$,
are independent of the choice of gauge (7), and are therefore uniquely
determined by the magnetic space-group type of the crystal.

(a) (b)

Figure 4: Extinctions of magnetic Bragg peaks in crystals with magnetic point
groups of type 1, containing time inversion $e^{\prime}$. Both figures show
the positions of expected magnetic Bragg peaks indexed by integers ranging
from -4 to 4. Filled circles indicate observed peaks and open circles indicate
positions of extinct peaks. The origin and the $D$ vectors used to generate
the patterns are indicated by an additional large circle. (a) corresponds to
the square crystal in Figure 2(a) with magnetic space group $p_{p}4mm$ and (b)
corresponds to the octagonal crystal in Figure 3(a) with magnetic space group
$p_{p}8mm$.

Particularly striking are the extinctions when the magnetic point group is of
type 1, containing time inversion $e^{\prime}$. The relation
$(e^{\prime})^{2}=e$ implies—through the group compatibility condition (6) and
the fact that $\Phi_{e}({\bf k})\equiv 0$—that $\Phi_{e^{\prime}}({\bf
k})\equiv 0$ or $1/2$. It then follows from the linearity of the phase
function that on exactly half of the wave vectors in $L$
$\Phi_{e^{\prime}}({\bf k})\equiv 1/2$ and on the remaining half, which from a
sublattice $L_{0}$ of index 2 in $L$, $\Phi_{e^{\prime}}({\bf k})\equiv 0$.
Thus, in magnetic crystals with magnetic point groups of type 1 at least half
of the diffraction peaks are missing. Figure 4 shows the positions of the
expected magnetic Bragg peaks corresponding to the square and octagonal
magnetic structures shown in Figures 2(a) and 3(a), illustrating this
phenomenon.

##  5 Generalizations of magnetic groups

There are two natural generalizations of magnetic groups. One is to color
groups with more than two colors, and the other is to spin groups where the
spins are viewed as classical axial vectors free to rotate continuously in any
direction.

An $n$-color point group $G_{C}$ is a subgroup of $O(d)\\!\times\\!S_{n}$,
where $S_{n}$ is the permutation group of $n$ colors. Elements of the color
point group are pairs $(g,\gamma)$ where $g$ is a $d$-dimensional (proper or
improper) rotation and $\gamma$ is a permutation of the $n$ colors. As before,
for $(g,\gamma)$ to be in the color point group of a finite object it must
leave it invariant, and for $(g,\gamma)$ to be in the color point group of a
crystal it must leave it indistinguishable, which in the special case of a
periodic crystal reduces to invariance to within a translation. To each
element $(g,\gamma)\in G_{C}$ corresponds a phase function
$\Phi_{g}^{\gamma}({\bf k})$, satisfying a generalized version of the group
compatibility condition (6). The color point group contains an important
subgroup of elements of the form $(e,\gamma)$ containing all the color
permutations that leave the crystal indistinguishable without requiring any
rotation $g$.

A spin point group $G_{S}$ is a subgroup of
$O(d)\\!\times\\!SO(d_{s})\\!\times\\!1^{\prime}$, where $SO(d_{s})$ is the
group of $d_{s}$-dimensional proper rotations operating on the spins, and
$1^{\prime}$ is the time inversion group as before. Note that the dimension of
the spins need not be equal to the dimension of space (e.g. one may consider a
planar arrangement of 3-dimensional spins). Also note that because the spins
are axial vectors there is no loss of generality by restricting their
rotations to being proper. Elements of the spin point group are pairs
$(g,\gamma)$ where $g$ is a $d$-dimensional (proper or improper) rotation and
$\gamma$ is a spin-space rotation possibly followed by time inversion. Here as
well, elements of the form $(e,\gamma)$ play a central role in the theory,
especially in determining the symmetry constraints imposed by the
corresponding phase functions $\Phi_{e}^{\gamma}({\bf k})$ on the patterns of
magnetic Bragg peaks, observed in elastic neutron diffraction experiments.

## Acknowledgments

I would like to thank David Mermin for teaching me the symmetry of crystals,
and Shahar Even-Dar Mandel for joining me in my investigation of the magnetic
symmetry of quasicrystals. My research in these matters is funded by the
Israel Science Foundation through grant No. 278/00.

## See also:

Crystals
Quasicrystals
Point groups
Space groups
Long-range order
Neutron diffraction
Magnetism

## Further reading:

On magnetic groups, color groups, and spin groups of periodic crystals:

  * 1.

W. Opechowski (1986), Crystallographic and metacrystallographic groups (North-
Holland, Amsterdam).

  * 2.

A. V. Shubnikov, N. V. Belov and others (1964), Colored symmetry (Pergamon
Press, New York).

  * 3.

W. Opechowski and R. Guccione (1965), in Magnetism, eds. G. T. Rado and H.
Suhl, Vol. 2A, Chap. 3 (Academic Press, New York).

  * 4.

R. L. E. Schwarzenberger (1984), “Colour symmetry,” Bull. London Math. Soc.
16, 209–240.

  * 5.

M. Senechal (1990), “Crystalline symmetries: An informal mathematical
introduction,” (Adam Hilger, Bristol).

  * 6.

D. B. Litvin and W. Opechowski (1974), “Spin groups,” Physica 76 538–554; D.
B. Litvin (1973), “Spin translation groups and neutron diffraction analysis,”
Acta Cryst. A29, 651–660; and D. B. Litvin (1977), “Spin point groups,” Acta
Cryst. A33, 279–287.

On the symmetry of quasiperiodic crystals:

  * 7.

For a review see N. D. Mermin (1992), “The space groups of icosahedral
quasicrystals and cubic, orthorhombic, monoclinic, and triclinic crystals,”
Rev. Mod. Phys. 64, 3–49.

  * 8.

For an elementary introduction see R. Lifshitz (1996), “The symmetry of
quasiperiodic crystals,” Physica A 232, 633–647.

On magnetic groups, color groups, and spin groups of quasiperiodic crystals:

  * 9.

R. Lifshitz (1997), “Theory of color symmetry for periodic and quasiperiodic
crystals,” Rev. Mod. Phys. 69, 1181.

  * 10.

R. Lifshitz and S. Even-Dar Mandel (2004), “Magnetically-ordered
quasicrystals: Enumeration of spin groups and calculation of magnetic
selection rules.” Acta Cryst. A60 167–178.