Probabilistic Cellular Automata Monte Carlo for the Maximum Clique Problem

1. Introduction

A graph is a combinatorial structure consisting of a set of objects called vertices or nodes together with a collection of links, called edges, between pairs of vertices. Graphs may be used to model the relationship between pairs of entities of a population such as those between individuals in a social network or between data points in a dataset.

A typical question in decision theory is: “Given a graph G with set of vertices V, is there a subset $A\subset V$ with cardinality k such that all pairs of vertices in A are connected by an edge?”. This problem is known as the clique problem and is known to be NP-complete (it is, actually, one of the 21 problems identified by Karp in [1]).

Its formulation as an optimization problem “what is the largest clique of a given graph G” is known as the Maximum Clique Problem and is NP-hard.

Though it is easy to verify that a collection of k nodes is indeed a clique for a graph G (where “easy” means that this verification can be done in polynomial time - this is actually the definition of an NP problem), no polynomial time algorithm is known to exist to determine the size of the largest clique of a graph. For this reason, in many situations where computation time is limited, finding an optimal clique using an exact algorithm is usually not a viable option. In these cases, having an algorithm capable of finding a good heuristic solution in a short time may be extremely valuable.

Beyond the exact algorithms, which, in the general case, are not suitable for large instances, several approaches have been used to tackle the clique problem such as greedy procedures, Monte Carlo methods, machine learning and artificial intelligence tools, and quantum computing algorithms. A recent overview can be found here [2].

Here, our focus is on Monte Carlo methods. The idea behind Monte Carlo methods consists in the ability to sample collections of nodes from a probability distribution which puts most of the weight on the largest cliques of a given graph. This probability distribution is usually approximated by the long-term distribution of a suitable Markov Chain (see [3] for an elementary introduction).

In the context of the Maximum Clique problem, Markov Chain Monte Carlo (MCMC) methods have been analyzed extensively from both the theoretical and numerical point of view. Even though in [4] some of the limitations of MCMC have been highlighted, more recent results such as those in [5,6,7] showed that this approach may be worth tackling practical instances.

This work aims to evaluate, numerically, a class of Markov Chain Monte Carlo algorithms, which here we refer to as Probabilistic Cellular Automata Monte Carlo (PCAMC), where the Markov chain lives on some multidimensional discrete space and at each step, all the components of the “current state” are updated independently one of the other. This is in contrast with the usual choice in MCMC of updating only one component at each step of the dynamics. The inherent parallel nature of PCAMC allows for a straightforward implementation that takes full advantage of parallel computing architecture giving PCAMC a significant advantage over standard MCMC methods from the computational point of view.

These algorithms have already been considered (by the author of this paper and several coauthors) to tackle certain instances of the Quadratic Unconstrained Binary Optimization Problem (QUBO) and proved to be rather effective. Here we will use them in the context of the Maximum Clique Problem and we will show that their performances are quite promising.

2. Materials and Methods

In this section, we give a mathematical definition of the maximum clique problem and see how it can be formulated in terms of the minimization of a suitable objective function. We, then, describe the Monte Carlo algorithms we want to evaluate, namely the one we refer to as the PCA (and that throughout the paper we denote by $\mathcal{A}$) and the Shaken dynamics (denoted by $\mathcal{S}$) and the algorithm we take as a reference, that is, the Metropolis “single spin flip” algorithm (that we denote by $\mathcal{M}$).

2.1. Preliminary Definitions and Statement of the Problem

Let $G=(V,E)$ be a non-oriented graph where V is the set of vertices and $E\subset V\times V$ is the set of edges. A graph G is said to be complete if $E=\left\{\right\{i,j\}\mathrm{for}\mathrm{all}i,j\in V,i\ne j\}$, that is if there as an edge between each pair of vertices of G. A graph $G^{\prime }=(V^{\prime },E^{\prime })$ is called a subgraph of G if $V^{\prime }\subset V$ and $E^{\prime }\subset E$. Let $A\subset V$ be a subset of the vertices of G. Then the graph $G\left[A\right]$ with set of vertices A and set of edges ${E|}_A=\{\{i,j\}:i,j\in A\}$ (that is the set of edges of G with both endpoints in A) is called the subgraph of Ginduced by A. The subset $C\subset V$ is called a clique for G if the induced subgraph $G\left[C\right]$ is complete. We denote by $\mathcal{K}\left(\mathcal{G}\right)$ the set of all cliques of G. A clique C is said to be maximal for G if, for all $C^{\prime }\subset V$ such that $C^{\prime }\supset C$, $C^{\prime }$ is not a clique. In words, this means that C is maximal if it can not be extended by adding to C an adjacent vertex. Further, a clique C is called a maximum clique of G if $\left|C\right|\ge |C^{\prime }|$ for all $C^{\prime }in\mathcal{K}\left(\mathcal{G}\right)$ where $\left|S\right|$ denotes the cardinality of set $\left|S\right|$. The cadinalitiy of a maximum clique of G is called the clique number of G and is denoted by $\omega \left(G\right)$...

Given a graph G, the maximum clique problem consists in determining the clique number $\omega \left(G\right)$, that is, the size of the largest clique of the graph. As mentioned in the introduction, this problem is NP-hard.

A graph $G=(V,E)$ can be encoded in several ways. One possibility is that of considering its adjacency matrix$A_G={\left(a_{ij}\right)}_{i,j=1\cdots N}$ defined by $a_{ij}=1$ if $\{i,j\}\in E$ and $a_{ij}=0$ otherwise. Note that $A_G$ is a symmetric matrix since G is non-oriented.

2.2. A Hamiltonian for the Clique Problem

A discrete optimization problem consists in finding the minimum of an objective function $H\left(\eta \right)$ where $\eta$ is an element of a suitable multidimensional space $\mathcal{X}$.

