2.1. Graphical framework

In this paper, we consider connected, edge-weighted, undirected graphs $G = (V, E, \omega )$ with a node set $V = \{1, \ldots, |V|\}$ , an edge set $E=\{ (i,j) \}^{|V|}_{i,j=1}$ , and edge weight $\omega _{ij}$ between node $i$ and node $j$ . The (weighted) adjacency matrix $W=( \omega _{ij})^{|V|}_{i,j=1}$ is a symmetric matrix whose entries $\omega _{ij}$ are zero if $i=j$ or $(i,j) \notin E$ , and positiveFootnote ² otherwise. Consequently, we can consider a given adjacency matrix $W$ (with non-negative entries and zeros on the diagonal) as defining the edge structure of the graph. Graphs with non-negative edge weights are known as unsigned graphs. We consider unweighted graphs as special cases of weighted graphs for which, for all $i,j\in V$ , $\omega _{ij}\in \{0,1\}$ .

We will reserve the notation $W$ for the adjacency matrix of the connected graph whose nodes we wish to cluster. Along the way we also encounter graphs defined by other adjacency matrices; therefore in our definitions of adjacency-matrix-dependent quantities, we will use the (dummy) symmetric matrix $C\in [0,\infty )^{|V|\times |V|}$ with entries $c_{ij}$ . We note that we do not require the diagonal entries $c_{ii}$ to be zero, as is sometimes required in other sources (i.e., we allow the graph defined by $C$ to have self-loops).

We denote the set of neighbours of node $i\in V$ with respect to the adjacency matrix $C$ by:

\begin{equation*} \mathcal {N}_C(i) \;:\!=\; \{j\in V \;:\; c_{ij} \gt 0\}. \end{equation*}

The (weighted) degree of node $i$ with respect to the adjacency matrix $C$ is

(1) \begin{equation} (d_C)_i\;:\!=\; \sum _{j\in V} c_{ij}. \end{equation}

The diagonal degree matrix $D_C$ of $C$ has entries $\left (D_C \right )_{ii}=(d_C)_i$ . The maximum and minimum node degrees with respect to $C$ are

\begin{equation*} d_{C,\textrm {max}} \;:\!=\; \max _{i\in V} (d_C)_i \quad \text {and} \quad d_{C,\textrm {min}} \;:\!=\; \min _{i\in V} (d_C)_i, \end{equation*}

respectively. The volume (i.e., total edge weight) of a subset $S\subset V$ with respect to $C$ is defined as:

\begin{equation*} \textrm {vol}_C(S)\;:\!=\;\sum _{i\in S} (d_C)_i = \sum _{i \in S, j \in V} c_{ij}. \end{equation*}

In particular, $\textrm{vol}_C(V)$ is also called the volume of the graph with adjacency matrix $C$ .

We define the set of real-valued node functions:

\begin{equation*} \mathcal {V} \;:\!=\; \left \{u\;:\; V \to \mathbb {R}\right \}. \end{equation*}

Since a function $u \in \mathcal{V}$ is fully determined by its values $u_1, \ldots, u_{|V|}$ at the (finitely many) nodes, we may equivalently interpret $u$ as a column vector $(u_1, \ldots, u_{|V|})^T \in \mathbb{R}^{|V|}$ . We will freely make use of both the interpretation as function and as vector and do not distinguish between the two in our notation. As a consequence, we can represent linear operators that act on functions $u\;:\; V \to \mathbb{R}$ by $|V|$ -by- $|V|$ real matrices. Also in this context, we will not distinguish between the operator and matrix interpretations in our notation.

The standard (Euclidean) inner product $\langle \cdot, \cdot \rangle$ on $\mathbb{R}^{|V|}$ is defined as $\langle u, v \rangle \;:\!=\;u^T v$ . The norms $||\cdot ||_1$ and $||\cdot ||_2$ are the taxicab norm (i.e., $1$ -norm) and the Euclidean norm (i.e., $2$ -norm), respectively; that is, for vectors $w\in \mathbb{R}^n$ ,Footnote ³

\begin{equation*} \|w\|_1 \;:\!=\; \sum _{i=1}^n |w_n| \quad \text {and} \quad \|w\|_2 \;:\!=\; \left (\sum _{i=1}^n w_n^2\right )^{\frac {1}{2}}. \end{equation*}

If $C$ does not contain a row with only zeros, that is, if all row sums $(d_C)_i$ are positive, we define the $C$ -degree-weighted inner product to be

(2) \begin{equation} \langle u, v \rangle _C \;:\!=\; \sum _{i\in V} u_i v_i (d_C)_i. \end{equation}

For future reference we note that, if $C$ and $\tilde C$ are two adjacency matrices with positive row sums, then

(3) \begin{equation} \langle u, v \rangle _C = \langle D_{\tilde C}^{-1} D_C u, v \rangle _{\tilde C}. \end{equation}

Given an adjacency matrix $C$ , and a function $u\;:\; V \rightarrow \mathbb{R}$ , we define the graph total variation ( $TV_C$ ) and graph signless total variation ( $TV_C^+$ ) as:

\begin{equation*} TV_C (u)\;:\!=\; \frac {1}{2} \sum _{i,j \in V} c_{ij} |u_i -u_j| \quad \text {and} \quad TV_C^+ (u)\;:\!=\; \frac {1}{2} \sum _{i,j \in V} c_{ij} |u_i + u_j|. \end{equation*}

A vector-valued node function $u=(u^{(1)}, \ldots, u^{(K)})\;:\;V \rightarrow \mathbb{R}^{K}$ , where $u^{(l)}$ is the $l^{\text{th}}$ component of $u$ , can be interpreted as a matrix $U\in \mathbb{R}^{|V|\times K}$ , with elements $U_{il} \;:\!=\; u_i^{(l)}$ . We write $U_{*l}$ for the $l^{\text{th}}$ column of $U$ ; thus, the column vector $U_{*l}$ is the vector representation of the function $u^{(l)}\;:\; V \to \mathbb{R}$ . We write $U_{j*}$ for $U$ ’s $j^{\text{th}}$ row. As with the vector interpretation of real-valued node functions above, we freely use the interpretation as function or as matrix, as suits our purpose.

For vector-valued node functions $u$ , we generalise the definition of graph total variation and graph signless total variation to

\begin{align*} \mathcal{TV}_C (u)\;:\!=\;\sum _{l=1}^{K} TV_C (u^{(l)}) = \frac{1}{2} \sum _{l=1}^{K} \sum _{i,j \in V} c_{ij} \left |u_i^{(l)} - u_j^{(l)}\right | = \frac{1}{2} \sum _{l=1}^{K} \sum _{i,j \in V} c_{ij} \left |U_{il} - U_{jl}\right |,\\[5pt] \mathcal{TV}_C^+ (u)\;:\!=\;\sum _{l=1}^{K} TV_C^+ (u^{(l)}) = \frac{1}{2} \sum _{l=1}^{K} \sum _{i,j \in V} c_{ij} \left |u_i^{(l)} + u_j^{(l)}\right | = \frac{1}{2} \sum _{l=1}^{K} \sum _{i,j \in V} c_{ij} \left |U_{il} + U_{jl}\right |. \end{align*}

For a matrix $U\in \mathbb{R}^{n\times p}$ with entries $U_{ij}$ , its $C$ -Frobenius norm (with $C\in \mathbb{R}^{n\times n}$ symmetric positive definiteFootnote ⁴ ) is given by:

\begin{equation*} \|U\|_{\textrm {Fr},C} \;:\!=\; \sqrt {\textrm {tr}(U^TCU)} = \sqrt {\sum _{j=1}^p \sum _{k, l=1}^n C_{kl} U_{kj} U_{lj}}, \end{equation*}

where $\textrm{tr}$ denotes the trace of a square matrix. If $C$ is the identity matrix $I$ , this reduces to the standard Frobenius norm $\|U\|_{\textrm{Fr}} \;:\!=\; \|U\|_{_{\textrm{Fr},I}}$ . Because the standard Frobenius norm is submultiplicative, that is, for all matrices $U_1 \in \mathbb{R}^{n\times p}$ and $U_2 \in \mathbb{R}^{p\times q}$ , $\|U_1U_2\|_{\textrm{Fr}} \leq \|U_1\|_{\textrm{Fr}} \|U_2\|_{\textrm{Fr}}$ , the $C$ -Frobenius norm satisfies the following property:

(4) \begin{equation} \|U_1 U_2\|_{\textrm{Fr},C} = \|C^{\frac 12} U_1 U_2\|_{\textrm{Fr}} \leq \|C^{\frac 12} U_1\|_{\textrm{Fr}} \|U_2\|_{\textrm{Fr}} = \|U_1\|_{\textrm{Fr},C} \|U_2\|_{\textrm{Fr}}. \end{equation}

The infinity operator norm of the matrix $U$ is given by:

\begin{equation*} \|U\|_\infty \;:\!=\; \max _{i\in \{1, \ldots, n\}} \sum _{j=1}^p |U_{ij}|. \end{equation*}

Given $K\in \mathbb{N}$ , let $\mathcal{A}=\{A_l\}_{l=1}^{K}$ be a multisetFootnote ⁵ of pairwise disjoint subsets of $V$ that partition the node set $V$ , that is, $V=\bigcup _{l=1}^{K} A_l$ and $A_{l_1} \cap A_{l_2} = \emptyset$ if $l_1 \neq l_2$ . Note that $A_l$ could be empty, so the number of non-empty sets in $\mathcal{A}$ is at most $K$ . We call two partitions (with possibly different numbers of elements) equivalent, if every element of their symmetric difference is $\emptyset$ . The canonical representative of an equivalence class of partitions is the unique partition in the equivalence class that does not contain any copy of $\emptyset$ . Each canonical representative $\mathcal{A}$ is in bijective correspondence to a node assignment $c\;:\; V \to \{1, \ldots, |\mathcal{A}|\}$ : given a canonical representative $\mathcal{A}=\{A_l\}_{l=1}^{K}$ , define $c_j=l$ if and only if $j\in A_l$ ; conversely, given a node assignment $c\;:\; V \to \{1, \ldots, K\}$ , define $A_l=\{j \in V:c_j=l\}$ for $l\in \{1, \ldots, K\}$ .

We will refer to the sets $A_l$ in a partition as the communities defined by that partition. Also the terms ‘clusters’ or ‘classes’ may be used. We use an indicator function $\delta$ with $\delta (c_i,c_j)=1$ if nodes $i$ and $j$ are in the same community and $\delta (c_i,c_j)=0$ otherwise.

2.2. Laplacians for unsigned graphs

In this section, we define graph Laplacians for unsigned graphs determined by an adjacency matrix $C$ .

For a graph with weighted adjacency matrix $C$ , we define the graph Laplacian matrix $L_C$ as [Reference Chung16]:

\begin{align*} (L_C)_{ij} \;:\!=\; \begin{cases} (d_C)_i - c_{ii}, &\text{if } i=j, \\[5pt] -c_{ij}, &\text{otherwise}. \end{cases} \end{align*}

We include the dependence on $C$ explicitly in the notation, since we require graph Laplacians for various different graphs.

The Laplacian matrix can be written as:

\begin{equation*} L_C = D_C - C. \end{equation*}

If $D_C$ is invertible, the random walk graph Laplacian matrix $L_{C_{\textrm{rw}}}$ and the symmetrically normalised graph Laplacian matrix $L_{C_{\textrm{sym}}}$ are given by:

\begin{align*} L_{C_{\textrm{rw}}} &\;:\!=\; D_C^{-1} L_C = I - D_C^{-1} C,\\[5pt] L_{C_{\textrm{sym}}} &\;:\!=\;D_C^{-\frac{1}{2}} L_C D_C^{-\frac{1}{2}}= I-D_C^{-\frac{1}{2}} C D_C^{-\frac{1}{2}}, \end{align*}

respectively. It is well known that the random walk graph Laplacian and symmetrically normalised graph Laplacian have the same eigenvalues [Reference von Luxburg74].

In the special case that $C=W$ , we recall that $G=(V,E, \omega )$ is connected, and in particular there is no isolated node. Thus, for all $i\in V$ , $(d_W)_i\gt 0$ and hence the matrix $D_W$ is invertible.

Let $u\in \mathbb{R}^{|V|}$ . We compute, for $i\in V$ ,

(5) \begin{equation} (L_C u)_i = \sum _{j\in V} c_{ij} \left ( u_i - u_j \right ) \quad \text{and} \quad \left ( L_{C_{\textrm{rw}}} u\right )_i = \frac{1}{(d_C)_i} \sum _{j\in V} c_{ij} (u_i -u_j), \end{equation}

where we require $(d_C)_i\gt 0$ for the second computation. We note that any self-loops (i.e., $c_{ii}\gt 0$ ) do not contribute to the images of the Laplacian $L_C$ but will contribute to the degree normalisation in $L_{C_{\textrm{rw}}}$ .

Proof. Let $u,v\in \mathcal{V}$ .

It follows from the symmetry of $C$ and (5) that $L_C$ is self-adjoint with respect to the Euclidean inner product:

\begin{equation*} \langle u, L_C v\rangle = \sum _{i,j\in V} c_{ij} u_i (v_i-v_j) = \sum _{i,j\in V} \left (c_{ij} u_i v_i - c_{ji} u_j v_i \right ) = \sum _{i,j\in V} c_{ij} v_i (u_i-u_j) = \langle L_C u, v\rangle . \end{equation*}

Interchanging the indices $i$ and $j$ in this calculation shows in a straightforward way that $L_C$ is positive semidefinite:

(6) \begin{equation} \langle L_C u, u \rangle =\frac 12 \big (\langle L_C u, u \rangle + \langle u, L_C u\rangle \big ) =\frac{1}{2} \sum _{i,j\in V} c_{ij} (u_i-u_j)^2 \geq 0. \end{equation}

Similarly, the symmetrically normalised graph Laplacian $L_{C_{\textrm{sym}}}$ is self-adjoint with respect to the Euclidean inner product, since

\begin{equation*} \langle u, L_{C_{\textrm {sym}}} v \rangle = \langle D_C^{-\frac 12} u, L_C D_C^{-\frac 12} v \rangle = \langle L_C D_C^{-\frac 12} u, D_C^{-\frac 12} v \rangle = \langle L_{C_{\textrm {sym}}} u, v \rangle, \end{equation*}

and it is positive semidefinite with respect to the same inner product:

(7) \begin{align} \langle L_{C_{\textrm{sym}}} u, u \rangle = \frac 12 \big (\langle L_{W_{\textrm{sym}}} u, u \rangle + \langle u, L_{W_{\textrm{sym}}} u\rangle \big ) =\frac{1}{2} \sum _{i,j\in V} c_{ij} \left (\frac{u_i}{\sqrt{(d_C)_i}} - \frac{u_j}{\sqrt{(d_C)_j}} \right )^2 \geq 0. \end{align}

We recall that, for all $i\in V$ , $(d_C)_i\gt 0$ , because the diagonal matrix $D_C$ is invertible.

Finally, using (5) we compute

\begin{equation*} \langle u, L_{C_{\textrm {rw}}} v\rangle _C = \langle u, D_C L_{C_{\textrm {rw}}} v \rangle = \langle u, L_C v \rangle = \langle L_C u, v \rangle = \langle D_C L_{W_{\textrm {rw}}} u, v \rangle = \langle L_{C_{\textrm {rw}}} u, v \rangle _C \end{equation*}

and

(8) \begin{equation} \langle L_{C_{\textrm{rw}}} u, u \rangle _C = \frac{1}{2} \sum _{i,j\in V} c_{ij} \left (u_i -u_j \right )^2 \geq 0. \end{equation}

The signless graph Laplacian (matrix) $Q_C$ for a graph with weighted adjacency matrix $C$ , and its random walk and symmetrically normalised variants $Q_{C_{\textrm{rw}}}$ and $Q_{C_{\textrm{sym}}}$ , respectively, are defined as:

\begin{align*} Q_C &\;:\!=\;D_C + C,\\[5pt] Q_{C_{\textrm{rw}}} &\;:\!=\; D_C^{-1} Q_C = I + D_C^{-1} C, \\[5pt] Q_{C_{\textrm{sym}}}&\;:\!=\; D_C^{-\frac{1}{2}} Q_C D_C^{-\frac{1}{2}}= I + D_C^{-\frac{1}{2}} C D_C^{-\frac{1}{2}}. \end{align*}

Proof. The proofs are analogous to the proofs of the statements in Lemma 2.1. For future reference, we do note that, for all $u\in \mathcal{V}$ ,

(9) \begin{equation} (Q_C u)_i = \sum _{j\in V} c_{ij} \left ( u_i + u_j \right ) \quad \text{and} \quad \left ( Q_{C_{\textrm{rw}}} u\right )_i = \frac{1}{(d_C)_i} \sum _{j\in V} c_{ij} (u_i +u_j), \end{equation}

and

(10) \begin{align} \langle Q_C u, u \rangle &= \frac{1}{2} \sum _{i,j\in V} c_{ij} (u_i+u_j)^2 \geq 0, \end{align}

(11) \begin{align} \langle Q_{C_{\textrm{sym}}} u, u \rangle &= \frac{1}{2} \sum _{i,j\in V} c_{ij} \left (\frac{u_i}{\sqrt{(d_C)_i}} + \frac{u_j}{\sqrt{(d_C)_j}} \right )^2 \geq 0, \end{align}

(12) \begin{align} \langle Q_{C_{\textrm{rw}}} u,u \rangle _C &= \frac{1}{2} \sum _{i,j\in V} c_{ij} \left (u_i +u_j \right )^2 \geq 0. \end{align}

2.3. Review of the modularity function

‘Community structure’ is not a single well-defined concept. It attempts to capture the notion of partitioning a collection of individualsFootnote ⁷ (set of nodes) into meaningful clusters (classes and communities). What is meaningful depends on the context. Attempts to quantify ‘community structure’ tend to fall into one of two categories: comparisons with an externally available reference clustering, which for easy reference we will call the ‘ground truth’ partition,Footnote ⁸ and applications of mathematically defined ‘measures’Footnote ⁹ of community structure that do not require a ground truth. The approaches and measures from these categories are called extrinsic and intrinsic, respectively.

Extrinsic comparisons are used mostly when testing new methods on benchmark data sets for which the ‘preferred’ community structure is known through other means, such as the MNIST data set of handwritten digits (see Section 7.2). It may also be useful when investigating if the features that determine the network structure can be used to detect a community structure which is defined in terms of other features, for example, when a collaboration network of scientists is used in an attempt to construct a clustering that agrees with the areas of expertise of the scientists.

Intrinsic approaches have the advantages that no known ground truth is needed and that the mathematical formulation in terms of an optimisation problem for a well-defined measure of community structure can help in algorithm development. Once the mathematical problem has been formulated, the problem also becomes independent of context, although it still depends on the context how useful (in a practical, real-world sense) a given mathematical formulation is.

This work mainly falls in the second category, as we will use modularity [Reference Newman and Girvan57, Reference Newman56] as a measure of community structure and develop a new algorithm to approach the modularity optimisation problem. We will, however, also dip our toes into the first category, when we judge the outcomes of our algorithm not only by their modularity scores, but also through comparisons with ground truth community structures where these are available.

In the most general form that we encounter in this work, the definition of the modularity of a partition $\mathcal{A}$ of the node set of a graph $G=(V,E,\omega )$ is

(13) \begin{equation} \mathcal{Q}(\mathcal{A};\; W,P) \;:\!=\; \frac{1}{\textrm{vol}_W(V)} \sum _{i,j \in V} \left ( \omega _{ij} - p_{ij} \right ) \delta (c_i,c_j). \end{equation}

The matrix $P=(p_{ij})_{i,j=1}^{|V|}$ encodes the expected edge weight $p_{ij}$ between nodes $i$ and $j$ under a given null model, that is, a random graph satisfying some constraints. Since we are interested in undirected graphs with non-negative edge weights in this work, we assume $P$ to be symmetric and to have non-negative entries. We note that there is no restriction for the diagonal elements of $P$ to be zero. The null model can be thought of as describing graph structures that one would expect to see if there were no community structure present.

The modularity optimisation problem consists of finding $\textrm{argmax}_{\mathcal{A}} \mathcal{Q}(\mathcal{A};\; W,P)$ . For simplicity of notation, where we write $\textrm{argmax}_{\mathcal{A}}$ or $\max _{\mathcal{A}}$ , it is implicitly assumed that the maximum is taken over all partitions $\mathcal{A}$ of $V$ . We emphasise that $K$ , that is, the number of sets in $\mathcal{A}$ , is not fixed: finding an optimal $K$ is part of the optimisation problem. Because we do not assume that the sets in $\mathcal{A}$ are non-empty, solutions of the optimisation problem are necessarily non-unique, as $\mathcal{Q}(\mathcal{A};\; W,P) = \mathcal{Q}(\mathcal{A}\cup \{\emptyset \}; W,P)$ .Footnote ¹⁰

One of the commonly used null models is the Newman–Girvan (NG) model [Reference Newman and Girvan57], which is a random unweighted graph in which edges are assigned unbiasedly at random to pairs of distinct nodes, under the constraint that the resulting node degrees are equal to the degrees $(d_W)_i$ in the graph $G$ . The resulting expected edge weights are (approximatelyFootnote ¹¹ ) $p_{ij}^{NG} = \frac{(d_W)_i (d_W)_j}{\textrm{vol}_W(V)}$ , which are used as entries for the matrix $P^{\textrm{NG}}$ . Modularity with the NG null model thus can be written as:

(14) \begin{equation} \mathcal{Q}(\mathcal{A};\; W,P^{\textrm{NG}}) =\frac{1}{\textrm{vol}_W(V)} \sum _{i,j \in V} \left ( \omega _{ij} - \frac{(d_W)_i (d_W)_j}{\textrm{vol}_W(V)} \right ) \delta (c_i,c_j). \end{equation}

We note indeed that the degrees under the null model are the same as those in $G$ : $(d_{P^{\textrm{NG}}})_i = \sum _{j\in V} P^{\textrm{NG}}_{ij} = (d_W)_i$ . This has some very useful consequences, for example, $\textrm{vol}_{P^{\textrm{NG}}}(V) = \sum _{i \in V} (d_{P^{\textrm{NG}}})_i = \sum _{i\in V} (d_W)_i = \textrm{vol}_W(V)$ and both $L_{P^{\textrm{NG}}_{\textrm{sym}}}$ and $Q_{P^{\textrm{NG}}_{\textrm{sym}}}$ are self-adjoint with respect to the inner product from (2) with $C=W$ .

With this choice of null model $\mathcal{Q}(\mathcal{A};\; W,P^{\textrm{NG}})=0$ if $K=1$ , that is, if the partition is $\mathcal{A}=\{V\}$ .

According to Fortunato and Barthelemy [Reference Fortunato and Barthélemy22], there is a drawback in optimising equation (14) to find community partitions: it is difficult to find a community partition in networks that contain many small communities. It is argued that (in an unweighted graph) the number of communities $K$ that produces the maximum modularity score is (approximately) equal to $\sqrt{\frac{\textrm{vol}_W(V)}{2}}$ . A partition with many small communities tends to have more communities than this optimal value.

To solve the above problem, Arenas [Reference Arenas, Fernández and Gómez2] proposed a generalised modularity based on the Reichardt and Bornholdt method [Reference Reichardt and Bornholdt62]:

(15) \begin{align} \mathcal{Q}_\gamma (\mathcal{A};\; W, P) &\;:\!=\; \frac{1}{\textrm{vol}_W(V)} \sum _{i,j \in V} \left (\omega _{ij} - \gamma p_{ij} \right ) \delta (c_i,c_j), \end{align}

(16) \begin{align} \mathcal{Q}_\gamma (\mathcal{A};\; W, P^{\textrm{NG}}) &= \frac{1}{\textrm{vol}_W(V)} \sum _{i,j \in V} \left (\omega _{ij} - \gamma \frac{(d_W)_i (d_W)_j}{\textrm{vol}_W(V)} \right ) \delta (c_i,c_j), \end{align}

where $\gamma \gt 0$ is a resolution parameter [Reference Reichardt and Bornholdt62]. The distinction between (15) and (13) is the parameterFootnote ¹² $\gamma$ , which allows (15) to be more flexible and find more network community partitions; we note that $\mathcal{Q}_1=\mathcal{Q}$ . Nevertheless, there are still some issues with modularity optimisation even in this case.

Lancichinetti and Fortunato [Reference Lancichinetti and Fortunato41] make the case that for large enough $\gamma$ , the partition with maximum modularity score will split random subgraphs, which is unwanted behaviour, since random subgraphs should not be identified as having a community structure. On the other hand, they argue that for small enough $\gamma$ , there will be communities that contain multiple subgraphs, even when the number of internal edges in these subgraphs is large and there is only one inter-subgraph edge. Again, this is unwanted behaviour. Moreover, they show that it is difficult and in some cases even impossible to select a value for $\gamma$ that eliminates both these biases. In brief, for smaller values of $\gamma$ , we expect fewer clusters than for larger values.

A strategy for sampling the range of possible resolutions is presented in Jeub et al. [Reference Jeub, Sporns and Fortunato36]. Another approach to deal with this shortcoming is to investigate the stability of communities over multiple resolution scales, as in Mucha et al. [Reference Mucha, Richardson, Macon, Porter and Onnela54]. Since here we are primarily interested in the ability of our new algorithms to optimise modularity, rather than the appropriateness of the chosen resolution scale, we will not pursue those approaches here, and the problem of how to determine a good value for $\gamma$ in any given context is a matter outside the scope of this paper.

3.1. Reformulation of modularity optimisation for binary segmentation

In this subsection, we derive a new expression for the modularity $\mathcal{Q}(\mathcal{A})$ restricted to partitions $\mathcal{A}=\{A_l\}_{l=1}^K$ with $K=2$ , transforming the maximisation problem into a minimisation problem.

Let $u$ be a real-valued function on the node set $V$ , with value $u_i$ on node $i$ . We define the set of $\{ -1,+1 \}$ -valued node functions as:

\begin{equation*} \mathcal {V}_{\textrm {bin}} \;:\!=\; \big \{u\;:\; V\to \{-1, +1\}\big \}. \end{equation*}

Specially, if $u \in \mathcal{V}_{\textrm{bin}}$ , we define the sets

\begin{equation*} V_1 =\{i \in V\;:\; u_i =1\} \quad \text {and} \quad V_{-1} =\{i \in V\;:\; u_i =-1\}. \end{equation*}

We consider a partition $\mathcal{A} = \{A_1, A_2\}$ , with $A_1=V_1$ , $A_2=V_{-1}$ , and corresponding node assignment $c$ .

If $u_i \in \mathcal{V}_{\textrm{bin}}$ , that is, $u_i \in \{ -1,+1\}$ , then $u_i^2 + u_j^2 =2$ , which implies that

(17) \begin{equation} (u_i u_j +1) = -\frac{1}{2} (u_i -u_j)^2 + 2 \quad \text{and} \quad -(u_i u_j +1) = -\frac{1}{2}(u_i +u_j)^2. \end{equation}

If $i,j \in V_1$ or $i,j \in V_{-1}$ , then $\delta (c_i,c_j)=1=\frac{1}{2}(u_i u_j +1)$ . Similarly, if $u_i \neq u_j$ , then $\delta (c_i,c_j)=0=$ $\frac{1}{2}(u_i u_j +1)$ . Hence, the modularity $\mathcal{Q}$ for $K=2$ clusters can be rewritten as:Footnote ¹³

(18) \begin{align} \mathcal{Q}_\gamma (\mathcal{A};\; W,P) &=\frac{1}{2\textrm{vol}_W(V)} \sum _{i,j \in V} \left ( \omega _{ij} - \gamma p_{ij} \right ) (u_i u_j +1) \nonumber \\[5pt] &=\frac{1}{2\textrm{vol}_W(V)} \sum _{i,j \in V} \omega _{ij}(u_i u_j +1) - \frac{\gamma }{2\textrm{vol}_W(V)} \sum _{i,j \in V} p_{ij} (u_i u_j +1) \nonumber \\[5pt] &= -\frac{1}{4\textrm{vol}_W(V)} \sum _{i,j \in V} \omega _{ij}(u_i - u_j)^2 + \frac{1}{\textrm{vol}_W(V)} \sum _{i,j \in V} \omega _{ij} -\frac{\gamma }{4\textrm{vol}_W(V)} \sum _{i,j \in V} p_{ij} (u_i + u_j)^2 \nonumber \\[5pt] &=-\frac{1}{\textrm{vol}_W(V)} \left [\frac 14\sum _{i,j \in V} \omega _{ij}(u_i - u_j)^2 + \frac \gamma 4 \sum _{i,j \in V} p_{ij} (u_i + u_j)^2 \right ] + \frac{1}{\textrm{vol}_W(V)} \sum _{i,j \in V} \omega _{ij}. \end{align}

For all $u \in \mathcal{V}_{\textrm{bin}}$ , one obtains

\begin{equation*} (u_i - u_j)^2 = \begin {cases} 0, &\text {if } u_i=u_j,\\[5pt] 4, &\text {if } u_i\neq u_j, \end {cases} \quad \text {and} \quad 2|u_i -u_j| = \begin {cases} 0, &\text {if } u_i=u_j,\\[5pt] 4, &\text {if } u_i\neq u_j. \end {cases} \end{equation*}

Therefore, if $u \in \mathcal{V}_{\textrm{bin}}$ , then $(u_i - u_j)^2 = 2|u_i -u_j|$ . Similarly, $(u_i + u_j)^2 = 2|u_i + u_j|$ , if $u \in \mathcal{V}_{\textrm{bin}}$ .

Since $\textrm{vol}_W(V)= \sum _{i,j \in V} \omega _{ij}$ , the third term in (18) equals one; in particular, it does not depend on $u$ . Thus,

\begin{equation*} \mathcal {Q}_\gamma (\mathcal {A};\; W, P) = -\frac {1}{\textrm {vol}_W(V)}\mathcal {Q}_{\textrm {bin},\gamma }(u;\; W, P) + 1, \end{equation*}

where

(19) \begin{align} \mathcal{Q}_{\textrm{bin},\gamma }(u;\; W, P) &\;:\!=\;\frac{1}{4} \sum _{i,j \in V} \omega _{ij}(u_i - u_j)^2 + \frac{\gamma }{4} \sum _{i,j \in V} p_{ij} (u_i + u_j)^2\nonumber \\[5pt] &= \frac 12 \sum _{i,j \in V} \omega _{ij}|u_i - u_j| + \frac \gamma 2 \sum _{i,j \in V} p_{ij} |u_i + u_j|\nonumber \\[5pt] &=TV_W(u) + \gamma TV_P^+(u). \end{align}

The maximisation of modularity $\mathcal{Q}_\gamma (\mathcal{A};\; W,P)$ from (15) over partitions $\mathcal{A}$ with $K=2$ is equivalent to the minimisation of $\mathcal{Q}_{\textrm{bin},\gamma } (u;\; W, P)$ (19) over all functions $u \in \mathcal{V}_{\textrm{bin}}$ (with the correspondence between $\mathcal{A}$ and $u$ introduced above).