Consider a graph $G=(V,E)$, with vertices $V=\{1,2,\cdots ,N\}$ and the space of boolean configurations with N components $\mathcal{X}={\{0,1\}}^N$. Then, for each subgraph $G\left[A\right]$, with $A\subset V$, there is one, and only one, ${\eta }^A\in \mathcal{X}$ defined by ${\eta }_i^A=1\iff i\in A$. We call the set of indices where $\eta$ is positive the support of $\eta$ end denoted by $supp\left(\eta \right)$ (in words $supp\left(\eta \right)$ is the set of vertices selected by the configuration $\eta$). We can, therefore, identify subgraphs of a graph with N vertices and configurations of ${\{0,1\}}^N$ and. With an abuse of notation we write “the subgraph $\eta$‘’ in place of the subgraph $G\left[A\right]$ induced by the set $A=\{i\in 1,\cdots ,N:{\eta }_i=1\}$. Similarly, if the vertices of subgraph $\eta$ are a clique of G we refer to $\eta$ as a clique of G. Moreover we write $\left|\eta \right|$ for the cardinalitiy of a configuration $\eta$ defined as $\left|\eta \right|=\left|\right\{i\in 1,\cdots ,N:{\eta }_i=1\left\}\right|$.

For any graph $G=(V,E)$ with N vertices it is possible to define the $N\times N$ matrix $J=J\left(G\right)$ as $J_{ij}={\displaystyle \frac{1}{2}}$ if $\{i,j\}\notin E$ and $J_{ij}=0$ otherwise. We call the matrix $J\left(G\right)$ defined in this way the matrix of missing edges of G. Note that the matrix J is symmetric and is closely related to the adjacency matrix of the graph. In particular, if $A=A\left(G\right)$ is the adjacency matrix of the graph, then $J_{ij}={\displaystyle \frac{1}{2}}(1-A_{ij})$ if $i\ne j$ and 0 otherwise. Therefore, any symmetric matrix J with zeroes on the diagonal and entries in $\left\{0,{\displaystyle \frac{1}{2}}\right\}$ uniquely determines a graph $G\left(J\right)$. .Moreover, when it does not give rise to confusion we drop the explicit dependence of J on G.

Given a graph $G=(V,E)$ consider the associated matrix of missing edges J and define the function

$$\begin{array}{c}\hfill H^{J,\lambda }\left(\eta \right)=\sum _{i,j\le N}{J}_{ij}{\eta }_i{\eta }_j-\lambda \sum _{i=1}^N{\eta }_i.\end{array}$$

(1)

In the remainder of the paper, to lighten the notation, we will write H instead of $H^{J,\lambda }$.

Lemma 1.

Let $G=(V,E)$ be a non-oriented graph. If $0<\lambda <1$ then the H is minimal if and only if η is a maximum clique of G.

Proof. 

Let us first show that the condition $0<\lambda <1$ is sufficient.

Let ${\eta }_k$ be a configuration with cardinality k. Then $H\left({\eta }_k\right)=-\lambda k$ if ${\eta }_k$ and $H\left({\eta }_k\right)>-\lambda k$ otherwise since the term ${\sum }_{i,j\le N}{J}_{ij}{\eta }_i{\eta }_j>0$ for all $\eta$’s that are not cliques (since there is at least one “missing edge”).

Let $\omega \left(G\right)$ be the clique number of G, that is, the cardinality of the largest clique of G. Then there exists a clique $\overline{\eta }$ with cardinality $\omega \left(G\right)$such that $H\left(\overline{\eta }\right)=-\lambda \omega \left(G\right)$. From the previous observation, it follows that $H\left(\eta \right)\ge -\lambda \omega \left(G\right)$ for all $\eta$’s with cardinality at most $\omega \left(G\right)$ with equality only if $\eta$ is a maximum clique for G.

It is left to show that, if the cardinality of $\eta$ is larger than $\omega \left(G\right)$, then $H\left(\eta \right)>-\lambda \omega \left(G\right)$. Let $\eta$ have cardinality $k>\omega \left(G\right)$. Write $S=supp\left(eta\right)$ and denote by A the maximum clique of $G\left[S\right]$. Further let $B=S\setminus A$ and call ${\eta }^A$ and ${\eta }^B$ the configurations whose support is, respctively, A and B. Note that A and B are disjoint by construction and $\left|B\right|\ge 1$ since, by assumption $\left|A\right|+\left|B\right|>\omega \left(G\right)$ and $\left|A\right|\le \omega \left(G\right)$. We have

$$\begin{array}{c}\hfill H\left(\eta \right)=H\left({\eta }^A\right)+H\left({\eta }^B\right)+\sum _{i\in A,j\in B}{J}_{ij}{\eta }_i{\eta }_j.\end{array}$$

(2)

Now observe that, by the assumed maximality of clique A, ${\sum }_{i\in A,j\in B}{J}_{ij}{\eta }_i{\eta }_j\ge \left|B\right|$ since there must be at least one missing edge between each vertex in B and the vertices of A. Moreover, it must be $H\left({\eta }^A\right)\ge -\lambda \omega \left(G\right)$ and $H\left({\eta }^B\right)\ge -\lambda \left|B\right|$. Hence

$$\begin{array}{c}\hfill H\left(\eta \right)\ge -\lambda \omega \left(G\right)+(1-\lambda )\left|B\right|\ge -\lambda \omega \left(G\right)\end{array}$$

(3)

as soon as $\lambda <1$.

To see that the condition $\lambda <1$ is necessary, consider the case where $\eta$ is a clique of size k. Then $H\left(\eta \right)=-\lambda k$. Let A be the set of indices such that ${\eta }_i=1$ (that is, the set of vertices of the clique). Consider the configuration $\tau$ obtained from $\eta$ by setting ${\eta }_j=1$ for some $j\notin A$ and suppose that there is an edge between vertex j and $k-1$ vertices in A. Then $H\left(\tau \right)=-\lambda (k+1)+1$. Then $H\left(\tau \right)\le H\left(\eta \right)$ as soon as $\lambda >=1$.