Similarly to Newman [Reference Newman56], we define the modularity matrix as:

\begin{equation*} B_\gamma \;:\!=\; W-\gamma P \end{equation*}

and denote its entries by $b_{ij}$ . Since $W$ and $P$ are assumed to be symmetric, so is $B_\gamma$ . We note that $B_\gamma$ does not need to have zeros on its diagonal.

If the condition $D_W=D_P$ is satisfied, as is the case, we recall from Section 2.3, if $P=P^{\textrm{NG}}$ is given by the NG null model, then

(20) \begin{equation} D_{B_\gamma } = (1-\gamma ) D_W. \end{equation}

In our considerations above, we have split the matrix $B_\gamma$ into the matrix $W$ with non-negative entries and the matrix $-\gamma P$ with non-positive entries, but we can also write $B_\gamma = B^+_\gamma - B^-_\gamma$ , where $B^+_\gamma$ and $B^-_\gamma$ are $|V|$ -by- $|V|$ matrices with entries

\begin{equation*} (b^+_\gamma )_{ij}\;:\!=\;\max \{ (b_\gamma )_{ij},0 \} \quad \text {and} \quad (b^-_\gamma )_{ij}\;:\!=\;-\min \{ (b_\gamma )_{ij},0 \}, \end{equation*}

respectively, with $(b_\gamma )_{ij} \;:\!=\; \omega _{ij}-\gamma p_{ij}$ the entries of $B_\gamma$ . Per definition $B^+_\gamma$ has non-negative entries and $-B^-_\gamma$ non-positive entries, yet in general $B^+_\gamma \neq W$ and $B^-_\gamma \neq \gamma P$ . Thus, this split will give another way to rewrite $\mathcal{Q}_\gamma$ analogously to (18):

\begin{align*} \mathcal{Q}_\gamma (\mathcal{A};\; W, P) &= -\frac{1}{2\textrm{vol}_W(V)} \left [ \frac{1}{2} \sum _{i,j \in V} (b^+_\gamma )_{ij}(u_i - u_j)^2 + \frac{1}{2} \sum _{i,j \in V} (b^-_\gamma )_{ij} (u_i + u_j)^2 \right ] + \frac{1}{\textrm{vol}_W(V)} \sum _{i,j \in V} (b^+_\gamma )_{ij}\\[5pt] &= -\frac{1}{\textrm{vol}_W(V)} \mathcal{Q}_{\textrm{bin},1}(u;\; B^+_\gamma, B^-_\gamma ) + \frac{1}{\textrm{vol}_W(V)} \sum _{i,j \in V} (b^+_\gamma )_{ij}. \end{align*}

Because the final term on the right-hand side does not depend on $u$ , we see that maximising $\mathcal{Q}_\gamma (\mathcal{A};\; W, P)$ over all bipartitions $\mathcal{A}$ is equivalent to minimising $\mathcal{Q}_{\textrm{bin},1}(u;\; B^+_\gamma, B^-_\gamma )$ over all $u\in \mathcal{V}_{\textrm{bin}}$ .

In the remainder of this paper, we would like to be able to consider graph Laplacians based on $B^+_\gamma$ and $B^-_\gamma$ . To be able to define random walk and symmetrically normalised (signless) graph Laplacians based on these matrices, we require the degree matrices $D_{B^+_\gamma }$ and $D_{B^-_\gamma }$ to be invertible. Additionally, for the symmetrically normalised (signless) graph Laplacians, we also require these matrices to be positive semidefiniteFootnote ¹⁴ so that their square roots are uniquely defined. This property is easily checked to hold, since both matrices are diagonal with non-negative entries.

The following lemma and corollary collect some useful results about their invertibility for easy reference.

Lemma 3.1. Let $i\in V$ and $\gamma \in (0,\infty )$ . Then $(d_{B^+_\gamma })_i \neq 0$ if and only if there exists a $j\in \mathcal{N}_W(i)$ such that

(21) \begin{equation} \gamma \lt p_{ij}^{-1} \omega _{ij}, \end{equation}

where we define $p_{ij}^{-1}\omega _{ij} \;:\!=\; +\infty$ if $p_{ij}=0$ .

Similarly, $(d_{B^-_\gamma })_i \neq 0$ if and only if there exists a $j\in \mathcal{N}_P(i)$ such that

(22) \begin{equation} \gamma \geq p_{ij}^{-1} \omega _{ij}. \end{equation}

Consequently, $D_{B^+_\gamma }$ is invertible if and only if

\begin{equation*} \gamma \lt \min _{i\in V} \max _{j\in \mathcal {N}_W(i)} p_{ij}^{-1} \omega _{ij} \end{equation*}

and $D_{B^-_\gamma }$ is invertible if and only if

\begin{equation*} \gamma \gt \max _{i\in V} \min _{j\in \mathcal {N}_P(i)} p_{ij}^{-1} \omega _{ij}. \end{equation*}

3.2. Generalisation to multiple clusters

Next, we extend the approach of Section 3.1 to identify appropriate partitions of the node set into multiple clusters. In this subsection, we consider $K$ to be fixed (but not necessarily equal to two). Then a partition $\mathcal{A}=\{A_l\}_{l=1}^K$ with $K$ parts is completely described by a function $u=(u^{(1)}, \ldots, u^{(K)})\;:\;V \rightarrow \mathbb{R}^{K}$ , where $u_i^{(l)}=1$ if and only if $i \in A_l$ and $u_i^{(l)}=-1$ otherwise. We can also encode this information in a matrix $U\in \mathbb{R}^{|V|\times K}$ , with elements $U_{il} \;:\!=\; u_i^{(l)}$ . We denote the set of all matrices corresponding to a partition in $K$ parts by:

(23) \begin{equation} Pt(K) \;:\!=\; \big \{ U \in \mathbb{R}^{|V| \times K} \;:\; \forall i\in V \, \, \forall l\in \{1, \ldots, K\} \, \, U_{il} \in \{-1, 1\} \text{ and } \sum _{k=1}^K U_{ik} = 2-K\big \}. \end{equation}

We may call such matrices partition matrices. The last condition in the definition of $Pt(K)$ guarantees that each node belongs to exactly one cluster. We recall that $\mathcal{A}$ may contain empty sets, thus, for all $l\in \{1, \ldots, K\}$ , $\sum _{j\in V} u_j^{(l)} \in [\!-\!|V|, |V|] \cap \mathbb{Z}$ .

Because we will require the use of (signless) total variation with respect to various different matrices, for the moment we consider an arbitrary real-valued matrix $C\in \mathbb{R}^{|V|\times |V|}$ with entries $c_{ij}$ .

We briefly compare the current set-up with that from Section 3.1 in the case $K=2$ with partition $\mathcal{A}=\{A_1, A_2\}$ . In the notation of the current section, this partition is encoded by a matrix $U\in Pt(2)$ ; in the set-up of Section 3.1 we encode the same partition by a function $v\in \mathcal{V}_{\textrm{bin}}$ . These two encodings are related via $U_{*1}=v$ and $U_{*2}=-v$ and thus

(24) \begin{align} \mathcal{TV}_C(U) &= \frac 12 \left (\sum _{i,j\in V} c_{ij} |U_{i1}-U_{j1}| + \sum _{i,j\in V} c_{ij} |U_{i2}-U_{j2}|\right )\nonumber \\[5pt] &= \frac 12 \left (\sum _{i,j\in V} c_{ij} |v_i-v_j| + \sum _{i,j\in V} c_{ij} |-v_i+v_j|\right ) = 2 TV_C(v) \end{align}

and similarly $\mathcal{TV}_C^+(U) = 2 TV_C^+(v)$ .

We return to the case of general $K$ . For notational convenience, in the following computations we write $\tilde U_{ijl} \;:\!=\; c_{ij}|U_{il}-U_{jl}|$ and $\hat U_{ijl} \;:\!=\; c_{ij}|U_{il}+U_{jl}|$ . For a matrix $U \in Pt(K)$ , we obtain that

(25) \begin{align} |U_{il} - U_{jl}| = \begin{cases} 0 &\text{if } i,j \in A_l \text{ or } i,j \in A_l^c,\\[5pt] 2 &\text{if } i \in A_l,j \in A_l^c \text{ or } i \in A_l^c,j \in A_l,\\[5pt] \end{cases} \end{align}

(26) \begin{align} |U_{il} + U_{jl}| = \begin{cases} 0 &\text{if } i\in A_l,j\in A_l^c \text{ or } i \in A_l^c,j \in A_l,\\[5pt] 2 &\text{if } i,j \in A_l \text{ or } i,j \in A_l^c. \end{cases} \end{align}

(We write $A_l^c$ for the complement $V\setminus A_l^c$ .) Thus, the generalised graph total variation and graph signless total variation of $U \in Pt(K)$ on a graph with adjacency matrix $C$ can be represented as:

\begin{align*} \mathcal{TV}_C(U) &= \frac{1}{2} \sum _{l=1}^K \left ( \sum _{i,j\in A_l} \tilde U_{ijl} + \sum _{i\in A_l, j \in A_l^c} \tilde U_{ijl} + \sum _{i\in A_l^c,j \in A_l} \tilde U_{ijl} + \sum _{i,j\in A_l^c} \tilde U_{ijl} \right )\\[5pt] &= \sum _{l=1}^K \left (\sum _{i\in A_l, j \in A_l^c} c_{ij} + \sum _{i\in A_l^c, j \in A_l} c_{ij}\right ) \end{align*}

and

\begin{align*} \mathcal{TV}^+_C(U) &= \frac{1}{2} \sum _{l=1}^K \left ( \sum _{i,j\in A_l} \hat U_{ijl} + \sum _{i\in A_l, j \in A_l^c} \hat U_{ijl} + \sum _{i\in A_l^c,j \in A_l} \hat U_{ijl} + \sum _{i,j\in A_l^c} \hat U_{ijl}\right )\\[5pt] &=\frac{1}{2} \sum _{l=1}^K \left ( \sum _{i,j\in A_l} \hat U_{ijl} + \sum _{i,j\in A_l^c} \hat U_{ijl} \right )\\[5pt] &=\sum _{l=1}^K \sum _{i,j\in A_l} c_{ij} + \sum _{l=1}^K \left (\sum _{i,j\in V} c_{ij} + \sum _{i,j\in A_l} c_{ij} - \sum _{i\in V, j\in A_l} c_{ij} - \sum _{i\in A_l, j\in V} c_{ij}\right )\\[5pt] &=\sum _{l=1}^K \sum _{i,j\in A_l} c_{ij} + \left ( \sum _{l=1}^K \sum _{i,j\in A_l} c_{ij} +(K-2)\textrm{vol}_C(V) \right )\\[5pt] &= 2\sum _{l=1}^K \sum _{i,j\in A_l} c_{ij} + (K-2)\textrm{vol}_C(V). \end{align*}

Using these expressions for the generalised total variation and generalised signless total variation, if $\mathcal{A}$ contains $K$ subsets, the modularity from (13) can be written as:Footnote ¹⁵

\begin{align*} \mathcal{Q}_\gamma (\mathcal{A};\; W, P) &=\frac{1}{\textrm{vol}_W(V)} \sum _{l=1}^K \sum _{i,j \in A_l} \omega _{ij} -\frac{\gamma }{\textrm{vol}_W(V)} \sum _{l=1}^K \sum _{i,j \in A_l} p_{ij} \nonumber \\[5pt] &=\frac{1}{\textrm{vol}_W(V)} \sum _{l=1}^K \left ( \sum _{i \in A_l, j\in V} \omega _{ij} - \sum _{i \in A_l, j\in A_l^c} \omega _{ij} \right ) - \frac{\gamma }{\textrm{vol}_W(V)} \sum _{l=1}^K \sum _{i,j \in A_l} p_{ij} \end{align*}

(27) \begin{align} &=\frac{1}{\textrm{vol}_W(V)} \sum _{i,j=1}^{|V|} \omega _{ij} - \frac{1}{\textrm{vol}_W(V)} \sum _{l=1}^K \sum _{i \in A_l, j\in A_l^C} \omega _{ij} - \frac{\gamma }{\textrm{vol}_W(V)} \sum _{l=1}^K \sum _{i,j \in A_l} p_{ij} \nonumber \\[5pt] &= 1 - \frac{1}{\textrm{vol}_W(V)} \left ( \frac 12 \sum _{l=1}^K \sum _{i \in A_l, j\in A_l^c} \omega _{ij} + \frac 12 \sum _{l=1}^K \sum _{i \in A_l^c, j\in A_l} \omega _{ij} + \gamma \sum _{l=1}^K \sum _{i,j \in A_l} p_{ij} \right )\nonumber \\[5pt] &= 1 - \frac{1}{\textrm{vol}_W(V)} \left ( \frac 12 \mathcal{TV}_W (U) + \frac \gamma 2 \mathcal{TV}_P^+(U) - \frac{\gamma }{2}(K-2)\textrm{vol}_P(V)\right ), \end{align}

where in the last line $U\in Pt(K)$ corresponds to the partition $\mathcal{A}$ .

Thus,

\begin{equation*} \mathcal {Q}_\gamma (\mathcal {A};\; W, P) = - \frac {1}{\textrm {vol}_W(V)} \mathcal {Q}_{\textrm {mul},\gamma }(U;\; W, P) + 1 + \frac {\gamma (K-2) \textrm {vol}_P(V)}{2 \textrm {vol}_W(V)}, \end{equation*}

if $U$ corresponds to $\mathcal{A}$ and

\begin{equation*} \mathcal {Q}_{\textrm {mul},\gamma }(U;\; W, P) \;:\!=\; \frac {1}{2} \mathcal {TV}_W (U) + \frac {\gamma }{2} \mathcal {TV}_P^+ (U). \end{equation*}

Thus, the maximisation of $\mathcal{Q}_\gamma (\mathcal{A};\; W, P)$ over all partitions $\mathcal{A}$ with fixed $K$ is equivalent to the minimisation of $\mathcal{Q}_{\textrm{mul},\gamma }(U;\; W, P)$ over all $U\in Pt(K)$ . If $K=2$ and $U$ and $v$ are related as in (24), then $\mathcal{Q}_{\textrm{mul},\gamma }(U;\; W, P) = \mathcal{Q}_{\textrm{bin},\gamma }(v;\; W, P)$ .

Similar to what we did in the case of bipartitions in Section 3.1, if we split $B_\gamma =W-\gamma P$ as $B_\gamma =B^+_\gamma -B^-_\gamma$ , we can do the analogue computation to (27) to find

\begin{align*} \mathcal{Q}_\gamma (\mathcal{A};\; W, P) &= \frac{1}{\textrm{vol}_W(V)} \sum _{i,j \in V} \left ( (b^+_\gamma )_{ij} - (b^-_\gamma )_{ij} \right ) \delta (c_i,c_j)\\[5pt] &= - \frac{1}{\textrm{vol}_W(V)} \mathcal{Q}_{mul,1}(U;\; B^+_\gamma, B^-_\gamma ) + \frac{2 \textrm{vol}_{B^+_\gamma }(V) + (K-2)\textrm{vol}_{B^-_\gamma }(V)}{2\textrm{vol}_W(V)}. \end{align*}

Hence, we can also maximise $\mathcal{Q}_\gamma (\mathcal{A};\; W, P)$ over all partitions $\mathcal{A}$ containing $K$ subsets by minimising $\mathcal{Q}_{mul,1}(U;\; B^+_\gamma, B^-_\gamma )$ over all $U\in Pt(K)$ . We emphasise that in $\mathcal{Q}_{mul,1}$ we choose $\gamma =1$ , since the influence of $\gamma$ is now in the matrices $B^+_\gamma$ and $B^-_\gamma$ , rather than in the structure of the function(al).

4.1. Binary classification with graph Ginzburg–Landau functionals

A central object in the diffuse-interface approach of [Reference Bertozzi and Flenner4] is the graph GL functional $f_\varepsilon \;:\; \mathcal{V} \to \mathbb{R}$ defined by:

(28) \begin{align} f_\varepsilon (u):&= \frac 14 \sum _{i,j\in V} \omega _{ij} \left ( u_i - u_j \right )^2 + \frac{1}{\varepsilon } \sum _{i\in V} \Phi \left ( u_i \right ) = \frac 12 \langle u, L_W u \rangle + \frac{1}{\varepsilon } \sum _{i\in V} \Phi \left ( u_i \right ) & \text{(by (6))}\nonumber \\[5pt] &= \frac 12 \langle u, L_{W_{\textrm{rw}}} u \rangle _W + \frac{1}{\varepsilon } \sum _{i\in V} \Phi \left ( u_i \right ) & \text{(by (8))}, \end{align}

where $\varepsilon \gt 0$ is a parameter and the function $\Phi (u)$ is a double-well potential with two minima. For example, a classical choice is the polynomial $\Phi (u)= \frac{1}{4}(u^2 -1)^2$ that has minima at $u=-1$ and $u=+1$ . The first of the two terms in (28) is called the graph Dirichlet energy.

In Van Gennip and Bertozzi [Reference van Gennip and Bertozzi72], it is proven that if $\varepsilon \downarrow 0$ , then the sequence of functionals $f_\varepsilon$ $\Gamma$ -converges to

\begin{equation*} f_0(u) \;:\!=\; \begin {cases} TV_W(u), &\text {if } u\in \mathcal {V}_{\textrm {bin}},\\[5pt] +\infty, &\text {otherwise}. \end {cases} \end{equation*}

For details about $\Gamma$ -convergence, we refer to Dal Maso [Reference Maso47] and Braides [Reference Braides8]. Here, it suffices to note that $\Gamma$ -convergence of $f_\varepsilon$ combined with an equicoercivity condition, which is also satisfied in this case (see [Reference van Gennip and Bertozzi72]), implies that minimisers of $f_\varepsilon$ converge to minimisers of $f_0$ as $\varepsilon \downarrow 0$ .

Similarly, in Keetch and Van Gennip [Reference Keetch and van Gennip39], the signless graph GL functional $f_0^+\;:\; \mathcal{V} \to \mathbb{R}$ , defined as

\begin{align*} f_\varepsilon ^+ (u):&= \frac 14 \sum _{i,j\in V} \omega _{ij} \left ( u_i + u_j \right )^2 + \frac{1}{\varepsilon } \sum _{i\in V} \Phi \left ( u_i \right ) = \frac 12 \langle u, Q_W u \rangle + \frac{1}{\varepsilon } \sum _{i\in V} \Phi \left ( u_i \right ) & \text{(by (10))} \\[5pt] &= \frac 12 \langle u, Q_{W_{\textrm{rw}}} u \rangle _W + \frac{1}{\varepsilon } \sum _{i\in V} \Phi \left ( u_i \right ) & \text{(by (12))}, \end{align*}

is introduced and it is proven that $f_\varepsilon ^+$ $\Gamma$ -converges to

\begin{equation*} f_0^+(u) \;:\!=\; \begin {cases} TV_W^+(u), &\text {if } u\in \mathcal {V}_{\textrm {bin}},\\[5pt] +\infty, &\text {otherwise}, \end {cases} \end{equation*}

as $\varepsilon \downarrow 0$ . Also in this case the equicoercivity condition that is required to conclude convergence of minimisers of $f_\varepsilon ^+$ to minimisers of $f_0^+$ is satisfied (see [Reference Keetch and van Gennip39]). We call the first of the two terms in $f_0^+$ the signless graph Dirichlet energy.

A straightforward adaptation of the proofs in [Reference van Gennip and Bertozzi72, Theorems 3.1 and 3.2] and [Reference Keetch and van Gennip39, Lemmas 4.3 and 4.4] shows that if we define

\begin{equation*} f_{\varepsilon, \gamma }^{\pm }(u;\; W, P) \;:\!=\; \frac 14 \sum _{i,j\in V} \omega _{ij} \left ( u_i - u_j \right )^2 + \frac \gamma 4 \sum _{i,j\in V} p_{ij} \left ( u_i + u_j \right )^2 + \frac {1}{\varepsilon } \sum _{i\in V} \Phi \left ( u_i \right ), \end{equation*}

then $f_{\varepsilon, \gamma }^\pm (\cdot ;\; W, P)$ $\Gamma$ -converges to

\begin{equation*} f_{0,\gamma }^\pm (u;\; W, P) \;:\!=\; \begin {cases} \mathcal {Q}_{\textrm {bin},\gamma }(u;\; W, P), &\text {if } u\in \mathcal {V}_{\textrm {bin}},\\[5pt] +\infty, &\text {otherwise}, \end {cases} \end{equation*}

and the required equicoercivity conditions are again satisfied that allow us to conclude that minimisers of $f_{\varepsilon, \gamma }^\pm (\cdot ;\; W, P)$ converge to minimisers of $\mathcal{Q}_{\textrm{bin},\gamma }(\cdot ;\; W, P)$ from (19). We provide more details about the proof when we consider the case with multiple clusters in Theorem 4.1.

Similarly, $f_{\varepsilon, \gamma }^\pm (\cdot ;\; B^+_\gamma, B^-_\gamma )$ $\Gamma$ -converges to $f_{0,\gamma }^\pm (u;\; B^+_\gamma, B^-_\gamma )$ and the required equicoercivity conditions are again satisfied.

For small $\varepsilon$ , minimisers of $f_{\varepsilon, \gamma }^\pm (\cdot ;\; W, P)$ (or $f_{\varepsilon, 1}^\pm (\cdot ;\; B^+_\gamma, B^-_\gamma )$ ) thus approximate (in the sense just describedFootnote ¹⁷ ) minimisers of $\mathcal{Q}_{\textrm{bin},\gamma }(\cdot ;\; W, P)$ (or $\mathcal{Q}_{\textrm{bin},1}(\cdot ;\; B^+_\gamma, B^-_\gamma )$ ), which we know to be equivalent to maximisers of modularity $\mathcal{Q}(\cdot ;\; W, P)$ restricted to bipartitions of $V$ .

Finding global minimisers of the non-convex function $f_\varepsilon ^\pm$ is a difficult task. Instead, we can focus on finding local minimisers and hope that these, in practice, are also good approximations for maximisers of $\mathcal{Q}_\gamma$ . To which extent this hope proves to be justified will be investigated numerically in Section 7. For now, we focus on the problem of finding such local minimisers.

One possible method is to compute a gradient flow of $f_\varepsilon ^\pm$ (see, e.g., Van Gennip et al. [Reference van Gennip, Guillen, Osting and Bertozzi73] and Keetch and Van Gennip [Reference Keetch and van Gennip39]). For $u,v \in \mathcal{V}$ and $s\in \mathbb{R}$ , using the self-adjointness of $L_W$ , $L_{W_{\textrm{rw}}}$ , $Q_W$ , and $Q_{W_{\textrm{rw}}}$ , we computeFootnote ¹⁸

\begin{align*} \left . \frac{d}{ds} f_{\varepsilon, \gamma }^\pm (u+sv;\; W, P) \right |_{s=0} &= \langle L_W u, v\rangle + \gamma \langle Q_P u, v\rangle + \frac 1\varepsilon \langle \Phi '\circ u, v\rangle \\[5pt] &= \langle L_{W_{\textrm{rw}}} u, v\rangle _W + \gamma \langle Q_{P_{\textrm{rw}}} u, v\rangle _W + \frac 1\varepsilon \langle D_W^{-1}\Phi ' \circ u, v\rangle _W. \end{align*}

Thus, the gradient flows of $f_{\varepsilon, \gamma }^\pm (\cdot ;\; W, P)$ with respect to the Euclidean inner product and the degree-weighted inner product are

(29) \begin{equation} \frac{du}{dt} = - L_W u - \gamma Q_P u - \frac 1\varepsilon \Phi ' \circ u \quad \text{and} \quad \frac{du}{dt} = -L_{W_{\textrm{rw}}} u - \gamma Q_{P_{\textrm{rw}}} u - \frac 1\varepsilon D_W^{-1} \Phi ' \circ u, \end{equation}

respectively. As is standard for gradient flows, a dependence on ‘time’ $t$ has been introduced so that we now may interpret $u$ as a function $u\;:\; \mathbb{R} \to \mathcal{V}$ .

In a completely analogous manner, we determine

\begin{equation*} \left . \frac {d}{ds} f_{\varepsilon, 1}^\pm (u+sv;\; B^+_\gamma, B^-_\gamma ) \right |_{s=0} = \langle L_{B^+} u, v\rangle + \langle Q_{B^-_\gamma } u, v\rangle + \frac 1\varepsilon \langle \Phi '\circ u, v\rangle \end{equation*}

and thus the gradient flow of $f_{\varepsilon, 1}^\pm (\cdot ;\; B^+_\gamma, B^-_\gamma )$ with respect to the Euclidean inner product is

(30) \begin{equation} \frac{du}{dt} = - L_{B^+} u - Q_{B^-} u - \frac 1\varepsilon \Phi ' \circ u. \end{equation}

For the gradient flow with respect to a degree-weighted inner product, the situation is more complicated. We recall the definition of the degree-weighted inner product in (2) and explicitly will be using the following versions:

(31) \begin{equation} \langle u, v \rangle _{B_\gamma } = \sum _{i\in V} u_i v_i (d_{B_\gamma })_i, \quad \langle u, v \rangle _{B^+_\gamma } = \sum _{i\in V} u_i v_i (d_{B^+_\gamma })_i, \quad \text{and} \quad \langle u, v \rangle _{B^-_\gamma } = \sum _{i\in V} u_i v_i (d_{B^-_\gamma })_i. \end{equation}

By (20), we know that, if $D_W=D_P$ and $\gamma =1$ , then $D_{B_1} = 0$ and therefore $D_{B^+_1}=D_{B^-_1}$ . Thus in this case, the $B^+_\gamma$ - and $B^-_\gamma$ -degree-weighted inner products are equal and the $B_\gamma$ -degree-weighted inner product is always zero.Footnote ¹⁹ If additionally $D_{B^+_1}=D_{B^-_1}$ is invertible (e.g., because the remaining assumption from Corollary 3.2 part 2 is satisfied), then we can write

\begin{equation*} \langle L_{B^+_1} u, v\rangle = \langle L_{{B^+_1}_{\textrm {rw}}} u, v\rangle _{B^+_1}, \quad \langle Q_{B^-_1} u, v\rangle = \langle Q_{{B^-_1}_{\textrm {rw}}} u, v\rangle _{B^+_1}, \quad \text {and} \quad \langle \Phi '(u), v\rangle = \langle D_{B^+_1}^{-1} \Phi '(u), v\rangle _{B^+_1}. \end{equation*}

Hence, if $\gamma =1$ , the gradient flow with respect to the $B^+$ -degree-weighted inner product is

(32) \begin{equation} \frac{du}{dt} = -L_{{B^+_1}_{\textrm{rw}}} u - Q_{{B^-_1}_{\textrm{rw}}} u - \frac 1\varepsilon D_{B^+_1}^{-1} \Phi ' \circ u. \end{equation}

If $\gamma \neq 1$ and $D_W=D_P$ , then the degree matrices with respect to $B^+_\gamma$ and $B^-_\gamma$ are no longer the same, but rather we have

(33) \begin{equation} D_{B^-_\gamma } = D_{B^+_\gamma } - D_{B_\gamma } = D_{B^+_\gamma } - (1-\gamma ) D_W. \end{equation}

This means we have to make a choice to use either the $B^+_\gamma$ - or $B^-_\gamma$ -degree-weighted inner product, as they are no longer identical. That choice will influence the resulting equations. We choose the former; calculations for the alternative choice are similar.

With this choice, and still assuming that $D_{B^-_\gamma }$ and $D_{B^+_\gamma }$ are invertible, for example because the assumptions from Corollary 3.2 part 1 hold, we still have

\begin{equation*} \langle L_{B^+_\gamma } u, v\rangle = \langle L_{{B^+_\gamma }_{\textrm {rw}}} u, v\rangle _{B^+_\gamma } \quad \text {and} \quad \langle \Phi '(u), v\rangle = \langle D_{B^+_\gamma }^{-1} \Phi '(u), v\rangle _{B^+_\gamma }, \end{equation*}

yet when rewriting the remaining term in $\left . \frac{d}{ds} f_{\varepsilon, 1}^\pm (u+sv;\; B^+_\gamma, B^-_\gamma ) \right |_{s=0}$ , the difference in $B^+_\gamma$ - and $B^-_\gamma$ -degree-weighted inner products is important:

\begin{align*} \langle Q_{B^-_\gamma } u, v\rangle &= \langle Q_{{B^-_\gamma }_{\textrm{rw}}} u, v\rangle _{B^-_\gamma } = \langle Q_{{B^-_\gamma }_{\textrm{rw}}} u, v\rangle _{B^+_\gamma } - \langle Q_{{B^-_\gamma }_{\textrm{rw}}} u, v\rangle _{B_\gamma } \\[5pt] &= \langle Q_{{B^-_\gamma }_{\textrm{rw}}} u, v\rangle _{B^+_\gamma } - \langle D_{B^+_\gamma }^{-1} D_{B_\gamma } Q_{{B^-_\gamma }_{\textrm{rw}}} u, v\rangle _{B^+_\gamma }. \end{align*}

To obtain the second equality above, we used the first equality in (33). For the first and third equalities, we used (3), with $C=I$ , $\tilde C = B_\gamma ^-$ , and $C=B_\gamma$ , $\tilde C = B_\gamma ^+$ , respectively. Thus, under the assumption that $D_W=D_P$ and that $D_{B^+_\gamma }$ and $D_{B^-_\gamma }$ are invertible, the gradient flow of $f_{\varepsilon, 1}^\pm (\cdot ;\; B^+_\gamma, B^-_\gamma )$ with respect to the $B^+_\gamma$ -degree-weighted inner product is

\begin{align*} \frac{du}{dt} &= -L_{{B^+_\gamma }_{\textrm{rw}}} u - Q_{{B^-_\gamma }_{\textrm{rw}}} u + D_{B^+_\gamma }^{-1} D_{B_\gamma } Q_{{B^-_\gamma }_{\textrm{rw}}} u - \frac 1\varepsilon D_{B^+_1}^{-1} \Phi ' \circ u\\[5pt] &= -L_{{B^+_\gamma }_{\textrm{rw}}} u - Q_{{B^-_\gamma }_{\textrm{rw}}} u + (1-\gamma ) D_{B^+_\gamma }^{-1} D_W Q_{{B^-_\gamma }_{\textrm{rw}}} u - \frac 1\varepsilon D_{B^+_1}^{-1} \Phi ' \circ u. \end{align*}

Here, we used (20) again. We note that this gradient flow indeed equals the one from (32) if $\gamma =1$ .

In fact, in this paper we do not solve these Allen–Cahn-type equations directly, which could be accomplished, for example, by a convex–concave splitting technique as in [Reference Bertozzi and Flenner4] (see Luo and Bertozzi [Reference Luo and Bertozzi46] for an analysis of the scheme), but we use a related MBO scheme, which we introduce in Section 5. Before doing that, we first consider an extension of the graph GL functional to multiple clusters.

4.2. Multiclass clustering with graph Ginzburg–Landau functionals