Finally observe that if $\lambda =0$, $H\left(\eta \right)=0$ for all $\eta$’s that are a clique and if $\lambda <0$ only the empty configuration $\overrightarrow{0}$ has energy $H\left(\overrightarrow{0}\right)=0$ whereas $H\left(\eta \right)>0$ for all $\eta \ne \overrightarrow{0}$.

   □

Thanks to the previous lemma, the maximum clique problem can be formulated in a standard optimization problem form as

$$\begin{array}{c}\hfill \underset{\eta \in {\{0,1\}}^N}{min}H\left(\eta \right)=\underset{\eta \in {\{0,1\}}^N}{min}\sum _{i,j\le N}{J}_{ij}{\eta }_i{\eta }_j-\lambda \sum _{i=1}^N{\eta }_i,0<\lambda <1.\end{array}$$

(4)

Note that, exploiting the fact that ${\eta }_i={\eta }_i^2$, It is possible to rewrite the previous formula in terms of a matrix $\tilde{J}=\tilde{J}\left(\lambda \right)$ as

$$\begin{array}{c}\hfill \underset{\eta \in {\{0,1\}}^N}{min}H\left(\eta \right)=\underset{\eta \in {\{0,1\}}^N}{min}\sum _{ij\le N}{\tilde{J}}_{ij}{\eta }_i{\eta }_j=\underset{\eta \in {\{0,1\}}^N}{min}{\eta }^T\tilde{J}\eta \end{array}$$

(5)