The previous sections dealt with partitioning of the node set $V$ into (at most) two subsets; now we turn our attention to the case where we allow partitions of up to and including $K$ subsets, where $K\geq 2$ is fixed. We recall that we allow empty subsets in the partition. We base our approach on the method described in Garcia-Cardona et al. [Reference Garcia-Cardona, Merkurjev, Bertozzi, Flenner and Percus24] and Merkurjev et al. [Reference Merkurjev, Garcia-Cardona, Bertozzi, Flenner and Percus48].

To generalise the GL functional $f_{\varepsilon, \gamma }^\pm$ to the multiclass context, we require multiclass generalisations of the (signless) graph Dirichlet energies and of the double-well potential. For the Dirichlet energy, we generalise the term $\frac 12 \langle u, L_W u\rangle$ to

\begin{equation*} \frac 12 \langle U, L_W U\rangle = \frac 12 \sum _{k=1}^K \langle U_{*k}, L_W U_{*k} \rangle = \frac 14 \sum _{k=1}^K \sum _{i,j\in V} \omega _{ij} (U_{ik}-U_{jk})^2, \end{equation*}

where, in a slight overload of the inner product notation, for matrices $U,V \in \mathbb{R}^{|V|\times K}$ we have defined

(34) \begin{equation} \langle U, V\rangle \;:\!=\; \sum _{k=1}^K \langle U_{*k}, V_{*k}\rangle . \end{equation}

Similarly, we extend the $W$ -degree-weighted inner product to matrices $U,V \in \mathbb{R}^{|V|\times K}$ by:

(35) \begin{equation} \langle U, V\rangle _W \;:\!=\; \sum _{k=1}^K \langle U_{*k}, V_{*k}\rangle _W. \end{equation}

Hence,

\begin{equation*} \frac 12 \langle U, L_W U\rangle = \frac 12 \langle U, L_{W_{\textrm {rw}}} U\rangle _W. \end{equation*}

In a similar way, the signless graph Dirichlet energy $\frac 12 \langle u, Q_W u\rangle$ is generalised to

\begin{equation*} \frac 12 \langle U, Q_P U\rangle = \frac 14 \sum _{k=1}^K \sum _{i,j\in V} p_{ij} (U_{ik}+U_{jk})^2 = \frac 12 \langle U, Q_{P_{\textrm {rw}}} U\rangle _{P}. \end{equation*}

To generalise the double-well potential to a multiple-well potential, we recall from Section 3.2 that a partition of $V$ into $K$ subsets can be described by a matrix $U\in Pt(K)$ . If we write $U_{i*}$ for the $i^{\text{th}}$ row of $U$ , then $U_{i*} \in \{-1, 1\}^K$ as row vector and there exists a unique $k\in \{1, \ldots, K\}$ such that $U_{ik}=1$ . We introduce a notation for such vectors. For $k\in \{1, \ldots, K\}$ , let $e^{(k)} \in \{-1,1\}^K$ be the row vector that satisfies $e^{(k)}_k=1$ and, for all $l\in \{1, \ldots, K\}\setminus \{k\}$ , $e^{(k)}_l=-1$ . Now we define the multiple-well potential for vectors $w\in \mathbb{R}^K$ by:

(36) \begin{equation} \Phi _{\textrm{mul}}(w)\;:\!=\; \frac 12\left ( \prod _{k=1}^K \frac{1}{4} ||w - e^{(k)}||_1^2 \right ). \end{equation}

Given a matrix $U\in \mathbb{R}^{V\times |K|}$ , $\sum _{i\in V}\Phi _{\textrm{mul}}(U_{i*})$ is non-negative, and it is zero if and only if, for all $i\in V$ , there exists a $k\in \{1, \ldots, K\}$ such that $U_{i*}=e^{(k)}$ . In other words, $U \mapsto \sum _{i\in V}\Phi _{\textrm{mul}}(U_{i*})$ achieves its global minimum exactly at each element of $Pt(K)$ and thus generalises the double-well potential term $\sum _{i\in V} \Phi (u)$ , which achieves its global minimum exactly at all elements of $\mathcal{V}_{\textrm{bin}}$ .

This brings us to the following multiclass variant of the functional $f_{\varepsilon, \gamma }^\pm (\cdot ;\; W, P)$ for matrices $U\in \mathbb{R}^{|V|\times K}$ :

\begin{align*} \mathcal{F}_{\varepsilon, \gamma }^\pm (U;\; W, P) &\;:\!=\; \frac 18 \sum _{k=1}^K \sum _{i,j\in V} \omega _{ij} (U_{ik}-U_{jk})^2 + \frac \gamma 8 \sum _{k=1}^K \sum _{i,j\in V} p_{ij} (U_{ik}+U_{jk})^2 + \frac{1}{\varepsilon } \sum _{i\in V} \Phi _{\textrm{mul}}(U_{i*})\\[5pt] &= \frac 14 \langle U, L_W U\rangle + \frac \gamma 4 \langle U, Q_P U\rangle + \frac{1}{\varepsilon } \sum _{i\in V} \Phi _{\textrm{mul}}(U_{i*})\\[5pt] &= \frac 14 \langle U, L_{W_{\textrm{rw}}} U\rangle _W + \frac \gamma 4 \langle U, Q_{P_{\textrm{rw}}} U\rangle _P + \frac{1}{\varepsilon } \sum _{i\in V} \Phi _{\textrm{mul}}(U_{i*}). \end{align*}

As in the case of binary classification, we have a $\Gamma$ -convergence result for $\mathcal{F}_\varepsilon ^\pm$ .

Theorem 4.1. Let $K\in \mathbb{N}$ with $K\geq 2$ and $\gamma \in (0,\infty )$ . Define $\mathcal{F}_0^\pm \;:\; \mathbb{R}^{|V|\times K} \to \mathbb{R}\cup \{+\infty \}$ by

\begin{equation*} \mathcal {F}_{0,\gamma }^\pm (U;\; W, P) \;:\!=\; \begin {cases} \mathcal {Q}_{\textrm {mul},\gamma }(U;\; W, P), &\text {if } U \in Pt(K),\\[5pt] +\infty, &\text {otherwise}. \end {cases} \end{equation*}

Then $\mathcal{F}_{\varepsilon, \gamma }^\pm (\cdot ;\; W, P)$ $\Gamma$ -converges Footnote ²⁰ to $\mathcal{F}_{0,\gamma }^\pm (\cdot ;\; W, P)$ as $\varepsilon \downarrow 0$ .

Moreover, if $(\varepsilon _n) \subset (0,\infty )$ is a sequence such that $\varepsilon _n \downarrow 0$ as $n\to \infty$ and $(U_n) \subset Pt(K)$ is a sequence for which there exists a $C\gt 0$ such that, for all $n\in \mathbb{N}$ , $\mathcal{F}_{\varepsilon _n,\gamma }^\pm (U_n;\; W, P) \leq C$ , then there exists a converging subsequence of $(U_n)$ with limit in $Pt(K)$ .

Proof. Our proof largely follows Van Gennip and Bertozzi [Reference van Gennip and Bertozzi72] and Boyd et al. [Reference Boyd, Bae, Tai and Bertozzi7]. The strategy to prove $\Gamma$ -convergence for $\mathcal{F}_\varepsilon ^\pm$ is to first prove $\Gamma$ -convergence of the functional $\varphi _\varepsilon (U) \;:\!=\; \frac 1\varepsilon \sum _{i=1}^{|V|} \Phi _{\textrm{mul}} (U_{i*})$ to

\begin{equation*} \varphi _0 (U) \;:\!=\; \begin {cases} 0, & \text {if } U \in Pt(K),\\[5pt] +\infty, & \text {if } U \in \mathbb {R}^{|V|\times K}\setminus Pt(K), \end {cases} \end{equation*}

and then use the fact that $\Gamma$ -convergence is preserved under addition of continuous terms (see Dal Maso [Reference Maso47] or Braides [Reference Braides8]) such as the (signless) graph Dirichlet terms. The proof then concludes with the observation that by (25) and (26), for all $U\in Pt(K)$ ,

\begin{equation*} \frac 18 \sum _{k=1}^K \sum _{i,j\in V} \omega _{ij} (U_{ik}-U_{jk})^2 = \frac 14 \sum _{k=1}^K \sum _{i,j\in V} \omega _{ij} |U_{ik}-U_{jk}| = \frac {1}{2} \mathcal {TV}_W (U) \end{equation*}

and

\begin{equation*} \frac 18 \sum _{k=1}^K \sum _{i,j\in V} p_{ij} (U_{ik}+U_{jk})^2 = \frac 14 \sum _{k=1}^K \sum _{i,j\in V} p_{ij} |U_{ik}-U_{jk}| = \frac {1}{2} \mathcal {TV}_P^+ (U). \end{equation*}

Per definition, to establish $\Gamma$ -convergence of $\varphi _\varepsilon$ to $\varphi _0$ , we have to prove two statements:

• • Lower bound. For all sequences $(\varepsilon _n) \subset (0,\infty )$ that converge to zero, for all $U\in \mathbb{R}^{|V|\times K}$ , and for all sequences $(U_n) \subset \mathbb{R}^{|V|\times K}$ that converge to $U$ , it holds that

  \begin{equation*} \varphi _0(U) \leq \underset {n \to \infty }{\liminf }\, \varphi _{\varepsilon _n}(U_n). \end{equation*}

• • Recovery sequence. For all sequences $(\varepsilon _n) \subset (0,\infty )$ that converge to zero and for all $U\in \mathbb{R}^{|V|\times K}$ , there exists a sequence $(U_n) \subset \mathbb{R}^{|V|\times K}$ that converges to $U$ and such that

  \begin{equation*} \varphi _0(U) \geq \underset {n \to \infty }{\limsup }\, \varphi _{\varepsilon _n}(U_n). \end{equation*}

Let $(\varepsilon _n)$ be a positive sequence such that $\varepsilon _n \downarrow 0$ as $n \to \infty$ and let $U\in \mathbb{R}^{|V|\times K}$ .

To prove the lower bound condition, first we assume that $U\in Pt(K)$ . Since, for all $n\in \mathbb{N}$ , $\varphi _{\varepsilon _n}$ is non-negative, we find that, for all $n\in \mathbb{N}$

\begin{equation*} \varphi _0(U) = 0 \leq \varphi _{\varepsilon _n}(U_n) \end{equation*}

and the required $\liminf$ -inequality is satisfied.

If $U \in \mathbb{R}^{|V|\times K}\setminus Pt(K)$ , then $\varphi _0(U) = +\infty$ . Moreover, there exists a $i\in V$ such that, for all $k\in \{1, \ldots, K\}$ , $U_{i*}\neq e^{(k)}$ . Since $(U_{\varepsilon _n})$ converges to $U$ , it follows that there exists a radius $r\gt 0$ such that for $n$ large enough $(U_{\varepsilon _n})_{i*} \in \mathbb{R}^K \setminus \bigcup _{k=1}^K B(e^{(k)}, r)$ , where $B(e^{(k)}, r)$ denotes the open ball with respect to the $1$ -norm in $\mathbb{R}^K$ centred at $e^{(k)}$ with radius $r$ . This in turn implies that there exists a $\tilde C\gt 0$ such that, for $n$ large enough, $\Phi _{\textrm{mul}}(U_{i*}) \geq \tilde C$ . Hence,

\begin{equation*} \underset {n \to \infty }{\liminf } \, \frac 1{\varepsilon _n} \sum _{i\in V} \Phi _{\textrm {mul}}(U_{i*}) \geq \underset {n \to \infty }{\liminf } \, \frac {\tilde C}{\varepsilon _n} = +\infty, \end{equation*}

thus,

\begin{equation*} \varphi _0(U) = +\infty = \underset {n \to \infty }{\liminf } \, \varphi _{\varepsilon _n}(U_n). \end{equation*}

To prove existence of a recovery sequence, we note that the $\limsup$ -inequality is trivially true if $U \in \mathbb{R}^{|V|\times K}\setminus Pt(K)$ ; hence, we assume that $U\in Pt(K)$ . Let $(U_n)$ be the constant sequence with, for all $n \in \mathbb{N}$ , $U_n=U$ . Then, for all $n\in \mathbb{N}$ , $\varphi _{\varepsilon _n} (U_n) = 0$ . Since $\varphi _0$ is non-negative, the sequence $(U_n)$ is indeed a recovery sequence.

Thus, $\varphi _\varepsilon$ $\Gamma$ -converges to $\varphi _0$ and, by our argument above, $\mathcal{F}_\varepsilon ^\pm$ $\Gamma$ -converges to $\mathcal{F}_0^\pm$ .

To prove the equicoercivity statement in the second part of the theorem, we assume that $(\varepsilon _n) \subset (0,\infty )$ is a sequence such that $\varepsilon _n \downarrow 0$ as $n\to \infty$ and $(U_n) \subset Pt(K)$ is a sequence for which there exists a $C\gt 0$ such that, for all $n\in \mathbb{N}$ , $\mathcal{F}_{\varepsilon _n}^\pm (U_n) \leq C$ . In particular, this implies that for all $n\in \mathcal{N}$ and for all $i\in V$ ,

(37) \begin{equation} 0 \leq \Phi _{\textrm{mul}}((U_n)_{i*}) \leq C. \end{equation}

Let $i\in V$ . If there exists a $k\in \{1, \ldots, K\}$ such that the sequence $\left (\|(U_n)_{i*}-e^{(k)}\|_1\right )$ is unbounded in $\mathbb{R}$ , then there must also exist an $l\in \{1, \ldots, K\}$ such that the sequence $\left (\|(U_n)_{i*}-e^{(l)}\|_1\right )$ converges to zero; otherwise, the sequence $\left (\Phi _{\textrm{mul}}((U_n)_{i*})\right )$ would be unbounded. This is a contradiction; hence, the sequence $\left ((U_n)_{i*}\right )$ is contained in a bounded subset of $\mathbb{R}^K$ and thus the sequence $(U_n)$ is bounded in $\mathbb{R}^{|V|\times K}$ . By Bolzano–Weierstraß, this sequence has a converging subsequence. Denote its limit by $U_\infty$ . If $U_\infty \in \mathbb{R}^{|V|\times K}\setminus Pt(K)$ , then by the same argument as in the proof of the lower bound, we know that

\begin{equation*} \underset {n \to \infty }{\liminf } \, \frac 1{\varepsilon _n} \sum _{i\in V} \Phi _{\textrm {mul}}(U_{i*}) = +\infty, \end{equation*}

which contradicts (37). Hence, $U_\infty \in Pt(K)$ , which completes the proof.

As discussed in Section 4.1, the $\Gamma$ -convergence and equicoercivity results of Theorem4.1 imply that minimisers of $\mathcal{F}_\varepsilon ^\pm$ converge to minimisers of $\mathcal{F}_0^\pm$ as $\varepsilon \downarrow 0$ .

In a completely analogous way to what we did at the end of Section 4.1, we may now derive gradient flows of $\mathcal{F}_\varepsilon ^\pm$ with respect to the Euclidean or $W$ -degree-weighted inner product for matrices. Because we will not actually use these gradient flows to find minimisers, but use MBO schemes (see Section 5.3) instead, we leave out the details of the derivation and only mention that formallyFootnote ²¹ we recover, for all $k\in \{1, \ldots, K\}$ , the Allen–Cahn-type equations:

\begin{equation*} \frac {dU_{*k}}{dt} = -L_W U_{*k} - Q_P U_{*k} - \frac 1\varepsilon (\mathcal {D}\Phi _{\textrm {mul}} \circ U)_{*k} \ \ \text {and} \ \ \frac {dU_{*k}}{dt} = -L_{W_{\textrm {rw}}} U_{*k} - Q_{P_{\textrm {rw}}} U_{*k} - \frac 1\varepsilon (D_W^{-1} \mathcal {D}\Phi _{\textrm {mul}} \circ U)_{*k}, \end{equation*}

where $(\mathcal{D}\Phi _{\textrm{mul}} \circ U)_{ik} \;:\!=\; \partial _k \Phi _{\textrm{mul}}(U_{i*})$ , with $\partial _k$ the partial derivative operator with respect to the $k^{\text{th}}$ variable.

Remark 4.2. As in Section 4.1, we can recover similar results as above, if we consider the functionals $\mathcal{F}_{\varepsilon, 1}^\pm (\cdot ;\; B^+_\gamma, B^-_\gamma )$ . Their $\Gamma$ -limit for $\varepsilon \downarrow 0$ is $\mathcal{F}_{0,1}^\pm (\cdot ;\; B^+_\gamma, B^-_\gamma )$ and their gradient flow with respect to the Euclidean inner product is, for all $k\in \{1, \ldots, L\}$ ,

\begin{equation*} \frac {dU_{*k}}{dt} = -L_{B^+} U_{*k} - Q_{B^-} U_{*k} - \frac 1\varepsilon (\mathcal {D}\Phi _{\textrm {mul}} \circ U)_{*k} . \end{equation*}

The gradient flow with respect to the $B^+_\gamma$ -degree-weighted inner product is, for all $k\in \{1, \ldots, L\}$ ,

\begin{align*} \frac{dU_{*k}}{dt} &= -L_{B^+_{\textrm{rw}}} U_{*k} - Q_{B^-_{\textrm{rw}}} U_{*k} + D_{B^+_\gamma }^{-1} D_{B_\gamma } Q_{{B^-_\gamma }_{\textrm{rw}}} u - \frac 1\varepsilon (D_{B^+}^{-1} \mathcal{D}\Phi _{\textrm{mul}} \circ U)_{*k}\\[5pt] &= -L_{B^+_{\textrm{rw}}} U_{*k} - Q_{B^-_{\textrm{rw}}} U_{*k} + (1-\gamma ) D_{B^+_\gamma }^{-1} D_W Q_{{B^-_\gamma }_{\textrm{rw}}} u - \frac 1\varepsilon (D_{B^+}^{-1} \mathcal{D}\Phi _{\textrm{mul}} \circ U)_{*k}. \end{align*}

We recall that we have used $D_W=D_P$ and invertibility of $D_{B^+_\gamma }$ and of $D_{B^-_\gamma }$ (as is guaranteed if assumptions of Corollary 3.2, part 1 are satisfied).

We also note that the $\Gamma$ -convergence and equicoercivity proofs in Theorem 4.1 mostly depended on the potential term and can thus easily be reproduced with the matrices $B^+$ and $B^-$ replacing $W$ and $P$ , respectively.

5.1. MBO schemes for binary community detection

The MBO scheme was adapted to graphs by Merkurjev et al. [Reference Merkurjev, Kostić and Bertozzi49] with the goal of (approximately) minimising the GL functional $f_\varepsilon$ . It consists of iteratively performing short-time graph diffusion by solving $\frac{du}{dt} = -L u$ , followed by a hard thresholding step to mimic the drive towards the wells of $\Phi$ from the non-linear reaction term in the gradient flow. Here, we take $L\in \{L_W, L_{W_{\textrm{rw}}}, L_{W_{\textrm{sym}}}\}$ to allow for different variants of the scheme. An MBO-type scheme is also employed by Hu et al. in [Reference Hu, Laurent, Porter and Bertozzi33] and Boyd et al. in [Reference Boyd, Bae, Tai and Bertozzi7] for modularity optimisation.

Unless specified differently, we assume that

(38) \begin{align} L_{\textrm{mix}}\in &\left \{L_W+\gamma Q_P, L_{W_{\textrm{sym}}}+ \gamma Q_{P_{\textrm{sym}}}, L_{W_{\textrm{rw}}}+\gamma Q_{P_{\textrm{rw}}},\right . \notag \\[5pt] &\hspace{0.3cm} \left . L_{B^+_\gamma }+Q_{B^-_\gamma }, L_{{B^+_\gamma }_{\textrm{sym}}}+Q_{{B^-_\gamma }_{\textrm{sym}}}, L_{{B^+_\gamma }_{\textrm{rw}}}+Q_{{B^-_\gamma }_{\textrm{rw}}}-D_{B^+_\gamma }^{-1} D_{B_\gamma } Q_{{B^-_\gamma }_{\textrm{rw}}}\right \}, \end{align}

This assumption is motivated by the linear operators that appear in the gradient flows of Sections 4.1 and 4.2.

For later reference, we note that

\begin{equation*} I-D_{B^+_\gamma }^{-1} D_{B_\gamma } = I - D_{B^+_\gamma }^{-1} (D_{B^+_\gamma } - D_{B^-_\gamma }) = D_{B^+_\gamma }^{-1}D_{B^-_\gamma } \end{equation*}

and thus for the last choice of $L_{\textrm{mix}}$ in (38), we can also write

(39) \begin{align} L_{{B^+_\gamma }_{\textrm{rw}}}+Q_{{B^-_\gamma }_{\textrm{rw}}}-D_{B^+_\gamma }^{-1} D_{B_\gamma } Q_{{B^-_\gamma }_{\textrm{rw}}} &= L_{{B^+_\gamma }_{\textrm{rw}}}+(I-D_{B^+_\gamma }^{-1} D_{B_\gamma }) Q_{{B^-_\gamma }_{\textrm{rw}}} = L_{{B^+_\gamma }_{\textrm{rw}}}+D_{B^+_\gamma }^{-1}D_{B^-_\gamma } Q_{{B^-_\gamma }_{\textrm{rw}}}\notag \\[5pt] &= L_{{B^+_\gamma }_{\textrm{rw}}}+D_{B^+_\gamma }^{-1} Q_{B^-_\gamma }. \end{align}

This form is often easier to work with, but it hides partly the fact that it is the random walk variant of $L_{\textrm{mix}}$ for the split of $B_\gamma$ into a positive and negative part.

We recall from Section 2.2 that connectedness of $G$ implies invertibility of $D_W$ and, if the null model is such that $D_W=D_P$ , also the invertibility of $D_P$ , so that the choices of $L_{\textrm{mix}}$ that involve inverses of those matrices are well defined. For the choices that require inverses of $D_{B^+_\gamma }$ or $D_{B^-_\gamma }$ we need the additional assumption that these matrices are invertible, which holds, for example, in the situation of Corollary 3.2.

Adapting the idea in [Reference Merkurjev, Kostić and Bertozzi49] to the equations in (29), we propose the following MBO-type scheme.

Binary $L_{\textrm{mix}}$ modularity MBO scheme___________________________________________________

• • Initialise. Choose an initial condition $u^0 \in \mathcal{V}_{\textrm{bin}}$ , a ‘time step’ $\tau \gt 0$ , and $L_{\textrm{mix}}$ as in (38).

• • Step $n+1$ : linear dynamics. Solve the equation $\frac{du}{dt} = -L_{\textrm{mix}} u$ on $(0,\tau ]$ with initial condition $u(0)=u^n$ .

• • Step $n+1$ : threshold. Define, for all $i\in V$ , $u^{n+1}_i \;:\!=\; \begin{cases} -1, &\text{if } u_i(\tau ) \lt 0,\\[5pt] 1, &\text{if } u_i(\tau ) \geq 0. \end{cases}$

• • Stop. Stop the scheme when a stopping condition or predetermined number of steps has been achieved.

__________________________________________________________________________________

To indicate explicitly the dependence on the choice of $L_{\textrm{mix}}$ , we call this the binary $L_{\textrm{mix}}$ modularity MBO (MMBO) scheme. We note that, besides the choices of $L_{\textrm{mix}}$ that follow from equations (29), (30), and (32), we also allow the variant choicesFootnote ²² $L_{\textrm{mix}}\in \{L_{W_{\textrm{sym}}}+\gamma Q_{P_{\textrm{sym}}}, L_{{B^+_\gamma }_{\textrm{sym}}}+Q_{{B^-_\gamma }_{\textrm{sym}}}\}$ .

We briefly mention here that we will choose the value for $\tau$ in the MMBO scheme via the method presented in [Reference Boyd, Bae, Tai and Bertozzi7]. It is an effective strategy and requires less manual adjustment of $\tau$ than other approaches. We give more details in Section 6.1.

5.2. Numerical schemes for binary MMBO

We employ two different numerical methods to implement the binary MMBO scheme, where the difference is found in how the methods solve the linear-dynamics step. One method uses (a truncation of) the closed-form solution [Reference Hale29, Reference Hall30], while the other uses an implicit finite-difference Euler scheme [Reference Butcher14]. In both methods, we need (the leading) eigenvalues and eigenvectors of $L_{\textrm{mix}}$ .

5.2.1. Closed-form matrix exponential solution

A closed-form solution at $t=\tau$ of the equation $\frac{du}{dt} = -L_{\textrm{mix}} u$ with initial condition $u(0)=u^n$ is given by [Reference Hale29]:

\begin{equation*} u(\tau ) = e^{-\tau L_{\textrm {mix}}} u^n. \end{equation*}

Here, $e^{-\tau L_{\textrm{mix}}}$ is the matrix exponential, which is defined by its series expansion:

\begin{equation*} e^{-\tau L_{\textrm {mix}}} \;:\!=\; I + \sum _{k=1}^\infty \frac 1{k!} (\!-\tau L_{\textrm {mix}})^k. \end{equation*}

If $L_{\textrm{mix}}$ has real eigenvalues $\lambda _1 \leq \ldots \leq \lambda _{|V|}$ with corresponding, linearly independent, eigenvectors $\xi _i$ , $i\in \{1, \ldots, |V|\}$ , then $e^{-\tau L_{\textrm{mix}}}$ has eigenvalues $e^{-\tau \lambda _i}$ with the same eigenvectors. Hence,

(40) \begin{equation} u(\tau ) = \sum _{i=1}^{|V|} c_i e^{-\tau \lambda _i} \xi _i, \end{equation}

for appropriate coefficients $c_i\in \mathbb{R}$ .

The following lemma shows that for each of our choices, $L_{\textrm{mix}}$ is indeed diagonalisable and thus has $|V|$ linearly independent eigenvectors.

Lemma 5.1. Assume $L_{\textrm{mix}}$ satisfies (38) and $D_P$ is invertible, where this is needed for $L_{\textrm{mix}}$ to be well defined. Then $L_{\textrm{mix}}$ has $|V|$ (possibly repeated) non-negative real eigenvalues with corresponding linearly independent normalised eigenvectors, where the normalisation is specified in each of the cases below.

Denote by $X$ a matrix having these eigenvectors as columns (in any order) and by $\Lambda$ , the diagonal matrix containing the corresponding eigenvalues in the same order. Then $L_{\textrm{mix}}$ is real diagonalisable (in the standard Euclidean inner product structure), that is, $L_{\textrm{mix}} = X \Lambda X^{-1}$ .

1. (a) Let $L_{\textrm{mix}}\in \{L_W+\gamma Q_P, L_{W_{\textrm{sym}}}+\gamma Q_{P_{\textrm{sym}}}, L_{B^+_\gamma }+Q_{B^-_\gamma }, L_{{B^+_\gamma }_{\textrm{sym}}}+Q_{{B^-_\gamma }_{\textrm{sym}}}\}$ . In the case that $L_{\textrm{mix}}=L_{{B^+_\gamma }_{\textrm{sym}}}+Q_{{B^-_\gamma }_{\textrm{sym}}}$ , assume that $D_{B^+_\gamma }$ and $D_{B^-_\gamma }$ are invertible. Take for the columns of $X$ eigenvectors of $L_{\textrm{mix}}$ with unit Euclidean norm. Then $X$ is orthogonal, that is, $X^{-1} = X^T$ .

2. (b) Let $L_{\textrm{mix}} = L_{W_{\textrm{rw}}}+\gamma Q_{P_{\textrm{rw}}}$ . Assume that the null model is such that $D_P=D_W$ . Moreover, take for the columns of $X$ eigenvectors of $L_{\textrm{mix}}$ with unit norm with respect to the $W$ -degree-weighted inner product from (2) . Then $X^{-1} = \tilde X^T D_W^{\frac 12}$ , where $\tilde X$ is the orthogonal matrix containing the Euclidean-normalised eigenvectors of $L_{W_{\textrm{sym}}}+\gamma Q_{P_{\textrm{sym}}}$ as columns, in the same order as the eigenvalues in $\Lambda$ .

3. (c) Assume that $D_{B^+_\gamma }$ and $D_{B^-_\gamma }$ are invertible. Let $L_{\textrm{mix}} = L_{{B^+_\gamma }_{\textrm{rw}}}+Q_{{B^-_\gamma }_{\textrm{rw}}}-D_{B^+_\gamma }^{-1} D_{B_\gamma } Q_{{B^-_\gamma }_{\textrm{rw}}}$ . Moreover, take for the columns of $X$ eigenvectors of $L_{\textrm{mix}}$ with unit norm with respect to the $B^+_\gamma$ -degree-weighted inner product from (31) . Then $X^{-1} = \tilde X^T D_{B^+_\gamma }^{\frac 12}$ , where $\tilde X$ is the orthogonal matrix containing the Euclidean-normalised eigenvectors of $L_{{B^+_\gamma }_{\textrm{sym}}} + D_{B^+_\gamma }^{-\frac 12} Q_{B^-_\gamma } D_{B^+_\gamma }^{-\frac 12}$ as columns, in the same order as the eigenvalues in $\Lambda$ .

Proof. The proof is given in Appendix B.

Remark 5.2. We note that the proof of Lemma 5.1, part (b) also establishes that $L_{W_{\textrm{rw}}}+\gamma Q_{P_{\textrm{rw}}}$ and $L_{W_{\textrm{sym}}}+\gamma Q_{P_{\textrm{sym}}}$ have the same eigenvalues. Moreover, the proof of Lemma 5.1, part (c) establishes that $L_{B^+_{\textrm{rw}}}+Q_{B^-_{\textrm{rw}}}-D_{B^+_\gamma }^{-1} D_{B_\gamma } Q_{{B^-_\gamma }_{\textrm{rw}}}$ and $L_{{B^+_\gamma }_{\textrm{sym}}} + D_{B^+_\gamma }^{-\frac 12} Q_{B^-_\gamma } D_{B^+_\gamma }^{-\frac 12}$ have the same eigenvalues.

Remark 5.3. If $\gamma =1$ and $D_W=D_P$ in case (c) of Lemma 5.1, the situation simplifies, since then $D_{B_1}=0$ and $D_{B^+_1} = D_{B^-_1}$ and thus

\begin{equation*} D_{B^+_1}^{\frac 12} L_{\textrm {mix}} D_{B^+_1}^{-\frac 12} = D_{B^+_1}^{\frac 12} \left (L_{{B^+_1}_{\textrm {rw}}}+Q_{{B^-_1}_{\textrm {rw}}}\right ) D_{B^+_1}^{-\frac 12} = L_{{B^+_1}_{\textrm {sym}}}+Q_{{B^-_1}_{\textrm {sym}}}. \end{equation*}

Part (a) of Lemma 5.1 then implies that $L_{\textrm{mix}}$ is real diagonalisable and $X=\tilde D_{B^+_1}^{-\frac 12} X$ , where $\tilde X$ is the matrix of Euclidean-normalised eigenvectors of $L_{{B^+_1}_{\textrm{sym}}}+Q_{{B^-_1}_{\textrm{sym}}}$ . Moreover, $L_{\textrm{mix}}$ has non-negative eigenvalues since both $L_{{B^+_1}_{\textrm{sym}}}$ and $Q_{{B^-_1}_{\textrm{sym}}}$ , and thus also their sum, are positive semidefinite with respect to the Euclidean inner product by Lemma 2.1 part (b) and Lemma 2.2 part (b). We note that this result is consistent with the result in part (c) of Lemma 5.1, since, in this case, $L_{{B^+_1}_{\textrm{sym}}} + D_{B^+_1}^{-\frac 12} Q_{B^-_1} D_{B^+_1}^{-\frac 12} = L_{{B^+_1}_{\textrm{sym}}}+Q_{{B^-_1}_{\textrm{sym}}}$ .