with ${\tilde{J}}_{ij}=J_{ij}$ if $i\ne j$ and ${\tilde{J}}_{ii}=-\lambda$. This formulation is known in the literature as Quadratic Unconstrained Binary Optimization Problem (QUBO): an NP-hard problem with a vast range of applications (see the text [8] for an overview) that allows for a straightforward embedding of (beyond the maximum clique problem) of several combinatorial optimization problems such as max-cut and graph coloring (see [9]. Further, QUBO is a model of interest in the field of quantum computing (see, e.g., [10,11,12]),

It is worth noting that the objective function of an optimization problem can be interpreted as the Hamiltonian energy function of a statistical mechanics model with state space $\mathcal{X}$. The statistical mechanics model’s ground states (minima of the Hamiltonian) correspond to configurations minimizing the objective function. The equilibrium distribution of a statistical mechanics model with Hamiltonian H is described by the Gibbs measure

$$\begin{array}{c}\hfill {\mu }_G\left(\eta \right)={\displaystyle \frac{1}{Z}}{e}^{-\beta H\left(\eta \right)},\eta \in \mathcal{X}\end{array}$$

(6)

where $\beta$ is a positive real parameter called the inverse temperature and Z is a normalizing constant called the partition function. As $\beta \to \infty$ (low-temperature regime) this probability measure concentrates on the ground states of the system. This fact suggests that a solution to the optimization problem can be found by sampling from the Gibbs measure for the corresponding statistical mechanics model at low temperature.

To this aim, it is possible to define a Markov chain on $\mathcal{X}$ with transition probability $P:\mathcal{X}\times \mathcal{X}\to [0,1]$ whose stationary distribution is the Gibbs measure ${\mu }_g$. If the chain is irreducible and aperiodic, then, irrespective of the initial starting configuration, the probability distribution of the chain after n steps tends to the stationary distribution as $n\to \infty$. The sampling from the Gibbs measure can, then, be performed by letting the chain evolve for a sufficiently long time. More precisely, let ${\mu }_n\left(\eta \right)$ be the probability distribution after n steps of the chain starting from configuration $\eta$ and let $\pi$ be the stationary measure for the chain with transition probability P. Then, denoting by d_TV$(·,·)$ the total variation distance.

$$\begin{array}{c}\hfill {\mathrm{d}}_{TV}({\mu }_n\left(\eta \right),\pi )\to 0\mathrm{as}n\to \infty .\end{array}$$

(7)

The time it takes for the chain to have a distribution that is within a distance $\epsilon$ from the equilibrium distribution is referred to as the mixing time of the chain. In formulae, the mixing time of a Markov Chain with state space $\mathcal{X}$ and stationary distribution $\pi$ is defined, for $\epsilon <{\displaystyle \frac{1}{2}}$, as

$$\begin{array}{c}\hfill t_{\mathrm{mix}}\left(\epsilon \right)=min\{n:{\mathrm{d}}_{TV}({\mu }_n\left(\eta \right),\pi )\le \epsilon \}\end{array}$$

(8)

In general, determining a priori bound on the mixing time of the chain that guarantees that the algorithm is useful in practice (that is it has a mixing time that is at most polynomial in the size of the problem) may be non-trivial. Actually, in the case of Erdős-Rényi random graphs, Jerrum ([4]) proved that the mixing time of the Metropolis process he considered for the Maximum Clique Problem has super-polynomial mixing time. However, for instances that are not too big (number of nodes up to a few thousand) Monte Carlo algorithms may still provide some good heuristic solutions (see, e.g., [5].

2.3. The Metropolis Single Spin Flip Algorithm

A reference choice among Markov Chain Monte Carlo algorithms is the Metropolis update rule with single spin flip whose transition probability matrix is as follows.

For a given graph $G=(V,E)$ with missing edges matrix J, consider the Hamiltonian function H defined in (1). Let ${\eta }^{\left(i\right)}$ be the configuration obtained from $\eta$ by flipping the occupation number of the i-th component, that is ${\eta }_i^{\left(i\right)}=1-{\eta }_i$ and ${\eta }_j^{\left(i\right)}={\eta }_j$ for all $j\ne i$.

$$\begin{array}{c}\hfill P_{\mathcal{M}}(\eta ,\tau )=\left\{\begin{array}{cc}{\displaystyle \frac{1}{\mathrm{N}}}{e}^{-\beta {[H\left(\tau \right)-H\left(\eta \right)]}_{+}}.\hfill & \mathrm{if}\tau ={\eta }^{\left(i\right)},i=1\cdots N\hfill \\ 1-\sum _i{\displaystyle \frac{1}{\mathrm{N}}}{e}^{-\beta {[H\left({\eta }^{\left(i\right)}\right)-H\left(\eta \right)]}_{+}}\hfill & \mathrm{if}\tau =\eta \hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}\right.\end{array}$$

(9)

where ${\left[·\right]}_{+}$ denotes the positive part. In words, at each step, an index $i\in 1,\cdots ,N$ is chosen a random and if the configuration obtained by flipping the value of the i-th component of $\eta$ has a Hamiltonian energy lower than that of $\eta$ the occupation number of ${\eta }_i$ is changed with probability one. If this is not the case, the occupation number of ${\eta }_i$ is changed with a probability that is exponentially small in the difference between the energy of the target configuration and the energy of the current one.

Let ${\pi }_{\mathcal{M}}\left(\eta \right)={\displaystyle \frac{1}{Z}}{e}^{-\beta H\left(\eta \right)}$. Then it is immediate to verify that $P_{\mathcal{M}}$ satisfies the detailed balance condition

$${\pi }_{\mathcal{M}}\left(\eta \right)P_{\mathcal{M}}(\eta ,\tau )={\pi }_{\mathcal{M}}\left(\tau \right)P_{\mathcal{M}}(\tau ,\eta )$$

ensuring that ${\pi }_{\mathcal{M}}$ is the stationary distribution of $P_{\mathcal{M}}$.

We remark that this algorithm is allowed to visit configurations that are not cliques of the graph G. However, in the limit $\beta \to \infty$, if $\lambda \propto {\displaystyle \frac{1}{\beta }}$ the algorithm is equivalent (see [5]) to the Metropolis algorithm considered by Jerrum in [4] which allows transitions only between cliques of G.

2.4. The “Pair Hamiltonian” Probabilistic Cellular Automton

A probabilistic cellular automaton on $\mathcal{X}$ is a Markov Chain whose transition probability matrix can be written, for any $\eta ,\tau \in \mathcal{X}$ as

$$\begin{array}{c}\hfill P_{\eta ,\tau }=\prod _{i}P({\tau }_i=·|\eta ).\end{array}$$

(10)

In words, this means that at each step, the value of each component ${\tau }_i$ of the target configuration $\tau$ can be sampled independently of the values of all other ${\tau }_j,j\ne i$. From a computational point of view this fact is particularly appealing, since, in principle, the update value of each component of $\tau$ could be demanded to a dedicated computing core (see [13]). With the advent of massively parallel processors such as GPUs, this is an actual possibility even on consumer hardware for configurations with several hundreds of components.

In [14] the PCA Monte Carlo algorithm has been introduced to study the minima of QUBO. When used to approach the maximum clique problem the algorithm is as follows.

Consider the pair Hamiltonian function defined on pairs of configurations

$$\begin{array}{c}\hfill H(\eta ,\tau )=-\left(\beta \sum _{i,j}{\tilde{J}}_{ij}{\eta }_i{\tau }_j+q\sum _i[{\eta }_i(1-{\tau }_i)+{\tau }_i(1-{\eta }_i)]\right),\end{array}$$

(11)

where $\tilde{J}$ is the matrix used in (5), and set the transition probability to be

$$\begin{array}{c}\hfill P_{\mathcal{A}}(\eta ,\tau )=e^{-H(\eta ,\tau )}={\displaystyle \frac{e^{-H(\eta ,\tau )}}{{\sum }_{\tau }{e}^{-H(\eta ,\tau )}}}\end{array}$$

(12)

Let $h_i\left(\eta \right)={\sum }_j{\tilde{J}}_{ij}{\eta }_j$ be the local energy field felt by the i-th component of $\eta$. Then $H(\eta ,\tau )$ can be rewritten as $H(\eta ,\tau )=\beta {\sum }_i{h}_i\left(\eta \right){\tau }_i+q\left([{\eta }_i(1-{\tau }_i)+{\tau }_i(1-{\eta }_i)]\right)$. Consequently, the transition probabilities can be rewritten as

$$\begin{array}{c}\hfill P_{\mathcal{A}}(\eta ,\tau )=\prod _i{\displaystyle \frac{e^{-\beta h_i\left(\eta \right){\tau }_i-q[{\eta }_i(1-{\tau }_i)+{\tau }_i(1-{\eta }_i)]}}{Z_i}}\end{array}$$

(13)

yielding

$$\begin{array}{c}\hfill P({\tau }_i=1|\eta )={\displaystyle \frac{e^{-\beta h_i\left(\eta \right)-q(1-{\eta }_i)}}{Z_i}}\mathrm{and}P({\tau }_i=0|\eta )={\displaystyle \frac{e^{-q{\eta }_i}}{Z_i}}\end{array}$$

(14)

with $Z_i=e^{-\beta h_i-q(1-{\eta }_i)}+e^{-q{\eta }_i}$. Therefore $P_{\mathcal{X}}$ defines, indeed, a PCA in the sense of (10).

By the symmetry of $\tilde{J}$, a straightforward computation shows that the detailed balance condition

$$\begin{array}{c}\hfill P_{\mathcal{A}}(\eta ,\tau ){\displaystyle \frac{{\sum }_{\tau }{e}^{-H(\eta ,\tau )}}{{\sum }_{\eta ,\tau }{e}^{-H(\eta ,\tau )}}}=P_{\mathcal{A}}(\tau ,\eta ){\displaystyle \frac{{\sum }_{\eta }{e}^{-H(\tau ,\eta )}}{{\sum }_{\eta ,\tau }{e}^{-H(\eta ,\tau )}}}\end{array}$$

(15)

holds and, hence,

$$\begin{array}{c}\hfill {\pi }_{\mathcal{A}}\left(\eta \right)={\displaystyle \frac{{\sum }_{\tau }{e}^{-H(\eta ,\tau )}}{{\sum }_{\eta ,\tau }{e}^{-H(\eta ,\tau )}}}\end{array}$$

(16)

is the stationary measure of $P_{\mathcal{A}}$. Note that also this algorithm may visit all configurations and not only those that are cliques.

A closer look at equation (13) shows that at each step, the dynamics tends to select, for the target configuration $\tau$, the vertices connected to all vertices selected by configuration $\eta$. Indeed, for these components, the local energy field is negative (equal to $-\lambda$). On the other hand, the probability of selecting, in $\tau$, components for which there is at least one missing edge with the components selected by $\eta$ is exponentially small (since $\lambda <1$). Further note that the term $e^{-q\left[{\eta }_i(1-{\tau }_i+{\tau }_i(1-{\eta }_i))\right]}$ reduces the probability that the occupation number of the i-th component of $\tau$ differs from the i-th component of $\eta$ so that the probability that $\eta$ and $\tau$ differ at too many sites stays small.

To get an idea of why this algorithm is expected to work, it is possible to argue as follows. At first observe that $H(\eta ,\eta )=\beta H\left(\eta \right)$ with H defined as in LABEL_ABOVE. Further note that, as q gets large, the weight of pairs $(\eta ,\tau )$ with $\eta \ne \tau$ in ${\pi }_{\mathcal{A}}\left(\eta \right)$ is exponentially depressed ${\pi }_{\mathcal{A}}\left(\eta \right)$ becomes closer to the Gibbs measure ${\displaystyle \frac{e^{-\beta H\left(\eta \right)}}{Z}}$ with Hamiltonian H and inverse temperature $\beta$.

In [14] and in [15] the authors considered QUBO instances where the coefficients of the matrix $\tilde{J}$ are realization of independent identically distributed random variables with expected value zero and showed that the PCA Monte Carlo performed quite well when compared to other heuristic techniques. Here, however, the elements of the matrix $\tilde{J}$ have a particular structure: they are negative on the diagonal and positive or null outside of it and it is not, a priori, obvious how this structure could have an impact on the practical performances of the algorithm.

2.5. The Shaken Dynamics

Shaken dynamics have been introduced in [16,17] in the process of finding efficient algorithms to sample from the Gibbs measure on spin systems (see [18,19]) and extend the class of PCA with transition probabilities of type $P_{\mathcal{A}}={\displaystyle \frac{e^{-\beta H(\eta ,\tau )}}{{\sum }_{\tau }{e}^{-\beta H(\eta ,\tau )}}}$ described above and defined in terms of a pair Hamiltonian.

The transition probabilities of the Shaken dynamics come as a combination of two steps and are of the following type:

$$\begin{array}{c}\hfill P_{\mathcal{S}}(\eta ,\tau )=\sum _{\sigma }\overleftarrow{P}(\eta ,\sigma )\overrightarrow{P}(\sigma ,\tau )\end{array}$$

(17)

with

$$\begin{array}{c}\hfill \overleftarrow{P}(\eta ,\sigma )={\displaystyle \frac{e^{-\overleftarrow{H}(\eta ,\sigma )}}{{\sum }_{\sigma }{e}^{-\overleftarrow{H}(\eta ,\sigma )}}}\mathrm{and}\overrightarrow{P}(\eta ,\sigma )={\displaystyle \frac{e^{-\overrightarrow{H}(\eta ,\sigma )}}{{\sum }_{\sigma }{e}^{-\overrightarrow{H}(\eta ,\sigma )}}}.\end{array}$$

(18)

The two functions $\overleftarrow{H}$ and $\overrightarrow{H}$ appearing in the definition of the transition probabilities must be such that $\overleftarrow{H}(\eta ,\sigma )=\overrightarrow{H}(\sigma ,\eta )$. If this is the case, then it follows that

$$\begin{array}{c}\hfill {\pi }_{\mathcal{S}}\left(\eta \right)={\displaystyle \frac{{\sum }_{\sigma }{e}^{-\overleftarrow{H}(\eta ,\sigma )}}{{\sum }_{\eta ,\sigma }{e}^{-\overleftarrow{H}(\eta ,\sigma )}}}\end{array}$$

(19)

is the stationary (reversible) measure for $P_{\mathcal{S}}$. To see this, observe that the detailed balance condition holds. Indeed, calling $Z_{\eta }={\sum }_{\sigma }{e}^{-\overleftarrow{H}(\eta ,\sigma )}$ and $Z={\sum }_{\eta ,\sigma }{e}^{-\overleftarrow{H}(\eta ,\sigma )}$, we have

$$\begin{array}{cc}\hfill {\pi }_{\mathcal{S}}\left(\eta \right)P_{\mathcal{S}}(\eta ,\tau )& ={\displaystyle \frac{{\sum }_{\sigma }{e}^{-\overleftarrow{H}(\eta ,\sigma )}}{{\sum }_{\eta ,\sigma }{e}^{-\overleftarrow{H}(\eta ,\sigma )}}}\sum _{\sigma }{\displaystyle \frac{e^{-\overleftarrow{H}(\eta ,\sigma )}}{{\sum }_{\sigma }{e}^{-\overleftarrow{H}(\eta ,\sigma )}}}{\displaystyle \frac{e^{-\overrightarrow{H}(\sigma ,\tau )}}{{\sum }_{\tau }{e}^{-\overrightarrow{H}(\sigma ,\tau )}}}\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill & ={\displaystyle \frac{{\sum }_{\sigma }{e}^{-\overleftarrow{H}(\eta ,\sigma )}}{{\sum }_{\eta ,\sigma }{e}^{-\overleftarrow{H}(\eta ,\sigma )}}}\sum _{\sigma }{\displaystyle \frac{e^{-\overrightarrow{H}(\sigma ,\eta )}}{{\sum }_{\sigma }{e}^{-\overleftarrow{H}(\eta ,\sigma )}}}{\displaystyle \frac{e^{-\overleftarrow{H}(\tau ,\sigma )}}{{\sum }_{\tau }{e}^{-\overrightarrow{H}(\sigma ,\tau )}}}\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill & ={\displaystyle \frac{{\sum }_{\tau }{e}^{-\overleftarrow{H}(\tau ,\sigma )}}{{\sum }_{\tau ,\sigma }{e}^{-\overleftarrow{H}(\tau ,\sigma )}}}\sum _{\sigma }{\displaystyle \frac{e^{-\overleftarrow{H}(\tau ,\sigma )}}{{\sum }_{\sigma }{e}^{-\overleftarrow{H}(\tau ,\sigma )}}}{\displaystyle \frac{e^{-\overrightarrow{H}(\sigma ,\eta )}}{{\sum }_{\sigma }{e}^{-\overrightarrow{H}(\sigma ,\eta )}}}={\pi }_{\mathcal{S}}\left(\tau \right)P_{\mathcal{S}}(\tau ,\sigma )\hfill \end{array}$$

(22)

Having defined the Shaken dynamics, we can describe how it can be used as a heuristic algorithm for the maximum clique problem.

Consider a graph $G=(V,E)$, with $\left|V\right|=N$, its associated missing edge matrix and J and the corresponding matrix $\tilde{J}=\tilde{J}\left(\lambda \right)$ used to formulate the maximum clique problem as QUBO (as in (5)).

Let $\overleftarrow{J}$ be such that $\overleftarrow{J}+\overleftarrow{J^T}=\tilde{J}$ and define

$$\begin{array}{c}\overleftarrow{H}(\eta ,\sigma )=-\beta \sum _{i,j\le N}{\overleftarrow{J}}_{ij}{\eta }_i{\sigma }_j-q\sum _i[{\eta }_i(1-{\sigma }_i)+{\sigma }_i(1-{\eta }_i)]\end{array}$$

(23)

$$\begin{array}{c}\overrightarrow{H}(\eta ,\sigma )=-\beta \sum _{i,j\le N}{\overleftarrow{J^T}}_{ij}{\eta }_i{\sigma }_j-q\sum _i[{\eta }_i(1-{\sigma }_i)+{\sigma }_i(1-{\eta }_i)].\end{array}$$

(24)

Then the condition $\overleftarrow{H}(\eta ,\sigma )=\overrightarrow{H}(\sigma ,\eta )$ is verified. Further, observe that

$$\begin{array}{c}\hfill \overleftarrow{H}(\eta ,\eta )={\eta }^T\tilde{J}\eta ={\eta }^T(\overleftarrow{J}+\overleftarrow{J^T})\eta =2\left({\eta }^T\overleftarrow{J}\eta \right)={\displaystyle \frac{1}{2}}\beta H\left(\eta \right)\end{array}$$

(25)

where $H\left(\eta \right)$ is the same the one defined in (5).

Strictly speaking, the Shaken dynamics is not a PCA. However the transition probabilities of each half-step are of the form (13) where, instead of the vector of field h, appear, alternatively, the two vectors of fields $\overleftarrow{h}=\overleftarrow{J}·\eta$ and $\overrightarrow{h}=\overleftarrow{J^T}·\sigma$.

Also in this case, if the parameter q is large, the stationary measure ${\pi }_S\left(\eta \right)$ is close to the Gibbs measure ${\displaystyle \frac{e^{-\beta H\left(\eta \right)}}{Z}}$ (the factor $e^{-1/2}$ appears both at the numerator and at the denominator and, therefore, cancels out) and, hence, the Shaken dynamics provides a heuristic algorithm for the maximum clique problem. Also this Markov chain lives on the whole set $\mathcal{X}$ and not only on the subsets corresponding to cliques.

In [17] the Shaken dynamics has been used, in the context of spin systems, to find the minima of a QUBO-like problem defined on a lattice that is a problem where the interaction $J_{ij}$ was different from zero only for pairs $i,j$ satisfying a certain relation. In that scenario, the Shaken dynamics appeared to be rather effective. However, in the case of the maximum clique problem, besides the requirement for J to be symmetric, there is no restriction on which (and how many!) pairs $J_{ij}$ can be different from zero as long as $i\ne j$. Therefore the effectiveness of this approach to the maximum clique problem has to be assessed.

3. Results

We tested our algorithms on a set of benchmark graphs commonly considered in the literature (DIMACS benchmarks). The graph considered have a different number of nodes N spanning from 125 to 4000. The graphs taken into account are of several types which can be retrieved from the graph id:

Erdős-Rényi graphs:

graphs with a fixed number of vertices and edges selected independently at random with uniform probability. These are graphs with id of type Cn.d (where n is the number of vertices and d/10 is the density (probability) of edges) and DSJCn_d (where n is the number of nodes and d/10 is the density of edges)

Steiner triple graphs:

with id MANN_aXX are graphs constructed to have a clique formulation of the Steiner Triple Problem

Brockington graphs:

with id brockN_d where N is the number of vertices and d is a parameter

Sanchis graphs:

with id genn_p0.x_y where n is the number of vertices, p0.x is the density of edges and y is the size of the planted clique

Hamming graphs:

with id hammingx-y with parameters x and y

Keller graphs:

with id kellern and parameters n = 4, 5, 6

P-hat graphs:

with id p-hatn-x where n is the number of nodes and x is an identifier; graphs generated with p-hat have wider node degree spread and larger cliques than uniform graphs.

More details can be found in [20].

We formulated the problem as in (5) and fixed, for all graphs and all algorithms, $\lambda =0.25$.

As far as the parameters $\beta$ and q, since we do not have an apriori argument for an optimal choice, we tested several pairs of values. However, to determine the values of q to perform our test, we kept into account the size N (number of nodes) of the graph and chose values of q so as to have $e^{-q}N=c_q\sqrt{N}$ for several values of $c_q$ as described below. The rationale for this choice is the fact that, for $\mathcal{A}$ and S, the term $e^{-q}$ is proportional to the probability of “flipping” a component neglecting the contribution depending on the “energy field” for that component. In this way we intended to change, at each iteration, a number of component of the order$\sqrt{N}$.

The values we used are the following:

Metropolis ($\mathcal{M}$):

$\beta$ in the range $1,1.5,2,\cdots ,4$

PCA ($\mathcal{A}$):

$\beta$ in the range $1,1.5,2,\cdots 4$ and $c_q$ in the set ${\displaystyle \frac{1}{8}},{\displaystyle \frac{1}{4}},{\displaystyle \frac{1}{2}},1,2,3,4$

Shaken dynamics ($\mathcal{S}$):

in the range $1,2,\cdots 7$ and $c_q$ in the set $\{{\displaystyle \frac{1}{64}},{\displaystyle \frac{1}{32}},{\displaystyle \frac{1}{16}},{\displaystyle \frac{1}{8}},{\displaystyle \frac{1}{4}},{\displaystyle \frac{1}{2}},{\displaystyle \frac{1}{1}}\}$

These sets of ranges for the parameters have been determined after some experimentation where we saw that certain ranges of parameters were not suitable for a given algorithm. In particular, we needed to find values of q big enough (that is $c_q$ sufficiently small) for the dynamics to retain a few components. Note that not all pairs of parameters $(\beta ,q)$ that have been taken into account are effective for a given algorithm.

For the Shaken dynamics $\mathcal{S}$, we also had to make a choice for the matrix $\overleftarrow{J}$. We proceeded as follows. For each $i\in 1,\cdots ,N$ we set

$$\begin{array}{c}\hfill {\overleftarrow{J}}_{ij}=\left\{\begin{array}{cc}{J}_{i{j}^{\prime }}\hfill & \mathrm{with}{j}^{\prime }=((j-1)modN)+1\mathrm{for}j\in (i+1,i+2,\cdots i+⌈{\displaystyle \frac{N}{2}}⌉-1)\hfill \\ {\displaystyle \frac{J_{ij}}{2}}\hfill & \mathrm{for}j=i+{\displaystyle \frac{N}{2}}\mathrm{if}\mathrm{N}\mathrm{is}\mathrm{even}\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}\right.\end{array}$$

(26)

Intuitively, for each row, the matrix $\overleftarrow{J}$ is obtained retaining the ${\displaystyle \frac{N}{2}}$ elements of $\tilde{J}$ “on the right” of the diagonal element. It is straightforward to verify that $\overleftarrow{J}+\overleftarrow{J^T}=\tilde{J}$.

This choice also provides a reason for why the values $\beta$ used for the algorithm $\mathcal{S}$ are in a range that is, roughly, twice as wide then that used for $\mathcal{M}$ and A: the fields used to update the configuration at each “half-step” are proportional to $\overleftarrow{J}·\eta$ whose components could be expected to be approximately one-half of the components of $J·eta$ (which are the fields determining the transitions of $\mathcal{A}$).

Simulations were run for a fixed number of iterations and the largest clique found in each run of any of the algorithms was recorded. For each set of parameters ($\beta$ for $\mathcal{M}$ and $(\beta ,q)$ for $\mathcal{A}$ and $\mathcal{S}$) we run 10 simulations. The number of iterations has been set in the following way:

Metropolis ($\mathcal{M}$):

100000 sweeps where a sweep corresponds to N steps of the Markov Chain (that is a total of $100000N$ “attempted” flips);

PCA ($\mathcal{A}$):

200000 steps of the Markov Chain;

Shaken dynamics ($\mathcal{S}$):

100000 complete steps of the Markov Chain (200000 half-steps).

Our intent was to perform a fair comparison of the three algorithms. However it is not obvious what the right criterion to establish fairness should be. For instance, one could run the chains for the same number of “nominal” steps. However, this approach would be rather penalizing for the Metropolis algorithm since it would take, on average, $O(NlogN)$ steps to give all vertices of the graphs at least one opportunity to be selected. Another opportunity would be to count the number of “attempted flips”. Since in both $\mathcal{A}$ and $\mathcal{S}$ all components are potentially updated, following this approach would require setting the number of steps for both $\mathcal{A}$ and $\mathcal{S}$ to be the same of the number of sweeps for $\mathcal{M}$. However, if run on a single core, the computing time for a step of $\mathcal{A}$ of $\mathcal{S}$ is significantly shorter than the time required for a sweep of $\mathcal{M}$ (see Table 1). Finally one could simply consider computation time. However, we feel that this criterion would be too dependent on the implementation details and the features of the computing machine. Because of these considerations, we decided to use a “hybrid” approach using as a primary criterion the total number of attempted flips but giving an extra quota of iterations to the parallel dynamics (note that, in each full step, the Shaken dynamics tries to update all component twice) to take into account (in a conservative way) their faster execution times.

The three algorithms were implemented in Julia (version 1.10, see [21]) and executed on an Nvidia DGX-1 system equipped with Intel Xeon CPU E5-2698 v4 @ 2.20GHz (80 cores). We considered only a CPU implementation of the algorithms. For linear algebra operations (matrix-vector products) we used BLAS libraries ([22,23]) with the default implementation shipped with the LinearAlgegra package of Julia. Random numbers were generated using the default random number generator available in Julia (which in version 1.10 implements the Xoshiro256++ algorithm [24]).

The Metropolis algorithm finds the best-known clique 30 out of 37 times proving to be, in practice, a rather effective tool for the maximum clique problem for instances with a size that is not too large. The PCA algorithm $\mathcal{A}$ allowed to find the maximum clique 27 out of 37 times whereas the Shaken dynamics found the maximum clique 23 out of 37 times. In those cases where neither of the algorithms found the maximum clique the Metropolis is the one which found the largest clique more often. It should be noted, though, that in the case of the graph MANN_a81, the PCA found a better solution of the Metropolis and in two more cases it found an equivalent solution.

In general, the performances of the two parallel algorithms, especially those of the PCA $\mathcal{A}$, are not much worse in terms of the value of the largest clique found. Nevertheless, it is possible to argue that the slightly worse performances of $\mathcal{A}$ and $\mathcal{S}$ are partly compensated by their computational effectiveness. Even using a single thread for the computation of the energy field, algorithm $\mathcal{A}$ was faster between 5 and 40 times. Running with 8 threads the speedup reached a factor beyond 100 for certain instances and it was more than 50 the majority of times. The situation is not too different for the Shaken dynamics $\mathcal{S}$. Here one has to take into account that, at each step, two configurations are visited: one for each half-step. A GPU implementation of the two algorithms is expected to yield an even greater speed up, especially on the larger instances. The increased speed could be advantageous in applications.

Looking at the two parallel algorithms, $\mathcal{A}$ appears to be better than $\mathcal{S}$ both in terms of the value of the solution provided and in terms of speed (even though the Shaken dynamics explores two configurations at each step the time to get a new configuration with the PCA appears to be less than ${\displaystyle \frac{1}{2}}$ of the time required by the shaken dynamics). Likely this can be explained by the fact that the energy fields of the PCA can be computed as the $\tilde{J}·\eta$ matrix-vector product where the matrix $\tilde{J}$ is symmetric, whereas the matrix $\overleftarrow{J}$ used to compute $\overleftarrow{J}·\eta$ is not.

4. Discussion

We believe that the results presented above foster a further investigation into the behavior of Probabilistic Cellular Automata. Indeed, they appear to provide good results while taking great advantage of the features of the computing architectures available nowadays. This investigation should start, in our opinion, from a better understanding of the relation between the stationary measure of the PCA and the Gibbs measure. In particular, it would be interesting to determine, for any $\overline{\beta }$, bounds for the parameters $(\beta ,q)$ so that the total variation distance between ${\pi }_{\mathcal{A}}$ (or ${\pi }_{\mathcal{S}}$) and the Gibbs measure at inverse temperature $\overline{\beta }$ is within a prescribed value $\epsilon$.

In the previous section, we observed that the values of the largest clique obtained with the Metropolis algorithm $\mathcal{M}$ are slightly better than those obtained with $\mathcal{A}$ and $\mathcal{S}$. It should be mentioned that the dependence on two parameters of the latter algorithms makes their tuning more challenging than the tuning of the single parameter of the Metropolis. For this reason, it would be interesting to study, from a theoretical point of view, the behavior of both dynamics in the $(\beta ,q)$ plane. In particular, it would be interesting to determine a “phase transition” curve identifying a “high temperature” region where the stationary measure is “close” to the uniform measure on $\mathcal{X}$ and a “low temperature” where the stationary measure is “concentrated” on the minima of the Hamiltonian 1. It is to be expected that this phase transition curve will depend on the law with which the graph has been generated An analog study would be interesting for the Metropolis case where it would be interesting to determine a critical value of $\beta$ separating the high and low-temperature regions.

As already mentioned, the Metropolis algorithm discussed above is equivalent, at large inverse temperature, to the algorithm described by Jerrum in [4]. In that paper, it has been shown that the mixing time of the chain, in the case of Erdős-Rényi graphs, is not polynomial in the number of vertices. Further, the limiting typical size of clique typically determined by the Metropolis process has been determined. It would be interesting to perform a similar analysis for the two parallel algorithms introduced here.

From the computational point of view, the most expensive task of the PCA and the Shaken dynamics is the computation of the vector of fields which is of the type $h\left(\eta \right)=J·\eta$. This matrix-vector product can be performed using highly optimized linear algebras libraries, which may take great advantage of multicore processors and GPUs.

As a consequence, the algorithm will benefit for free from improvements in the performances of linear algebra libraries. Since the use of these libraries is ubiquitous in computer applications, improvements in this direction will be driven by the efforts of communities coming from the most diverse domains.

Note this type of advantage could be exploited even further. Indeed, it is possible to let k PCA evolve in parallel. Denoting by $\overrightarrow{E}$ the matrix whose columns are the configurations of the k chains, the collection of the k vector of fields${\left\{h_i\left({\eta }^{\left(m\right)}\right)\right\}}_{i=1,\cdots ,N;m=1,\cdots ,k}$ can be computed as the matrix matrix product $\tilde{J}·\overrightarrow{E}$. Performing this operation is, generally, significantly faster than the computation of k matrix-vector products. Note that, for each chain, a different pair of parameters $\beta$ and q could be used without significantly affecting the computation time.

The possibility to run multiple chains at the time with different inverse temperatures $\beta$ makes it possible to consider both dynamics $\mathcal{A}$ and $\mathcal{S}$ for algorithms like parallel tempering [25] which has been applied to the clique problem using standard single spin Metropolis dynamics (see, e.g. [7]).

Note that Markov Chain Monte Carlo algorithms defined in terms of a “pair Hamiltonian” were already considered in [5,26,27]. Also, the algorithm considered in those works (referred to as “Cavity”) had the feature of updating more than one component of the configuration at each step. However, in these works, the cardinality of the configuration stayed fixed throughout the whole evolution of the chain and, therefore, its components could not be updated independently limiting, in this way, the possibility to exploit parallel computing architectures at their fullest.

As a final comment, we would like to mention that both the parallel algorithms discussed in this paper may have applications outside combinatorial optimization. As far as the PCA is concerned, it has been suggested in [28] that this Markov Chain could be used to sample Exponential Random Graphs in the generation of Synthetic Power Grids. For what concern the shaken dynamics, it has been hinted in [17,29] has a microscopic for tidal dissipation.