Remark 5.4. In case (b) of Lemma 5.1, we know that $L_{\textrm{mix}}$ is self-adjoint with respect to the inner product $\langle \cdot, \cdot \rangle _W$ (see Section 2.2) and thus we expect the columns of $X$ to be orthonormal with respect to this inner product. Indeed, we compute

\begin{equation*} X^T D_W X = \tilde X^T D_W^{-\frac 12} D_W^{\frac 12} D_W^{-\frac 12} \tilde X = \tilde X^T \tilde X = I. \end{equation*}

A similar conclusion holds in case (c) of the lemma, yet now with respect to the $B^+_\gamma$ -degree-weighted inner product from (31) as can easily be seen by replacing $W$ by $B^+_\gamma$ in the computation above. This is not surprising, as we know from the proof of the lemma, that $D_{B^+_\gamma }^{\frac 12} L_{\textrm{mix}} D_{B^+_\gamma }^{-\frac 12}$ is symmetric and thus is self-adjoint with respect to the Euclidean inner product; hence, $L_{\textrm{mix}}$ is self-adjoint with respect to the $B^+_\gamma$ -degree-weighted inner product, since, for all $u,v\in \mathcal{V}$ ,

\begin{align*} \langle L_{\textrm{mix}} u, v \rangle _{B^+_\gamma } &= \langle D_{B^+_\gamma } L_{\textrm{mix}} u, v \rangle = \langle D_{B^+_\gamma }^{\frac 12} (D_{B^+_\gamma }^{\frac 12} L_{\textrm{mix}} D_{B^+_\gamma }^{-\frac 12}) (D_{B^+_\gamma }^{\frac 12} u),v\rangle \\[5pt] &= \langle D_{B^+_\gamma }^{\frac 12} L_{\textrm{mix}} D_{B^+_\gamma }^{-\frac 12} (D_{B^+_\gamma }^{\frac 12} u), D_{B^+_\gamma }^{\frac 12} v\rangle = \langle D_{B^+_\gamma }^{\frac 12} u, D_{B^+_\gamma }^{\frac 12} L_{\textrm{mix}} D_{B^+_\gamma }^{-\frac 12} (D_{B^+_\gamma }^{\frac 12} v)\rangle = \langle u, L_{\textrm{mix}} v \rangle _{B^+_\gamma }. \end{align*}

For the next lemma, we recall that $\gamma \gt 0$ .

Lemma 5.5. $\quad$

1. (a) Let $L_{\textrm{mix}}\in \big \{L_W+\gamma Q_P, L_{W_{\textrm{sym}}}+ \gamma Q_{P_{\textrm{sym}}}\big \}$ . Assume $D_P$ is invertible, where this is needed to define $L_{\textrm{mix}}$ . If the null model is such that the matrix $P$ has at least one positive entry, then the eigenvalues of $L_{\textrm{mix}}$ are positive.

2. (b) Let $L_{\textrm{mix}} = L_{W_{\textrm{rw}}}+\gamma Q_{P_{\textrm{rw}}}$ . If the null model is such that $D_W=D_P$ , then the eigenvalues of $L_{\textrm{mix}}$ are positive.

3. (c) Let $L_{\textrm{mix}} \in \big \{L_{B^+_\gamma }+Q_{B^-_\gamma }, L_{{B^+_\gamma }_{\textrm{sym}}}+Q_{{B^-_\gamma }_{\textrm{sym}}}, L_{{B^+_\gamma }_{\textrm{rw}}}+Q_{{B^-_\gamma }_{\textrm{rw}}}-D_{B^+_\gamma }^{-1} D_B Q_{{B^-_\gamma }_{\textrm{rw}}}\big \}$ and assume that $D_{B^+_\gamma }$ and $D_{B^-_\gamma }$ are invertible, in those cases where this is needed to define $L_{\textrm{mix}}$ . Assume one of the following conditions is satisfied:

   1. (i) the graph with adjacency matrix $B^+_\gamma$ is connected and the matrix $B^-_\gamma$ has at least one positive entry;

   2. (ii) the graph with adjacency matrix $B^-_\gamma$ is connected and the matrix $B^+_\gamma$ has at least one positive off-diagonal entry $(b^+_\gamma )_{ij}$ ; moreover, there exists a path along an odd number of edges from $i$ to $j$ in the graph with adjacency matrix $B^-_\gamma$ ; or

   3. (iii) the matrix $B^-_\gamma$ has positive diagonal entries.

   Then the eigenvalues of $L_{\textrm{mix}}$ are positive.

Proof. The proof is given in Appendix B.

Remark 5.6. The assumptions in the first two cases of Lemma 5.5 are satisfied for the NG null model.

If the graph defined by the adjacency matrix $W$ is connected and has no self loops and $P$ has positive diagonal elements, as is again the case for the NG null model, then, for all $i\in V$ , $(b^-_\gamma )_i = \gamma p_{ii}\gt 0$ . So (if also $D_{B^+_\gamma }$ and $D_{B^-_\gamma }$ are invertible), the conditions in the third part of the lemma are satisfied.

Remark 5.7. In general, it cannot be expected that the eigenvalues of $L_{\textrm{mix}}$ are positive without additional assumptions on the graph. For example, if $L_{\textrm{mix}} = L_{B^+_\gamma } + Q_{B^-_\gamma }$ where

\begin{equation*} B = \begin {pmatrix} 0\;\;\;\; & 1\;\;\;\; & -1\;\;\;\; & 0\\[5pt] 1\;\;\;\; & 0\;\;\;\; & -1\;\;\;\; & -1\\[5pt] -1\;\;\;\; & -1\;\;\;\; & 0\;\;\;\; & 1\\[5pt] 0\;\;\;\; & -1\;\;\;\; & 1\;\;\;\; & 0 \end {pmatrix}, \end{equation*}

it can be checked that $u\in \mathcal{V}$ with $u_1=u_2=1$ and $u_3=u_4=-1$ is an eigenfunction (eigenvector) with eigenvalue zero.

We now return to expression (40). When we are in case (a) of Lemma 5.1, the normalised eigenvectors $\xi _i$ of $L_{\textrm{mix}}$ are orthonormal, and thus for the coefficients in (40) we compute $ c_i = \langle \xi _i, u^n\rangle .$

In case (b) of the lemma, we know by Remark5.4 that $L_{\textrm{mix}}$ has eigenvectors $\xi _i$ that are orthonormal with respect to the inner product $\langle \cdot, \cdot \rangle _W$ , and thus $c_i = \langle \xi _i, u^n\rangle _W$ instead.

Finally, in case (c) Remark5.4 tells us that a similar conclusion holds, if we use the $B^+_\gamma$ -degree-weighted inner product from (31): $ c_i = \langle \xi _i, u^n\rangle _{B^+_\gamma }.$

Writing (40) fully in matrix form, we thus obtain

(41) \begin{align} u(\tau ) &= \left \{\begin{array}{lr} X e^{-\tau \Lambda } X^T u^n, &\text{in case\,(a) of Lemma\,5.1},\\[5pt] X e^{-\tau \Lambda } X^T D_W u^n = D_W^{-\frac 12} \tilde X e^{-\tau \Lambda } \tilde X^T D_W^{\frac 12} u^n, &\text{in case\,(b) of Lemma\,5.1},\\[5pt] X e^{-\tau \Lambda } X^T D_{B^+_\gamma } u^n = D_{B^+_\gamma }^{-\frac 12} \tilde X e^{-\tau \Lambda } \tilde X^T D_{B^+_\gamma }^{\frac 12} u^n, &\text{in case\,(c) of Lemma\,5.1} \end{array}\right \}\notag \\[5pt] &= X e^{-\tau \Lambda } X^{-1} u^n. \end{align}

We recall that the matrices $X$ are not the same in each case, and neither are the matrices $\tilde X$ .

5.2.2. Implicit Euler finite-difference discretisation

Next, we take a look at the implicit Euler finite-difference discretisation. We discretise the time domain $(0, \tau ]$ into $N_t\in \mathbb{N}$ intervals $(t_{k-1}, t_k]$ ( $k\in \{1, \ldots, N_t\}$ ) of equal length, thus $t_k = \frac{k}{N_t} \tau = k \delta t$ , where $\delta t \;:\!=\; \frac \tau{N_t}\gt 0$ . Approximating $u$ at the points $t_k$ by $u(t_k) \approx u^k \in \mathbb{R}$ (with $u^0=u(0)$ ), the discretisation of the equation $\frac{du}{dt} = -L_{\textrm{mix}} u$ is given by

\begin{equation*} \frac {u^k - u^{k-1}}{\delta t}=-L_{\textrm {mix}} u^k. \end{equation*}

This we can rewrite as

(42) \begin{equation} u^k = \left ( I + \delta t L_{\textrm{mix}} \right )^{-1} u^{k-1}, \end{equation}

where $I$ is the identity matrix of the appropriate size.

From Lemma 5.1, we recall that $X$ is a matrix containing the normalisedFootnote ²³ eigenvectors of $L_{\textrm{mix}}$ as columns and $\Lambda$ a diagonal matrix containing the corresponding eigenvalues of $L_{\textrm{mix}}$ (in the same order as the eigenvectors). By the same lemma, we know that $L_{\textrm{mix}}$ is diagonalisable as $L_{\textrm{mix}} = X \Lambda X^{-1}$ . Hence, (42) can be written as

\begin{equation*} u^k = \left [ X(I + \delta t \Lambda ) X^{-1} \right ]^{-1} u^{k-1} = \left [ X (I + N_t^{-1}\tau \Lambda )^{-1} X^{-1} \right ]^k u^0. \end{equation*}

In case (a) of Lemma 5.1, $X$ is orthogonal and thus its inverse can be computed as $X^{-1}=X^T$ . In cases (b) and (c), we have seen that $X^{-1}= \tilde X^T D^{\frac 12}$ , for the appropriate orthogonal eigenvector matrix $\tilde X$ and appropriate degree matrix $D$ .

In the linear-dynamics step of the MMBO scheme, we are interested in $u(\tau )$ , for which we find

(43) \begin{equation} u(\tau ) \approx u^{N_t} = \left [ X (I + N_t^{-1}\tau \Lambda )^{-1} X^{-1} \right ]^{N_t} u^n, \end{equation}

where (in a slight abuse of superscript notation) we recall that $u^0=u(0)=u^n$ , with superscript $0$ indicating the initial condition for the Euler finite-difference scheme, but superscript $n$ indicating the iteration number of the MMBO scheme. We recognise an approximation of the closed-form solution from (41) in (43).

Similarly to what we did at the end of Section 5.2.1, $u^{N_t}$ can be written as

\begin{align*} u^{N_t} = \begin{cases} \left [X \left (I+ \frac{\tau }{N_t} \Lambda \right )^{-1} X^T \right ]^{N_t} u^n, &\text{in case\,(a) of Lemma\,5.1},\\[5pt] \left [D_W^{-\frac 12} \tilde X \left (I+ \frac{\tau }{N_t} \Lambda \right )^{-1} \tilde X^T D_W^{\frac 12} \right ]^{N_t} u^n, &\text{in case\,(b) of Lemma\,5.1},\\[5pt] \left [D_{B^+}^{-\frac 12} \tilde X \left (I+ \frac{\tau }{N_t} \Lambda \right )^{-1} \tilde X^T D_{B^+}^{\frac 12} \right ]^{N_t} u^n, &\text{in case\,(c) of Lemma\,5.1}, \end{cases} \end{align*}

where $X$ and $\tilde X$ are the appropriate orthogonal eigenvector matrices for each case of Lemma 5.1.

5.3. Multiclass MMBO scheme

Just as the binary $L_{\textrm{mix}}$ MMBO scheme from Section 5.1 was inspired by the Allen–Cahn-type equations from Section 4.1, so we can base a multiclass $L_{\textrm{mix}}$ MMBO scheme on the multiclass Allen–Cahn-type equations from Section 4.2. We recall that $K\geq 2$ is fixed.

Multiclass $L_{\textrm{mix}}$ modularity MBO scheme________________________________________________

• • Initialise. Choose an initial condition $U^0 \in Pt(K)$ , a ‘time step’ $\tau \gt 0$ , and $L_{\textrm{mix}}$ as in (38).

• • Step $n+1$ : linear dynamics. Solve the equation $\frac{dU}{dt} = -L_{\textrm{mix}} U$ on $(0,\tau ]$ with initial condition $U(0)=U^n$ .

• • Step $n+1$ : threshold. Define, for all $i\in V$ , $U^{n+1}_{i*} \;:\!=\; e^{(k^*)}$ , whereFootnote ²⁵

  \begin{equation*} k^* \in \operatorname{argmax}_{k\in \{1, \ldots, K\}} U_{ik}(\tau). \end{equation*}

• • Stop. Stop the scheme when a stopping condition or predetermined number of steps has been achieved.

__________________________________________________________________________________

The linear-dynamics step is a straightforward generalisation of the analogous step in the binary algorithm. The threshold step now needs to take into account that there are $K\geq 2$ clusters to choose from. We assign each node to that cluster which is represented by the highest value in the row vector $U(\tau )_{i*}$ . In [Reference Cucuringu, Pizzoferrato and van Gennip20, Appendix A.3], it is proven that this procedure corresponds to first projecting the vector $U(\tau )_{i*}$ onto the $(K-1)$ -simplex

(44) \begin{equation} \mathfrak{S}(K) \;:\!=\; \left \{ w \in [\!-1,1]^K \;:\; \sum _{k=1}^K w_k = 2-K\right \} \end{equation}

and then determining the nearest vertex (or vertices) of the simplex (i.e., nearest vector $e^{(k)}$ ) to this projected vector.

5.4. Numerical schemes for multiclass MMBO

We briefly discuss the generalisations of the numerical schemes that we encountered in Section 5.2 to the multiclass case. An in-depth look at the resulting algorithms follows in Section 6.

5.4.1. Closed-form matrix exponential solution

If the relevant (for the case at hand) assumptions from Lemma 5.1 are satisfied, then the closed-form solution of the multiclass-linear-dynamics step is a straightforward generalisation of the solution in the binary case, since each column of $U$ satisfies the same linear system of ODEs and there is no coupling between the dynamics of different columns. Thus, $U(\tau ) = e^{-\tau L_{\textrm{mix}}} U^n$ and from there the same arguments as in Section 5.2.1 lead to

(45) \begin{equation} U(\tau ) = X e^{-\tau \Lambda } X^{-1} U^n, \end{equation}

where we recall that the meaning of $X$ depends on which of the three cases of Lemma 5.1 $L_{\textrm{mix}}$ satisfies.

5.4.2. Implicit Euler finite-difference discretisation

For the same reasons as mentioned in Section 5.4.1, the Euler finite-difference scheme from Section 5.2.2 also straightforwardly generalises from the binary to the multiclass setting, assuming that the relevant assumptions of Lemma 5.1 are satisfied. Hence, we find, from (43), that

(46) \begin{equation} U(\tau ) \approx \left [ X (I + N_t^{-1}\tau \Lambda )^{-1} X^{-1} \right ]^{N_t} U^n, \end{equation}

where we recall that $N_t\in \mathbb{N}$ is the number of steps in the finite-difference discretisation of the interval $[0,\tau ]$ .

6.1. Choice of the time step

To avoid having to manually fine-tune the value of the time step $\tau$ in our algorithms, we follow the method proposed by Boyd et al. [Reference Boyd, Bae, Tai and Bertozzi7, Section 4.5].

Two main considerations are of importance in the choice of $\tau$ . If the value is chosen too small, then the linear-dynamics step of the MBO scheme will change the initial condition very little and the threshold step will return the initial condition again. In other words, the initial condition will be stationary under the MBO scheme. If the value of $\tau$ is chosen too large, then (because all eigenvalues of $L_{\textrm{mix}}$ are positive) the result of one application of the linear-dynamics step will be a function which is close to zero on all nodes of the graph. Only the structure contained in the mode(s) with the smallest eigenvalue(s) is retained, which is typically not sufficient to find optimal communities. Moreover, to have a threshold step which is robust to noise, we prefer the values on the nodes to be clearly separated and not all clustered together near zero. Thus, we need to choose a value of $\tau$ which is neither too small nor too large. The details of what is ‘too small’ and ‘too large’ depend on the structure of the graph and the choice of initial condition. In [Reference Boyd, Bae, Tai and Bertozzi7], a specific choice of $\tau$ is suggested for the particular MBO-type scheme that is employed in that paper. Adapting a method from Van Gennip et al. [Reference van Gennip, Guillen, Osting and Bertozzi73], upper ( $\tau _{\textrm{upp}}$ ) and lower ( $\tau _{\textrm{low}}$ ) bounds are establishedFootnote ²⁶ by [Reference Boyd, Bae, Tai and Bertozzi7] and their geometric mean $\sqrt{\tau _{\textrm{low}}\tau _{\textrm{upp}}}$ is used as the value of $\tau$ . For the lower bound, $\tau _{\textrm{low}}$ is computed for $K=2$ (even in cases where $K\neq 2$ ). An explicit numerical value for general $K$ is harder to obtain and [Reference Boyd, Bae, Tai and Bertozzi7] expects (without proof) that the case $K=2$ presents the worst-case scenario. In the next lemma, we adapt the method from [Reference Boyd, Bae, Tai and Bertozzi7] to our MBO scheme and give upper and lower bounds on $\tau$ (which are not expected to be sharp but are efficiently computable).

Algorithm 1. The MMBO scheme using the closed-form solution of the linear-dynamics step

We follow [Reference Boyd, Bae, Tai and Bertozzi7] in choosing $\tau = \sqrt{\tau _{\textrm{low}}\tau _{\textrm{upp}}}$ .

6.2. MMBO scheme using the closed-form solution of the linear-dynamics step

We recall from Section 5.4.1 that, to employ the closed-form solution of the linear-dynamics step to compute the multiclass $L_{\textrm{mix}}$ MMBO scheme, we need to compute the expression in (45). The resulting algorithm is summarised in Algorithm1. The command $random.sample(K, N)$ returns a list of length $N$ , with uniformly random sampling from $1$ to $K$ of $N$ distinct elements ( $N \leq K$ ).

6.3. Alternative variant of the MMBO scheme

From Section 5.4.2, we recall that for the MMBO scheme with the Euler finite-difference scheme, we have to solve equation (46).

We use the same thresholding step and stopping conditions for this alternative variant as we do for the MMBO scheme in Algorithm1. The MMBO scheme using the Euler finite-difference discretisation is summarised in Algorithm2. The main difference between Algorithm1 and Algorithm2 is the diffusion step, as seen in (45) and (46).

Algorithm 2. The MMBO scheme using the Euler finite-difference discretisation

6.4. Nyström extension with QR decomposition

Both Algorithm1 and Algorithm2 contain steps that require us to compute the $m$ leading eigenvalues and corresponding eigenvectors of $L_{\textrm{mix}}$ . In the examples we consider in Section 7, the sizes of the graphs go up to tens of thousands of nodes, making it time-consuming to perform operations on the matrix.

The Nyström approximation [Reference Nyström58], which generates a low-rank approximation of the original matrix from a subset of its columns, is an effective way to tackle this issue. The choice of sampling method could affect the Nyström approximation performance, since different samples provide different approximations of the original adjacency matrix $W$ .

Before applying the Nyström method to our particular choices for the matrix $L_{\textrm{mix}}$ , we first give a brief explanation of the method for a general, symmetric matrix $C\in \mathbb{R}^{|V|\times |V|}$ .

We sample $k$ distinct pointsFootnote ²⁹ uniformly at random from $|V|$ points and partition the matrix $C$ as

(50) \begin{equation} C = \begin{pmatrix} C_{11}\;\;\;\; & C_{21}^T\\[5pt] C_{21}\;\;\;\; & C_{22} \end{pmatrix}, \end{equation}

where $C_{11} \in \mathbb{R}^{k \times k}$ , $C_{21} \in \mathbb{R}^{(|V|-k) \times k}$ , and $C_{22} \in \mathbb{R}^{(|V|-k) \times (|V|-k)}$ . We have relabelled the points such that our $k$ sampled points are the first $k$ points. We note that $C_{11}$ and $C_{22}$ are symmetric.

Because $C_{11}$ is a real and symmetric matrix, we can perform an eigenvalue decomposition to obtain $C_{11} = U \Lambda _k U^T$ , whereFootnote ³⁰ $\Lambda _k\;:\!=\;\textrm{diag}(\lambda _1, \ldots, \lambda _k)\in \mathbb{R}^{k\times k}$ is a diagonal matrix with the $k$ eigenvalues (counted according to multiplicity) $\lambda _i$ of $C_{11}$ on its diagonal and $U \in \mathbb{R}^{k\times k}$ is an orthogonal matrix which has the corresponding eigenvectors of $C_{11}$ as columns (in the order corresponding to the order of the eigenvalues in $\Lambda _k$ ).

We write $C_{11}^\dagger$ for the Moore–Penrose pseudoinverseFootnote ³¹ of $C_{11}$ . If $C_{11}$ is invertible, then the pseudoinverse is equal to the inverse of $C_{11}$ . If $C_{11}$ is not invertible, then the relationships $C_{11}^\dagger C_{11} C_{11}^\dagger = C_{11}^\dagger$ and $C_{11} C_{11}^\dagger C_{11} = C_{11}$ still hold. The pseudoinverse can be computed as $C_{11}^\dagger = U \Lambda _k^\dagger U^T$ , where $\Lambda _k^\dagger$ is the pseudoinverse. Specifically, this means that $\Lambda _k^\dagger$ is a diagonal matrix whose diagonal elements are the reciprocals of the diagonal elements of $\Lambda _k$ when those are non-zero, and zero when those are zero.

We would like to use the columns of $U$ to approximate $k$ eigenvectors of $C$ . Eigenvectors of $C$ are in $\mathbb{R}^{|V|}$ , whereas the columns of $U$ are in $\mathbb{R}^k$ ; the main idea of the Nyström extension is to approximate the ‘missing’ $|V|-k$ entries by $C_{21}U\Lambda _k^\dagger$ . This approximation is inspired by a quadrature rule; for more details, we refer to Fowlkes et al. [Reference Fowlkes, Belongie, Fan and Malik23] and Bertozzi and Flenner [Reference Bertozzi and Flenner4]. Thus,

\begin{align*} U_{C} \;:\!=\; \begin{pmatrix} U \\[5pt] C_{21} U \Lambda _k^\dagger \end{pmatrix} \in \mathbb{R}^{|V|\times k} \end{align*}

is an approximation of $k$ eigenvectors of $C$ , which in turn gives us an approximation of the full matrix $C$ :

(51) \begin{align} \bar{C} \;:\!=\; U_{C} \Lambda _k U_{C}^T &= \begin{pmatrix} U \\[5pt] C_{21} U \Lambda _k^\dagger \end{pmatrix} \Lambda _k \begin{pmatrix} U \\[5pt] C_{21} U \Lambda _k^\dagger \end{pmatrix}^T \nonumber \\[5pt] &= \begin{pmatrix} U\Lambda _k U^T & U \Lambda _k \Lambda _k^\dagger U^T C_{21}^T\\[5pt] C_{21} U \Lambda _k^\dagger \Lambda _K U^T & C_{21} U \Lambda _k^\dagger \Lambda _k \Lambda _k^\dagger U^T C_{21}^T\\[5pt] \end{pmatrix} = \begin{pmatrix} C_{11} & C_{11} C_{11}^\dagger C_{21}^T \\[5pt] C_{21} C_{11}^\dagger C_{11} & C_{21} C_{11}^\dagger C_{21}^T \end{pmatrix}\nonumber \\[5pt] &= \begin{pmatrix} C_{11}\\[5pt] C_{21} \end{pmatrix} C_{11}^\dagger \begin{pmatrix} C_{11} & C_{21}^T\end{pmatrix}. \end{align}

Comparing (50) and (51), we obtain that $C_{22}$ is approximated as

\begin{align*} C_{22} \approx C_{21} C_{11}^\dagger C_{21}^T. \end{align*}

We note that in many sources, the off-diagonal blocks are often given as $C_{21}$ and $C_{21}^T$ , rather than $C_{21} C_{11}^\dagger C_{11}$ and $C_{11} C_{11}^\dagger C_{21}^T$ , respectively. If $C_{11}$ is invertible, these are equivalent, but in general they need not be.Footnote ³² In our numerical tests in Section 7, the matrices that are used for $C_{11}$ are always invertible.

Now we wish to find the eigenvalue decomposition of $\bar{C}$ .Footnote ³³ We follow the QR decomposition method introduced in Budd et al. [Reference Budd, van Gennip and Latz12], which uses the thin QR decomposition [Reference Golub and Van Loan27, Theorems 5.2.2 and 5.2.3]:

\begin{equation*} \begin {pmatrix} C_{11} \\[5pt] C_{21}\end {pmatrix} = \tilde Q R, \end{equation*}

whereFootnote ³⁴ $\tilde Q\in \mathbb{R}^{|V|\times k}$ has orthonormal columns (i.e., $\tilde Q^T\tilde Q=I$ ) and $R\in \mathbb{R}^{k\times k}$ is an upper triangular matrix with positive diagonal entries. This decomposition is possible if $\begin{pmatrix} C_{11} \\ C_{21}\end{pmatrix}$ has full column rank, that is, if its column rank is $k$ . We note that this is guaranteed to be the caseFootnote ³⁵ if $C_{11}$ is invertible, and thus in particular in our numerical studies in Section 7.

By (51), we see that $\bar{C} = QR C_{11}^\dagger (QR)^T$ . The matrix $RC_{11}^\dagger R^T$ is real and symmetric and thus admits an eigendecomposition:

\begin{equation*} RC_{11}^\dagger R^T = \Upsilon \Sigma \Upsilon ^T, \end{equation*}

with $\Upsilon \in \mathbb{R}^{k\times k}$ orthogonal and $\Sigma \in \mathbb{R}^{k\times k}$ diagonal. Then

(52) \begin{equation} \bar{C} = \tilde QR C_{11}^\dagger R^T \tilde Q^T = \tilde Q \Upsilon \Sigma \Upsilon ^T \tilde Q^T = \tilde Q \Upsilon \Sigma (\tilde Q \Upsilon )^T. \end{equation}

Since $(\tilde Q \Upsilon )^T \tilde Q \Upsilon = I$ , we have that $\tilde Q\Upsilon \in \mathbb{R}^{|V|\times k }$ has orthonormal columns. Thus, we can view the $k$ diagonal entries of $\Sigma$ as approximate eigenvalues of $\bar{C}$ (and thus of $C$ ) with the columns of $\tilde Q\Upsilon$ the corresponding approximate eigenvectors.

Remark 6.2. The method explained above can be used to estimate $k$ eigenvalues and corresponding eigenvectors of the matrix $C$ based on sampling the submatrices $C_{11}$ and $C_{21}$ . In our applications, we require eigenvalues and eigenvectors of the matrix $D_C+C$ or of a normalised matrix $D_C^{-\frac 12} C D^{-\frac 12}$ or $D_C^{-1} C$ , where $D_C = \textrm{diag}(d_C)$ is the diagonal degree matrix based on $C$ with diagonal entries $(d_C)_i$ as in (1). Because we cannot compute $D_C$ exactly if we do not have access to all entries of $C$ , we first have to compute an approximation to $D_C$ . We set $\bar{D}_C\;:\!=\;\textrm{diag}(\bar{d}_C)$ with

\begin{align*} \bar{d}_C &\;:\!=\; \bar C \textbf{1}_{|V|} = \bar C \begin{pmatrix} \textbf{1}_k \\[5pt] \textbf{1}_{|V|-k} \end{pmatrix} = \begin{pmatrix} C_{11} \textbf{1}_k + C_{11} C_{11}^\dagger C_{21}^T \textbf{1}_{|V|-k} \\[5pt] C_{21} C_{11}^\dagger C_{11} \textbf{1}_k + C_{21} C_{11}^\dagger C_{21}^T \textbf{1}_{|V|-k} \end{pmatrix}, \end{align*}

where $\textbf{1}_k \in \mathbb{R}^k$ is the $k$ -dimensional column vector whose entries are all $\boldsymbol{1}$ (see [Reference Garcia-Cardona, Merkurjev, Bertozzi, Flenner and Percus24]).Footnote ³⁶

Now we approximate $D_C+C$ by $\bar{D}_C+\bar{C}$ and $D_C^{-1}C$ by $\bar{D}_C^\dagger \bar{C}$ . If $\bar{D}_C$ has non-zero diagonal elements, then $\bar{D}_C^\dagger =\bar{D}_C^{-1}$ . Moreover, we approximate $D_C^{-\frac 12} C D^{-\frac 12}$ by $(\bar{D}_C^\dagger )^{\frac 12} \bar{C} (\bar{D}_C^\dagger )^{\frac 12}$ , where we require $\bar{D}_C^\dagger$ (and thus also $\bar{D}_C$ ) to have non-negative diagonal elements, for the square root to be well defined.

In general, there is no guarantee that $\bar{D}_C$ has non-zero or non-negative diagonal elements, but all matrices of type $\bar{D}_C$ that we use in our numerical tests in Section 7 do. The following observation may be of use establishing these properties in specific situations. If we denote by $d_{C_{1:k}} \;:\!=\; C_{11} \textbf{1}_k + C_{21}^T \textbf{1}_{|V|-k}\in \mathbb{R}^k$ , the column vector containing the first $k$ entries of the actual degree vector $d_C\in \mathbb{R}^{|V|}$ , then

\begin{equation*} \bar {d}_C = \begin {pmatrix} C_{11} C_{11}^\dagger d_{C_{1:k}}\\[5pt] C_{21} C_{11}^\dagger d_{C_{1:k}} \end {pmatrix}. \end{equation*}

If $C$ is the adjacency matrix of a graph with positive degrees for all nodes (e.g., a connected graph), then all entries of $\bar{d}_C$ will be positive if the matrices $C_{11}C_{11}^\dagger$ and $C_{21} C_{11}^\dagger$ both preserve entrywise positivity of vectors. If $C_{11}$ is invertible, this is clearly the case for $C_{11}C_{11}^\dagger$ . In fact, in that case the first $k$ entries of $\bar{d}_C$ are exactly equal to the first $k$ entries of the actual degree vector $d_C$ .

To obtain approximate eigenvalues and eigenvectors of $D_C+C$ , $D_C^{-\frac 12} C D^{-\frac 12}$ , or $D_C^{-1} C$ , we can now apply the decomposition from (52) to $D_C+C$ , $(\bar{D}_C^\dagger )^{\frac 12} \bar{C} (\bar{D}_C^\dagger )^{\frac 12}$ or $\bar{D}_C^\dagger \bar{C}$ , respectively, instead of to $\bar{C}$ .