Measuring the stability of spectral clustering

TL;DR

This paper introduces new measures for the stability of spectral clustering based on the minimal perturbation needed to cause ambiguity in the spectral gap, and proposes algorithms to compute these measures.

Contribution

It defines structured and unstructured distances to ambiguity for spectral clustering stability and develops a two-level iterative algorithm to compute these measures.

Findings

Structured and unstructured stability indicators can differ significantly.

Maximal stability criterion may suggest different cluster counts depending on the measure.

Numerical experiments demonstrate the effectiveness of the proposed stability measures.

Abstract

As an indicator of the stability of spectral clustering of an undirected weighted graph into $k$ clusters, the $k$th spectral gap of the graph Laplacian is often considered. The $k$th spectral gap is characterized in this paper as an unstructured distance to ambiguity, namely as the minimal distance of the Laplacian to arbitrary symmetric matrices with vanishing $k$th spectral gap. As a conceptually more appropriate measure of stability, the structured distance to ambiguity of the $k$-clustering is introduced as the minimal distance of the Laplacian to Laplacians of graphs with the same vertices and edges but with weights that are perturbed such that the $k$th spectral gap vanishes. To compute a solution to this matrix nearness problem, a two-level iterative algorithm is proposed that uses a constrained gradient system of matrix differential equations in the inner iteration and a one-dimensional optimization of the perturbation size in the outer iteration. The structured and unstructured distances to ambiguity are compared on some example graphs. The numerical experiments show, in particular, that selecting the number $k$ of clusters according to the criterion of maximal stability can lead to different results for the structured and unstructured stability indicators.

Tables

Table 1. Table 1: Frequency of k 𝑘 k being the optimal number of clusters provided by the spectral gaps and the SDA in Example  5.2 ( α = 1 / 2 𝛼 1 2 \alpha=1/2 ).

$k$	$4$	$5$	$6$	$7$	$8$
$g_{\mathrm{opt}}$	0.0%	0.0%	69.6%	20.4%	0.0%
${\delta }_{\mathrm{opt}}$	0.0%	0.0%	78.0%	5.2%	6.8%

Table 2. Table 2: Frequency of k 𝑘 k being the optimal number of clusters provided by the spectral gaps and the SDA in Example  5.2 ( α = 1 / 4 𝛼 1 4 \alpha=1/4 ).

$k$	$4$	$5$	$6$	$7$	$8$
$g_{\mathrm{opt}}$	0.0%	0.0%	89.6%	7.6%	2.8%
${\delta }_{\mathrm{opt}}$	0.0%	0.0%	95.6%	3.2%	1.2%

Equations

$$W=(w_{ij})\in {\mathbb{R}}^{n\times n}.$$

$$D=\mathrm{diag}(d_i),d_i=(W1{)}_i,where1:=(1,\dots ,1{)}^T\in {\mathbb{R}}^n.$$

$$L=D-W,\text{ i.e.,}L(W)=\mathrm{diag}(W1)-W.$$

$$C_i\cap C_j=\varnothing \text{mbox}fori,j=1,\dots ,k,i\ne j\text{mbox}and\bigcup _{i=1}^k{C}_i=\mathcal{V}.$$

$$L_{\mathrm{Sym}}:=D^{-1/2}L{D}^{-1/2}=I-D^{-1/2}W{D}^{-1/2}.$$

$$\frac{{\lambda }_{k+1}-{\lambda }_k}{\sqrt{2}}=\min \left\{\Vert L(W)-\hat{L}{\Vert }_F:\hat{L}\in \mathrm{Sym}(n),{\lambda }_k(\hat{L})={\lambda }_{k+1}(\hat{L})\right\},$$

$$\Vert L-\tilde{L}{\Vert }_F^2$$

$$\displaystyle=2\Bigl{(}\frac{\lambda_{k+1}-\lambda_{k}}{2}\Bigr{)}^{2}=\tfrac{1}{2}(\lambda_{k+1}-\lambda_{k})^{2}$$

$$\frac{{\lambda }_{k+1}-{\lambda }_k}{\sqrt{2}}=\Vert L-\tilde{L}{\Vert }_F.$$

$$\tfrac{1}{2}(\lambda_{k+1}-\lambda_{k})^{2}=\min_{\lambda}\Bigl{(}(\lambda_{k}-\lambda)^{2}+(\lambda_{k+1}-\lambda)^{2}\Bigr{)}\leq(\lambda_{k}-\widehat{\lambda}_{k})^{2}+(\lambda_{k+1}-\widehat{\lambda}_{k+1})^{2}.$$

$$\sum _i({\lambda }_i-{\hat{\lambda }}_i{)}^2\le \Vert L-\hat{L}{\Vert }_F^2,$$

$$\frac{1}{2}({\lambda }_{k+1}-{\lambda }_k{)}^2\le \Vert L-\hat{L}{\Vert }_F^2,$$

$$\tilde{L}=L(W)+\frac{1}{2}({\lambda }_{k+1}-{\lambda }_k)x_k{x}_k^T-\frac{1}{2}({\lambda }_{k+1}-{\lambda }_k)x_{k+1}{x}_{k+1}^T,$$

$$\displaystyle\begin{aligned} \delta_{k}(W)&=\min\big{\{}&&\|L(W)-L(\widehat{W})\|_{F}:\,\widehat{W}\in\mathbb{R}^{n\times n},\,\lambda_{k}(L(\widehat{W}))=\lambda_{k+1}(L(\widehat{W}))\,,\\ &&&\widehat{w}_{ij}=\widehat{w}_{ji}\geq 0\text{ for }(i,j)\in\mathcal{E}\text{ and }\widehat{w}_{ij}=0\text{ for }(i,j)\notin\mathcal{E}\big{\}}\,.\end{aligned}$$

$${\delta }_k(W)\ge ({\lambda }_{k+1}-{\lambda }_k)/\sqrt{2},$$

$$k_{\mathrm{opt}}(W)=\arg \underset{k}{\max }{\delta }_k(W),$$

$$F_{\varepsilon}(E)=\lambda_{k+1}\bigl{(}L(W+\varepsilon E)\bigr{)}-\lambda_{k}\bigl{(}L(W+\varepsilon E)\bigr{)},$$

$${\delta }_k(W)={\epsilon }^{\star }=\Vert L(W)-L(W^{\star }){\Vert }_F.$$

$$P_{\mathcal{E}}(A)\big{|}_{ij}:=\begin{cases}\displaystyle{\tfrac{1}{2}}\left(a_{ij}+a_{ji}\right)\,&\text{if }(i,j)\in\mathcal{E}\\[2.84526pt] 0\,&\text{otherwise}.\end{cases}$$

$$⟨P_{\mathcal{E}}(A),W⟩=⟨A,W⟩\text{ for all }W\in \mathrm{Sym}(\mathcal{E}).$$

$$L:\mathrm{Sym}(\mathcal{E})\to {\mathbb{R}}^{n\times n}.$$

$$L^{\ast }:{\mathbb{R}}^{n\times n}\to \mathrm{Sym}(\mathcal{E}),$$

$$L^{\ast }(V)=P_{\mathcal{E}}(\mathrm{diagvec}(V)1^T-V),$$

$$L^{\ast }(L(W))=P_{\mathcal{E}}(d{1}^T)+W\text{with }d=W1.$$

$$⟨V,L(W)⟩$$

$$=⟨\mathrm{diagvec}(V)1^T,W⟩-⟨V,W⟩=⟨P_{\mathcal{E}}(\mathrm{diagvec}(V)1^T-V),W⟩,$$

$$E\in \mathrm{Sym}(\mathcal{E}),\Vert L(E){\Vert }_F=1,W+\epsilon E\ge 0.$$

$$\dot{\lambda }=⟨L^{\ast }(x{x}^T),\epsilon \dot{E}⟩.$$

$$\frac{1}{\epsilon }\frac{d}{dt}{F}_{\epsilon }(E(t))=⟨G_{\epsilon }(E(t)),\dot{E}(t)⟩$$

$$G_{\epsilon }(E)$$

$$=P_{\mathcal{E}}\left((x_{k+1}\bullet x_{k+1}-x_k\bullet x_k)1^T-\left(x_{k+1}{x}_{k+1}^T-x_k{x}_k^T\right)\right),$$

$$L^{\ast }(x{x}^T-y{y}^T)\ne 0.$$

Full text

AMS classification: 15A18, 65K05, 15A60

1 Introduction

Clustering — identifying groups with similarities in a large data set — is a fundamental task in data analysis. In order to cluster an undirected weighted graph, several algorithms are considered in the literature. Spectral clustering (originating with Fiedler [5]; see [19] and [21] for more recent accounts) offers significant advantages over more traditional combinatorial techniques such as direct $k$-means or single linkage. Quoting [21], “Results obtained by spectral clustering often outperform the traditional approaches; spectral clustering is very simple to implement and can be solved efficiently by linear algebra methods.”

The stability of spectral clustering algorithms is heuristically — and very conveniently — associated with the spectral gap in the graph Laplacian, that is the difference between the $k$th and $(k+1)$-st smallest eigenvalues of the Laplacian matrix in the case of clusters built from $k$ eigenvectors. If these two eigenvalues coalesce, then the spectral $k$-clustering becomes ambiguous, and arbitrarily small perturbations may yield different clusterings.

In this paper we measure the stability of spectral $k$-clustering, for any number $k$ of clusters, by the structured distance to ambiguity, which is introduced here as the distance (w.r.t. the Frobenius norm) of the given Laplacian matrix to the set of Laplacian matrices that correspond to a weighted graph with the same vertices and edges but perturbed weights such that the $k$th and $(k+1)$-st smallest eigenvalues of the graph Laplacian coalesce — and for which therefore spectral clustering with $k$ clusters becomes ambiguous. On the other hand, the $k$th spectral gap divided by $\sqrt{2}$ will be characterized as the distance of the Laplacian matrix to the much larger set of arbitrary, unstructured symmetric matrices with coalescing $k$th and $(k+1)$-st eigenvalues.

A stability indicator is useful in selecting an appropriate number $k$ of clusters, by choosing $k$ such that the stability indicator is maximized. The structured distance to ambiguity, whose computation will be discussed in this paper, may be too costly computationally to be used as a general-purpose stability indicator in comparison with the spectral gap, which is available without further computations in spectral clustering. Nevertheless, comparing the structured distance to ambiguity and the spectral gap for representative examples within a class of graphs of interest gives insight into the reliability (or otherwise) of the spectral gap as a stability indicator.

The structured distance to ambiguity considered in this paper is similar in spirit to other structured robustness measures that arise in matrix and control theory, such as stability radii, passivity distances, distance to singularity, etc., see [13, 11, 12, 4, 7, 10, 17]. The notion of stability considered here is, however, different from [22], where a statistical perspective on clustering stability is developed.

A main objective of this paper is to show how the structured distance to ambiguity can be computed. The proposed algorithm is an iterative algorithm, where in each step a pair of adjacent eigenvalues and associated eigenvectors of the Laplacian of a graph with perturbed weights are computed. For a large sparse graph (where the number of edges leaving any vertex is moderately bounded), these computations can typically be done with a complexity that is linear in the number of vertices.

A feature of this algorithm in common with recent algorithms for eigenvalue optimization as given in [8, 7, 10, 9, 3] is a two-level procedure for matrix nearness problems, where in an inner iteration a gradient flow drives perturbations to the original matrix of a fixed size into a (local) minimum of a nonnegative functional that depends on eigenvalues and eigenvectors, and in an outer iteration the perturbation size is optimized such that the functional becomes zero.

The paper is organized as follows. In Section 2 we recall basic concepts of spectral clustering. In Section 3 we introduce the structured distance to ambiguity (SDA) as a measure of stability for spectral $k$-clustering and we characterize the $k$th spectral gap as an unstructured distance to ambiguity. In Section 4 we describe a numerical method to compute the SDA, which requires the solution of a structured matrix nearness problem. We propose a two-level iterative method, which is based on a gradient system of matrix differential equations in an inner iteration and a one-dimensional optimization of the perturbation size in an outer iteration, and we discuss algorithmic aspects. In Section 5 we present the results of numerical experiments where the spectral gap (i.e., the unstructured distance to ambiguity) and the theoretically more appropriate measure of stability discussed here (i.e., the structured distance to ambiguity) are compared on some classes of graphs. We give a conclusion in Section 6.

2 Spectral clustering

Let $\mathcal{G}=(\mathcal{V},\mathcal{E})$ be an undirected graph with vertex set $\mathcal{V}=\{1,\dots,n\}$ and edge set $\mathcal{E}\subset\mathcal{V}\times\mathcal{V}$. We assume that the graph is weighted, that is, a non-negative weight $w_{ij}=w_{ji}\geq 0$ is associated with each edge $(i,j)\in\mathcal{E}$ between two vertices $i$ and $j$. We set $w_{ij}=0$ for $(i,j)\notin\mathcal{E}$. The weighted adjacency matrix of the graph is the matrix

[TABLE]

The degrees $d_{i}=\sum_{j=1}^{n}w_{ij}$ are the elements of the diagonal matrix

[TABLE]

The (unnormalized) Laplacian of the graph is given by the matrix $L=L(W)$,

[TABLE]

The Laplacian matrix $L$ is symmetric and positive semi-definite; since - by construction - $L\mathbf{1}=0$, $\lambda_{1}=0$ is the smallest eigenvalue of $L$. Note that the matrix $L$ does not depend on (possible) diagonal elements of the matrix $W$, which means that self-edges do not change the graph Laplacian. The graph Laplacian and its eigenvectors provide important instruments to spectral clustering, as stated by the following theorems.

Theorem 2.1** (Bi-partition; Fiedler [5])**

Let $W\in\mathbb{R}^{n\times n}$ be the weight matrix of an undirected graph and $L$ the corresponding Laplacian matrix. Let $0=\lambda_{1}\leq\lambda_{2}\leq\ldots\leq\lambda_{n}$ be the eigenvalues of $L$. Then, the graph is disconnected if and only if $\lambda_{2}=0$. Moreover, if $0=\lambda_{2}<\lambda_{3}$, then the entries of the corresponding eigenvector orthogonal to $\mathbf{1}$ assume only two different values, of different sign, which mark the membership to the two connected components.

If $\lambda_{2}$ is a simple eigenvalue, then the corresponding eigenvector is known as the Fiedler vector. In spectral graph theory, inspired by Theorem 2.1, it is common to compute the second smallest eigenvalue of $L$ and label the positive components of the Fiedler vector as belonging to one subset and the negative ones to another subset, and in this way obtaining a natural partition of the graph. However, this becomes unreliable when a small perturbation of the weights yields a coalescence of the eigenvalues $\lambda_{2}$ and $\lambda_{3}$.

More generally we have the following result (see e.g. [21]). For a subset of vertices $C\subset\mathcal{V}$, we here denote the indicator vector $\mathbf{1}_{C}$ as the vector whose $\ell$th entry is equal to $1$ if $v_{\ell}\in C$ and is equal to zero otherwise.

Theorem 2.2** ($k$-partition)**

Let $W\in\mathbb{R}^{n\times n}$ be the weight matrix of an undirected graph and $L$ the corresponding Laplacian matrix. Then the multiplicity $k$ of the eigenvalue [math] (the dimension of $\ker(L)$) equals the number of connected components $C_{1},\ldots,C_{k}$ in the graph. The eigenspace of the eigenvalue [math] is spanned by the indicator vectors $\mathbf{1}_{C_{1}},\ldots,\mathbf{1}_{C_{k}}$.

Nonempty sets $C_{1},\ldots,C_{k}$ form a clustering of the graph if

[TABLE]

Similar to the case of two clusters, this result motivates an algorithm for clustering a graph into $k$ components, which is reviewed in Algorithm 1; see [21]. For the classical $k$-means algorithm we refer the reader e.g. to [16]. Analogous results and algorithms are extended to the normalized Laplacian

[TABLE]

Algorithm 1 Unnormalized spectral clustering algorithm

1:weight matrix $W$, number $k$ of clusters

2:Clusters $C_{1},\ldots,C_{k}$

3:Begin

4:     Compute the Laplacian matrix $L=D-W$

5:     Compute the eigenvectors $x_{1},\ldots,x_{k}$ associated with the $k$ smallest nonzero eigenvalues of $L$

6:     Set $X=\left[x_{1}\ |\ x_{2}\ |\ \ldots\ |\ x_{k}\right]$

7:     For $i=1,\ldots,n$, define $r_{i}\in\mathbb{R}^{k}$ the vector given by the $i$-th row of $X$

8:     Cluster the points $\left(r_{i}\right)_{i=1,\ldots,n}$ by the $k$-means algorithm into $k$ clusters $C_{1},\ldots,C_{k}$

9:Return $C_{1},\ldots,C_{k}$

3 Structured distance to ambiguity as a stability measure for spectral clustering

In order to evaluate the robustness of the spectral clustering algorithm it is essential to quantify the sensitivity to perturbations of the invariant subspace associated with the $k$ eigenvectors used by the algorithm. Suppose that a small perturbation to the weights of the graph makes the eigenvalues $\lambda_{k}$ and $\lambda_{k+1}$ of the Laplacian coalesce. Then the eigenvector associated with the computed $k$th smallest eigenvalue can change completely and hence can yield a different clustering.

It has been suggested in the literature to use the spectral gap $\lambda_{k+1}-\lambda_{k}$ as an indicator of the stability of the $k$-clustering; cf. [21]. For unstructured perturbations of the graph Laplacian, this is usually motivated by the Davis–Kahan theorem (see, e.g., [20] and [6]), which tells us that the distance between the eigenspaces of the Laplacian of a graph and of any perturbed symmetric matrix has a bound that is proportional to the perturbation size and inversely proportional to the spectral gap $\lambda_{k+1}-\lambda_{k}$ of the Laplacian $L(W)$.

In another direction, which is more related to the concepts studied here, the spectral gap can be characterized as an unstructured distance to ambiguity (up to the scaling factor $\sqrt{2}$) by the following result, where $\|A\|_{F}=\sqrt{\sum_{i,j}a_{ij}^{2}}$ denotes the Frobenius norm of a matrix $A=(a_{ij})$.

Theorem 3.1

The $k$th spectral gap, divided by $\sqrt{2}$, is characterized as

[TABLE]

where $\mathrm{Sym}(n)$ denotes the set of all symmetric real $n\times n$-matrices.

Proof. Let here $L=L(W)$ for brevity.

(a) We diagonalize $L=Q\Lambda Q^{T}$ with the diagonal matrix $\Lambda=\text{diag}(\lambda_{i})$ of eigenvalues and the orthogonal matrix $Q$ formed of the eigenvectors. If we set $\widetilde{L}=Q\widetilde{\Lambda}Q^{T}$ with $\widetilde{\Lambda}=\text{diag}(\widetilde{\lambda}_{i})$ and with $\widetilde{\lambda}_{i}=\lambda_{i}$ for $i\notin\{k,k+1\}$ and coalescing $k$th and $(k+1)$-st eigenvalues $\widetilde{\lambda}_{k}=\widetilde{\lambda}_{k+1}=\tfrac{1}{2}(\lambda_{k}+\lambda_{k+1})$, then we obtain

[TABLE]

so that

[TABLE]

(b) Let now $\widehat{L}$ be an arbitrary symmetric matrix with ordered eigenvalues $\widehat{\lambda}_{i}$ such that the $k$th and $(k+1)$-st eigenvalues coalesce: $\widehat{\lambda}_{k}=\widehat{\lambda}_{k+1}$. We then have

[TABLE]

By the Hoffman–Wielandt theorem [15, p. 368],

[TABLE]

and so we find

[TABLE]

which concludes the proof. ∎

The interpretation is that applying Algorithm 1 is justified when the gap is relatively large, otherwise a small perturbation of the weight matrix may yield coalescence of $\lambda_{k}$ and $\lambda_{k+1}$ and significantly change the clustering.

The conceptual weakness of this approach is that it considers unstructured perturbations of the Laplacian, as opposed to the admissible perturbations that are Laplacians $L(\widehat{W})$ of a perturbed weight matrix that preserves the symmetry and non-negativity, i.e., $\widehat{w}_{ij}=\widehat{w}_{ji}\geq 0$, and the sparsity pattern of $W$, i.e., $\widehat{w}_{ij}=0$ for $(i,j)\notin\mathcal{E}$. In (1), the minimizer is given by

[TABLE]

where $x_{k}$ and $x_{k+1}$ denote normalized eigenvectors of $L(W)$ to the eigenvalues $\lambda_{k}$ and $\lambda_{k+1}$, respectively. Apart from exceptional cases, $\widetilde{L}$ is not the Laplacian of a graph.

We therefore propose the following stability measure for $k$-clustering, which is given by a structured distance to ambiguity,

[TABLE]

Here, the Laplacian $L(W)$ is not compared with arbitrary symmetric matrices with vanishing $k$th spectral gap as in (1), but with the Laplacians of weight matrices on the same graph with vanishing $k$th spectral gap. In view of (1), it is clear that

[TABLE]

and often $\delta_{k}(W)$ is substantially larger than the scaled spectral gap. In view of Theorem 2.1, we further note that $\delta_{1}(W)$ is the Frobenius-norm distance of the Laplacian $L(W)$ to that of a nearest disconnected graph.

For clustering a graph it is usually not known beforehand what the best choice of the number $k$ of clusters is. Asking for the most stable clustering of a graph (most stable with respect to admissible perturbations of the weights) determines the optimal number $k$ as

[TABLE]

where $k$ can be limited to $k\leq k_{\max}$ and $k_{\max}$ is given a priori or chosen such that $\lambda_{k_{\max}}$ is smaller than some threshold. This criterion for the selection of $k$ is considered here instead of choosing $k$ such that the spectral gap $\lambda_{k+1}-\lambda_{k}$ is maximized. The latter is computationally cheaper, but appears conceptually less meaningful. In our numerical experiments in Section 5 we will compare the two criteria for some families of graphs.

Remark 3.1

There are obvious alternatives to the above stability measure:

• •

Instead of the unnormalized Laplacian $L(W)=D-W$ with $D=\mathrm{diag}(W\mathbf{1})$ we might work with the normalized Laplacian $I-D^{-1/2}WD^{-1/2}$.

• •

Instead of minimizing the perturbation of the Laplacian we might minimize the perturbation of the weights, $\|W-\widehat{W}\|_{F}$.

• •

Instead of allowing for perturbations in every entry of the weight matrix $W$, perturbations might be restricted to selected entries.

In the following we concentrate on $\delta_{k}(W)$ as given above, but our algorithm for computing $\delta_{k}(W)$ is readily extended to the other cases.

4 Computation of the structured distance to ambiguity

In this section we describe an approach to compute the structured distance to ambiguity $\delta_{k}(W)$ defined in the previous section.

4.1 Outline of the computational approach

Our approach is summarized by the following two-level method:

• •

Inner iteration: Given $\varepsilon>0$, we aim to compute a symmetric matrix $E=(e_{ij})\in\mathbb{R}^{n\times n}$ with the same sparsity pattern as $W$ (i.e., $e_{ij}=0$ if $w_{ij}=0$), with $L(E)=\mathrm{diag}(E\mathbf{1})-E$ of unit Frobenius norm, with $W+\varepsilon E\geq 0$ (with entrywise inequality) such that the difference between the $(k+1)$-st and $k$th smallest eigenvalues of $L(W+\varepsilon E)$ is minimal. The obtained minimizer is denoted by $E_{\varepsilon}$.

• •

Outer iteration: We compute the smallest value of $\varepsilon$ such that the $k$th and $(k+1)$-st eigenvalues of $L(W+\varepsilon E_{\varepsilon})$ coalesce.

In order to compute $E_{\varepsilon}$ for a given $\varepsilon>0$, we make use of a constrained gradient system for the functional

[TABLE]

under the constraints of unit Frobenius norm of $L(E)=\mathrm{diag}(E\mathbf{1})-E$, the non-negativity $W+\varepsilon E\geq 0$ and symmetry and the sparsity pattern of $E$.

In the outer iteration we compute the smallest $\varepsilon$ with $F_{\varepsilon}(E_{\varepsilon})=0$, denoted $\varepsilon^{\star}$, by a combined Newton-bisection method. The algorithm computes the perturbed weight matrix $W^{\star}=W+\varepsilon^{\star}E_{\varepsilon^{\star}}$ whose Laplacian has coalescent $k$th and $(k+1)$-st eigenvalues. We then have

[TABLE]

Remark 4.1

The constrained minimization of $F_{\varepsilon}$ for a given $\varepsilon$ is a nonconvex and nonsmooth optimization problem, which may have several local minima. By following a trajectory of the constrained gradient system into a stationary point, as we do here, we may thus end up in a local (instead of global) minimizer $E_{\varepsilon}$. In our numerical experiments this situation was observed only in contrived examples [3, Example 2]. Using several trajectories starting from different random initial values can mitigate this potential problem. In any case, the algorithm computes an $\varepsilon^{\star}=\|L(W)-L(W^{\star})\|_{F}$ that is an upper bound to $\delta_{k}(W)$, and the bound can be made tighter by running more trajectories.

4.2 Gradient of the functional $F_{\varepsilon}$

We denote by $\langle X,Y\rangle=\mathrm{trace}(X^{T}Y)$ the inner product of matrices that corresponds to the Frobenius norm $\|X\|_{F}=\langle X,X\rangle^{1/2}$ on $\mathbb{R}^{n\times n}$.

For the edge set $\mathcal{E}$ and for an arbitrary matrix $A=(a_{ij})\in\mathbb{R}^{n\times n}$, we define the symmetric projection onto the sparsity pattern given by $\mathcal{E}$ as

[TABLE]

We denote by $\mathrm{Sym}(\mathcal{E})$ the space of symmetric matrices with sparsity pattern $\mathcal{E}$ and note that $P_{\mathcal{E}}$ is the orthogonal projection from $\mathbb{R}^{n\times n}$ to $\mathrm{Sym}(\mathcal{E})$:

[TABLE]

The Laplacian operator is a linear map

[TABLE]

Its adjoint with respect to the Frobenius inner product,

[TABLE]

is given in the following lemma.

Lemma 4.1

For $V\in\mathbb{R}^{n\times n}$, let $L^{*}(V)\in\mathrm{Sym}(\mathcal{E})$ be defined by $\langle L^{*}(V),W\rangle=\langle V,L(W)\rangle$ for all $W\in\mathrm{Sym}(\mathcal{E})$. Then,

[TABLE]

where $\mathrm{diagvec}(V)\in\mathbb{R}^{n}$ is the vector of the diagonal entries of $V$. Furthermore, if $\mathcal{E}$ contains no auto-loops, i.e., $(i,i)\notin\mathcal{E}$ for all $i=1,\dots,n$, then we have for $W\in\mathrm{Sym}(\mathcal{E})$,

[TABLE]

Proof. We have, for all $V\in\mathbb{R}^{n\times n}$ and $W\in\mathrm{Sym}(\mathcal{E})$,

[TABLE]

which yields the formula for $L^{*}(V)$. The formula for $L^{*}(L(W))$ then follows on inserting $V=L(W)=\mathrm{diag}(W\mathbf{1})-W$ into this formula.∎

The following lemma is basic for the derivation of the gradient system of matrix differential equations that we use to minimize the functional $F_{\varepsilon}$ of (2). Here we call a matrix $E\in\mathbb{R}^{n\times n}$ $\varepsilon$-feasible if

[TABLE]

Lemma 4.2

[3, Lemma 3.2]** Let $\varepsilon>0$ and consider a differentiable path of $\varepsilon$-feasible matrices $E(t)$ $(t_{0}<t<t_{1})$. Let $\lambda(t)$ be a continuous path of simple eigenvalues of the Laplacian matrix $L(W+\varepsilon E(t))$ with the associated eigenvector $x(t)$ of unit Euclidean norm. Then, $\lambda(\cdot)$ is differentiable and (denoting the derivative w.r.t. $t$ by a dot and omitting the argument $t$)

[TABLE]

We note that $L^{*}(xx^{T})=P_{\mathcal{E}}\left((x\bullet x)\mathbf{1}^{T}-xx^{T}\right)$, where $x\bullet x=(x_{i}^{2})\in\mathbb{R}^{n}$ denotes the vector of squares of the entries of the vector $x=(x_{i})\in\mathbb{R}^{n}$.

An immediate consequence of Lemma 4.2 applied to the eigenvalues in the functional $F_{\varepsilon}$ of (2) is the following result.

Theorem 4.1

Let $\varepsilon>0$. Along a differentiable path $E(t)$ of $\varepsilon$-feasible matrices such that the $k$th and $(k+1)$-st eigenvalues $\lambda_{k}(t)$ and $\lambda_{k+1}(t)$ of the Laplacian $L(W+\varepsilon E(t))$ are simple eigenvalues, we have for the functional $F_{\varepsilon}$ of (2) that

[TABLE]

with the (rescaled) gradient

[TABLE]

where $x_{k}$ and $x_{k+1}$ are the normalized eigenvectors associated with the eigenvalues $\lambda_{k}$ and $\lambda_{k+1}$, respectively, of $L(W+\varepsilon E)$.∎

The following result shows that the gradient $G_{\varepsilon}(E)$ cannot be zero. This fact will be important later.

Lemma 4.3

Let $x,y\in\mathbb{R}^{n}$ be normalized eigenvectors of a graph Laplacian $\widehat{L}=L(\widehat{W})$ corresponding to different eigenvalues $\lambda\neq\mu$. Then,

[TABLE]

Proof. We show that assuming $L^{*}(xx^{T}-yy^{T})=0$ leads to a contradiction. We would then have

[TABLE]

i.e.,

[TABLE]

Multiplying with the weights $w_{ij}$ and summing over all $(i,j)\in\mathcal{E}$ yields, however,

[TABLE]

and in the same way

[TABLE]

which is in contradiction to (5). ∎

4.3 Gradient flow under the unit-norm constraint

We consider now the constrained gradient system associated with the functional $F_{\varepsilon}$ with the constraint that $L(E(t))$ has unit Frobenius norm. We then get the following matrix differential equation:

[TABLE]

Here we note that the denominator cannot vanish, because $L^{*}(L(E))=0$ would imply that $0=\langle L^{*}(L(E)),E\rangle=\langle L(E),L(E)\rangle=\|L(E)\|_{F}^{2}$ in contradiction to $\|L(E)\|_{F}=1$. The role of $\kappa$ is that of a Lagrange multiplier that ensures that the constraint $\|L(E)\|_{F}=1$ is satisfied. The given formula for $\kappa$ is obtained as follows: Differentiating the constraint $\|L(E)\|_{F}^{2}=\langle L(E),L(E)\rangle=1$ gives $\langle\dot{E},L^{*}(L(E))\rangle=\langle L(\dot{E}),L(E)\rangle=\frac{1}{2}\frac{d}{dt}\|L(E)\|_{F}^{2}=0$. Taking the inner product of both sides of the differential equation with $L^{*}(L(E))$ yields

[TABLE]

which gives $\kappa$.

4.4 Non-negativity constraint

It may happen that along the solution trajectory of (6) some positive entry of $W+\varepsilon E$ becomes negative. In our experiments we never observed such a situation although this could potentially occur. In this subsection we explain how to deal with the further constraint $W+\varepsilon E(t)\geq 0$ for all $t$.

One possibility would be to follow the lines of [3] and consider KKT conditions by managing active constraints and integrating a piecewise smooth system of differential equations. Another possibility is to add a penalization term such as the following:

[TABLE]

where $(a)_{-}=\min\left(a,0\right)$, and to minimize the functional

[TABLE]

for increasing $c$. We have

[TABLE]

giving an extra term to the gradient system (6).

With the notation

[TABLE]

the differential equation (6) is then replaced by

[TABLE]

4.5 Monotonicity of the spectral gap and stationary points

The following monotonicity result follows directly from the construction of the gradient system.

Theorem 4.2

Let $E(t)$ with $L(E(t))$ of unit Frobenius norm satisfy the differential equation (8) with $G_{\varepsilon,c}(E)$ of (7). Then, $F_{\varepsilon,c}(E(t))$ decreases monotonically with $t$. Moreover, if $W+\varepsilon E(t)\geq 0$, then

[TABLE]

where $\lambda_{k}(t)$ and $\lambda_{k+1}(t)$ are the $k$th and $(k+1)$-st eigenvalues of $L(W+\varepsilon E(t))$.

Proof. Using (8) we obtain, with $G=G_{\varepsilon,c}(E)$,

[TABLE]

where the final inequality follows directly from the Cauchy–Schwarz inequality. This yields the monotonic decay of $F_{\varepsilon,c}(E(t))$ and hence also the second statement, since $F_{\varepsilon,c}(E(t))=\lambda_{k+1}(t)-\lambda_{k}(t)$ if $W+\varepsilon E(t)\geq 0$. ∎

Stationary points of (8) are characterized as follows.

Theorem 4.3

The following statements are equivalent along solutions of (8).

$\frac{d}{dt}F_{\varepsilon,c}\left(E(t)\right)=0$. 2. 2.

$\dot{E}=0$. 3. 3.

$G_{\varepsilon,c}(E)$* is a real multiple of $L^{*}(L(E))$.*

Proof. By the strict Cauchy–Schwarz inequality in (9), $\frac{d}{dt}F_{\varepsilon,c}\left(E(t)\right)$ can be zero only if $G$ is a multiple of $L^{*}(L(E))$, and then $\langle G,\dot{E}\rangle=0$ with $\dot{E}=-G+\kappa L^{*}(L(E))$ implies $\dot{E}=0$. ∎

The inner iteration is formulated as pseudocode in Algorithm 2, where a projected explicit Euler method is used for the discretization of the differential equation (8).

[h!]

Algorithm 2 Inner iteration

1:$W$, $k$, $\varepsilon$, $c$ (penalization parameter), $E_{0}$ (initial perturbation), tol (tolerance), $h$ (step size for gradient flow)

2:$E_{\varepsilon,c}$, $F_{\varepsilon,c}(E_{\varepsilon,c})$

3:Begin

4:     Compute $L(W+\varepsilon E_{0})$, $F_{\varepsilon,c}(E_{0})$

5:     Initialize $F^{(new)}_{\varepsilon,c}=F_{\varepsilon,c}(E_{0})$, $F^{(old)}_{\varepsilon,c}=F^{(new)}_{\varepsilon,c}+1$, $E=E_{0}$

6:     while $|F^{(new)}-F^{(old)}|/|F^{(new)}|>\mathrm{tol}$ do

7:         Compute $G_{\varepsilon,c}(E)$, $L^{*}(L(E))$ and $\kappa$

8:         Compute $\dot{E}=-G_{\varepsilon,c}(E)+\kappa L^{*}(L(E))$

9:         Update $E:=E+h\dot{E}$

10:         Normalize $E$ such that $\|L(E)\|_{F}=1$

11:         Compute $L(W+\varepsilon E)$ and $F_{\varepsilon,c}(E)$

12:         Update $F^{(old)}_{\varepsilon,c}=F^{(new)}_{\varepsilon,c}$, $F^{(new)}_{\varepsilon,c}=F_{\varepsilon,c}(E)$

13:     end while

14:Return $E_{\varepsilon,c}=E$ and $F_{\varepsilon,c}(E_{\varepsilon,c})=F_{\varepsilon,c}^{(new)}$

Remark 4.2

This algorithm can be refined by using an adaptive stepsize strategy as in [3, Section 5.3]. In particular, when the gradient flow is close to the minimum, there can be oscillations in the functional $F_{\varepsilon,c}(E)$, which are due to the discrete steps in the gradient system. This can be avoided by repeating lines 7 to 9 in Algorithm 2 with a smaller stepsize $h$ if the functional $F_{\varepsilon,c}(E)$ grows.

4.6 One-dimensional optimization over $\varepsilon$

Let $E_{\varepsilon,c}$ denote the minimizer of the functional $F_{\varepsilon,c}$, and let $c$ be so large that $W+\varepsilon E_{\varepsilon,c}\geq 0$, which implies that $F_{\varepsilon,c}(E_{\varepsilon,c})=F_{\varepsilon}(E_{\varepsilon,c})=:f(\varepsilon)$. Let $\varepsilon^{\star}>0$ be the minimal zero of $f(\varepsilon)$. Generically we expect that for a given perturbation size $\varepsilon<\varepsilon^{\star}$, the eigenvalues ${\lambda_{k}(W+\varepsilon E_{\varepsilon,c})>0}$ and $\lambda_{k+1}(W+\varepsilon E_{\varepsilon,c})>0$ are simple. If so, then $f_{c}(\varepsilon):=F_{\varepsilon,c}(E_{\varepsilon,c})$ is a smooth function of $\varepsilon$ and we can exploit its regularity to obtain a fast iterative method to converge to $\varepsilon^{\star}$ from the left. Otherwise we can use a bisection technique to approach $\varepsilon^{\star}$.

Assumption 4.1

We assume that the $k$th and $(k+1)$-st smallest eigenvalues of the Laplacian $L(W+\varepsilon E_{\varepsilon,c})$ are simple. Moreover, $E_{\varepsilon,c}$ is assumed to be a smooth function of $\varepsilon$ in some interval.

We then have the following result.

Lemma 4.4

Under Assumption 4.1, the function $f_{c}(\varepsilon)=F_{\varepsilon,c}(E_{\varepsilon,c})$ is differentiable and its derivative equals (with ${\phantom{a}{}^{\prime}}=d/d\varepsilon$)

[TABLE]

Proof. Differentiating $f_{c}(\varepsilon)=F_{\varepsilon,c}(E_{\varepsilon,c})$ with respect to $\varepsilon$ we obtain as in Theorem 4.1 and noting that $\frac{d}{d\varepsilon}(\varepsilon E_{\varepsilon,c})=E_{\varepsilon,c}+\varepsilon E^{\prime}_{\varepsilon,c}$ with $E^{\prime}_{\varepsilon,c}=\frac{d}{d\varepsilon}E_{\varepsilon,c}$, and setting $G(\varepsilon)=G_{\varepsilon,c}(E_{\varepsilon,c})$, that

[TABLE]

Now we use the property of minimizers, stated by Theorem 4.3, and conclude

[TABLE]

which yields that, with

[TABLE]

since $L(E_{\varepsilon,c})$ is of unit norm for all $\varepsilon$. So we have $f_{c}^{\prime}(\varepsilon)=\langle G(\varepsilon),E_{\varepsilon,c}\rangle$. ∎

Using Theorem 4.3 and Lemmas 4.3 and 4.4, we obtain the following result.

Theorem 4.4

If Assumption 4.1 is satisfied for $\varepsilon\in(0,\varepsilon^{\star})$ and $c$ is sufficiently large so that $W+\varepsilon E_{\varepsilon,c}\geq 0$ for all $\varepsilon$, then $E_{\varepsilon}:=E_{\varepsilon,c}$ is independent of $c$ and $f(\varepsilon):=F_{\varepsilon}(E_{\varepsilon})=f_{c}(\varepsilon)$ has the following properties:

$f(\varepsilon)$* is positive and strictly monotonically decreasing for $\varepsilon\in(0,\varepsilon^{\star})$, and*

[TABLE] 2. 2.

$f(\varepsilon)=0$* for $\varepsilon\geq\varepsilon^{\star}$.*

Proof. The independence of $c$ is obvious if $W+\varepsilon E_{\varepsilon,c}\geq 0$, since then the penalty term $Q_{\varepsilon}(E_{\varepsilon,c})=0$.

(1) We consider first small $\varepsilon$ near 0. For any $E$ with $\|L(E)\|_{F}=1$ we have

[TABLE]

where $F_{0}$ and $G_{0}$ are defined as $F_{\varepsilon}$ and $G_{\varepsilon}$ for $\varepsilon=0$ (they are then independent of $E$). With the choice

[TABLE]

we clearly have $\langle G_{0},E\rangle=-\|G_{0}\|_{F}^{2}/\|L(G_{0})\|_{F}<0$. Since $f(\varepsilon)=F_{\varepsilon}(E_{\varepsilon})\leq F_{\varepsilon}(E_{0})$ and $G(\varepsilon)=G_{\varepsilon}(E_{\varepsilon})=G_{0}+O(\varepsilon)$, we conclude that for sufficiently small $\varepsilon>0$ we have (with $E(\varepsilon)=E_{\varepsilon}$)

[TABLE]

We now show that this inequality actually holds also for all $\varepsilon\in(0,\varepsilon^{\star})$. By Theorem 4.3 we have

[TABLE]

where $\kappa(\varepsilon)\neq 0$ by Lemma 4.3. We note that this implies

[TABLE]

We then have

[TABLE]

Hence $\langle G(\varepsilon),E(\varepsilon)\rangle$ does not change sign, and so we have (11) for all $\varepsilon\in(0,\varepsilon^{\star})$. By Lemma 4.4, this further yields

[TABLE]

This proves statement (1) of the theorem.

(2) Once we have reached the minimum value $\varepsilon^{\star}$ such that $f(\varepsilon^{\star})=0$ the function remains identically zero. In order to see this recall that $L(W+\varepsilon^{\star}E^{\star})$ (with $E^{\star}=E(\varepsilon^{\star})$ the extremizer) has coalescent eigenvalues $\lambda_{k}=\lambda_{k+1}$. Then consider the perturbation, for arbitrary $\delta>-1$,

[TABLE]

where $\varepsilon>0$ is determined by the normalization condition $\|L(E)\|_{F}=1$. Clearly, $L(W+\varepsilon E)=(1+\delta)L(W+\varepsilon^{\star}E^{\star})$ still has coalescent eigenvalues $\lambda_{k}=\lambda_{k+1}$, so that we have $f(\varepsilon)=F_{\varepsilon}(E)=0$. This implies $\varepsilon\geq\varepsilon^{\star}$, since $\varepsilon^{\star}$ is the smallest solution $\varepsilon$ of the problem. Evidently when $\delta\rightarrow\infty$ then $\varepsilon\rightarrow\infty$, so that $f(\varepsilon)$ vanishes on the whole half-line $\varepsilon\geq\varepsilon^{\star}$. ∎

Remark 4.3

Forming the Laplacian and taking the square of the Frobenius norm on both sides of (12), noting the normalizations $\|L(E)\|_{F}=\|L(E^{\star})\|_{F}=1$ and introducing $W^{\star}=W+\varepsilon^{\star}E^{\star}$, we have

[TABLE]

This inequality holds true for all $\delta>-1$ if and only if the orthogonality relation

[TABLE]

is satisfied. This necessary condition for optimality provides a simple check for a numerically obtained extremizer.

4.7 A Newton bisection method

For $\varepsilon<\varepsilon^{\star}$, we make use of the standard Newton iteration

[TABLE]

to get an updated value $\hat{\varepsilon}$. In a practical algorithm it is useful to couple the Newton iteration (14) with a bisection technique, in the same way as for the method presented in [3]. We adopt a tolerance tol which allows us to distinguish whether $\varepsilon<\varepsilon^{\star}$, in which case we may use the derivative formula and perform the Newton step, or $\varepsilon>\varepsilon^{\star}$, so that we have to make use of bisection. The method is formulated in Algorithm 3.

[ht!]

Algorithm 3 Newton-bisection method

1:$W$, $m_{\max}$ (max number of iterations), tol (tolerance),  $\varepsilon_{\rm lb}$ and $\varepsilon_{\rm ub}$ (starting values for the lower and upper bounds for $\varepsilon^{\star}$), $\varepsilon_{0}\in[\varepsilon_{lb},\varepsilon_{ub}]$, $c$ (regularization parameter)

2:$\varepsilon^{\star}$ (upper bound for the distance), $E(\varepsilon^{\star})$

3:Begin

4:     Compute $E(\varepsilon_{0})$ by applying Algorithm 2  

5:     Set $m=0$ 

6:     while $m\leq m_{\max}$ do

7:         if $f_{c}(\varepsilon_{m})<{\rm tol}$ then

8:              Set $\varepsilon_{\rm ub}=\min(\varepsilon_{\rm ub},\varepsilon_{m})$ 

9:              Set $\varepsilon_{m+1}=(\varepsilon_{\rm lb}+\varepsilon_{\rm ub})/2$  (bisection step)

10:         else

11:              Set $\varepsilon_{\rm lb}=\max(\varepsilon_{\rm lb},\varepsilon_{m})$ 

12:              Compute $f_{c}(\varepsilon_{m})$ and $f_{c}^{\prime}({\varepsilon_{m}})$ 

13:              Compute $\varepsilon_{m+1}=\displaystyle{\varepsilon_{m}-\frac{f_{c}(\varepsilon_{m})}{f_{c}^{\prime}(\varepsilon_{m})}}$  (Newton step)

14:         end if

15:         if $\varepsilon_{m+1}\not\in(\varepsilon_{\rm lb},\varepsilon_{\rm ub})$ then

16:              Set $\varepsilon_{m+1}=(\varepsilon_{\rm lb}+\varepsilon_{\rm ub})/2$

17:         end if

18:         if $m=m_{\max}$ or $\varepsilon_{\rm ub}-\varepsilon_{\rm lb}<{\rm tol}$ then

19:              Return $\varepsilon_{m+1}$ and the interval $[\varepsilon_{\rm lb},\varepsilon_{\rm ub}]$ 

20:              Stop

21:         else

22:              Set $m=m+1$

23:         end if

24:         Compute $E(\varepsilon_{m})$ by applying Algorithm 2  

25:     end while

26:Return $\varepsilon^{\star}=\varepsilon_{m}$

The overall method is formulated in Algorithm 4. In our experience, usually it already terminates with $c_{0}=0$.

[h!]

Algorithm 4 The overall algorithm

1:Matrix $W$ and increasing values $c_{0}=0,c_{1},\ldots,c_{m}$ of penalization parameters

2:$\varepsilon^{\star}$ (upper bound for the distance), $E^{\star}=E(\varepsilon^{\star})$

3:Begin

4:     Set $\ell=0$ 

5:     while $\ell\leq m$ do

6:         Find $\varepsilon_{\ell}=\min_{\varepsilon}\{\varepsilon\,|\,F_{\varepsilon,c_{\ell}}(E)=0\}$ by applying Algorithm 3 

7:         Set $E_{\ell}$ the minimizer (normalized by $\|L(E_{\ell})\|_{F}=1$) 

8:         if $W+\varepsilon_{\ell}E_{\ell}\geq 0$ then

9:              Return $\varepsilon^{\star}=\varepsilon_{\ell}$, $E^{\star}=E_{\ell}$ 

10:              Stop

11:         end if

12:         Set $\ell=\ell+1$

13:     end while

4.8 Effective monotonicity with respect to $\varepsilon$

Assume we have integrated equation (8) with $\varepsilon=\varepsilon_{0}$ and we increase to a new value $\varepsilon_{1}>\varepsilon_{0}$. It may not be possible to take the minimizer $E_{\varepsilon_{0}}$ (with $\|L(E_{\varepsilon_{0}})\|_{F}=1$ and $W+\varepsilon_{0}E_{\varepsilon_{0}}\geq 0$) computed in the previous step as the initial value for the constrained gradient flow (8), since we might no longer have $W+\varepsilon_{1}E_{\varepsilon_{0}}\geq 0$ for $\varepsilon_{1}>\varepsilon_{0}$. We therefore first integrate the free gradient system for $F_{\varepsilon_{1},c}$ with initial value proportional to $E_{\varepsilon_{0}}$,

[TABLE]

that is, the differential equation (8) without the Lagrange multiplier $\kappa$ that enforces the norm constraint. Note that $\|L(E(0))\|_{F}=\varepsilon_{0}/\varepsilon_{1}<1$. Along the solution $E(t)$, the functional $F_{\varepsilon_{1},c}(E(t))$ decreases by Theorem 4.1. For the norm of $L(E(t))$ we have

[TABLE]

Since $W+\varepsilon_{0}E_{\varepsilon_{0}}\geq 0$, we have $G_{\varepsilon_{0},c}(E_{\varepsilon_{0}})=G_{\varepsilon_{0}}(E_{\varepsilon_{0}})$, and from Theorem 4.3 and the proof of Theorem 4.4 we know that

[TABLE]

Hence we obtain

[TABLE]

and so the norm of $L(E(t))$ increases with $t$, at least initially. If $\varepsilon_{1}$ is sufficiently close to $\varepsilon_{0}$, we are thus guaranteed to reach an instant $\bar{t}$ where

[TABLE]

We then continue with the norm-preserving differential equation (8),

[TABLE]

with the initial value $E(\bar{t})$, which then has $L(E(t))$ of unit norm for all $t$.

In this way, alternating between the free and the norm-constrained gradient flow, the computed functional $F_{\varepsilon}(E(t))$ is made decreasing with respect to both $t$ and $\varepsilon$.

4.9 Choice of the inital perturbation size $\varepsilon_{0}$ and the initial perturbation matrix

While one might just start with a random perturbation, an educated guess starts from the normalized free gradient $E^{0}=-G_{\varepsilon}(0)/\|L(G_{\varepsilon}(0))\|_{F}$. We propose to choose as a first guess the scaled spectral gap \varepsilon_{0}=\bigl{(}\lambda_{k+1}\left(L(W)\right)-\lambda_{k}\left(L(W)\right)\bigr{)}/\sqrt{2}, which underestimates the optimal $\varepsilon^{\star}=\delta_{k}(W)$ and therefore leads us into the regime where Newton rather than bisection is used for updating $\varepsilon$ in the outer iteration.

5 Numerical experiments

In this section we compare the behavior of the spectral gap and the structured distance to ambiguity (SDA) as stability indicators. First, we determine the optimal numbers of clusters by the criterion of maximal stability with both stability indicators in a family of stochastic block models with varying edge probabilities, alongside a reduced model with similar behavior. Then we examine their performance on a model with randomly generated numbers near given centers, where the optimal number of clusters is given a priori but is here determined by the criterion of maximal stability with both stability indicators.

5.1 Stochastic block model and reduced model

The stochastic block model (SBM) is a model of generating random graphs that tend to have communities. It is an important model in a wide range of fields ranging from sociology to physics (see [14] and, e.g., [1, 2]).

Here we consider a stochastic block model with the following parameters:

• •

the number $n$ of vertices,

• •

a partition of the vertex set $\{1,\ldots,n\}$ into $r$ disjoint subsets $C_{1},\ldots,C_{r}$ called communities,

• •

a symmetric $(r\times r)$-matrix $P=(p_{ij})$ containing edge probabilities.

The graph is then sampled in the way that any two vertices $v_{i}\in C_{i}$ and $v_{j}\in C_{j}$ are connected with probability $p_{ij}$.

For the stability question considered in this paper, this model can be reduced to a model with $2r$ vertices. Assume that the probability matrix $P$ of the SBM is such that $p_{ii}=1$ for all $1\leq i\leq r$. In this case, any two vertices $v$ and $w$ belonging to the same community are connected. If $P=I$, the SBM thus generates a disconnected graph with $r$ communities, which are complete graphs of size $|C_{1}|,\ldots,|C_{r}|$. We represent this graph by a graph with vertex set $\{1,\ldots,2r\}$ such that the vertices $2k-1$ and $2k$ are connected with weight $|C_{k}|$. The edge probabilites $p_{ij}$ in case $P\neq I$ are then represented by inserting matrices $\mu_{ij}I_{2}$ in the respective part of the weight matrix in the reduced model, where $\mu_{ij}$ is an appropriate function of $p_{ij}$ that takes the community sizes into account.

To illustrate the construction of this reduced model, consider an SBM with $n=300$ vertices, three communities of size $100$ and the probability matrix

[TABLE]

Figure 1 shows the adjacency matrix of one sample of the corresponding SBM. This $(300\times 300)$-matrix is represented by the $(6\times 6)$-matrix

[TABLE]

With regards to clustering stability, we observe a similar behavior for the full and the reduced model, as the following example shows.

Example 5.1

(Chain SBM) We measure the clustering stability when applied to an SBM as described above.

• •

We use $r=8$ communities of size $100$ in the reduced model and use $\mu_{1}I_{2},\ldots,\mu_{r-1}I_{2}$ on the off-diagonal, where

[TABLE]

and $\mu_{1}\in\{2,4,\ldots,100\}$ (matrix (16) shows such a model with $r=3$ and $\mu_{1}=20$). For small values of $\mu_{1}$, we expect a clustering into $r$ communities to be most stable, whereas for increasing $\mu_{1}$, the optimal number of clusters should decrease. We compute the optimal number of clusters $k_{\mathrm{opt}}(\delta)$ provided by the SDA and the optimal number of clusters $k_{\mathrm{opt}}(g)$ as provided by the spectral gaps. Figure 2 shows the results. We observe the expected behavior in both robustness measures, but the SDA tends to opt for a lower number of clusters for smaller values of $\mu_{1}$. Figure 3 shows the measures $\delta_{6},\delta_{7},\delta_{8}$ and $g_{6},g_{7},g_{8}$ for different values of $\mu_{1}$. As we expect, $\delta_{8}$ and $g_{8}$ are decreasing, $\delta_{7}$ and $g_{7}$ are increasing up to a certain point and then decreasing again.

• •

To compare the behavior of the above reduced model to the full SBM, we compute the same values for a SBM with $r=6$ communities of size $30$, edge probabilities $p_{k}=\mu_{k}/100$, where $\mu_{1}\in[0,50]$. Figure 4 shows the resulting optimal number of clusters, Figure 5 shows the distances to ambiguity and the spectral gaps for different values of $p_{1}$. We see how the results are affected by randomness in the graph generation but conclude that the behavior of both models is similar.

It is remarkable that the optimal number of communities provided by both stability indicators does not differ by more than one. This suggests that the spectral gap is a reasonable stability indicator in this case, even though it does not take the special structure of weight matrices into account. We observe a similar behavior in different stochastic block models with a similar structure.

5.2 Randomly distributed numbers around centers

In the previous section we have compared the spectral gap to the distance to ambiguity as a stability measure. We have seen a similar behavior yet we cannot conclude which measure works better or worse. We now present an example where we can measure a success rate of the stability indicators.

Consider a given set of centers $m_{1},\ldots,m_{r}$ with the intention to create clusters $C_{1},\ldots,C_{r}$ around them. We generate a set of $n$ random numbers $x_{i}$ such that for each $i\in\{1,\ldots,n\}$, we first choose an index $j\in\{1,\ldots,r\}$ randomly (uniformly distributed) and then the random number $x_{i}$ which is normally distributed with expectation value $m_{j}$ and variance $1$. If the centers $m_{1},\ldots,m_{r}$ are well separated, we expect that the numbers $x_{1},\ldots,x_{n}$ are naturally clustered into $r$ communities.

For a graph with weights $w_{ij}$ related to the distances $|x_{i}-x_{j}|$, we then compute the SDA $\delta_{k}$ and the spectral gap $g_{k}$ for $k\in\{r_{\min},r_{\max}\}$. Since the data set is constructed such that it should naturally have $r$ communities, we expect that $\arg\max_{k}\delta_{k}=r$ in most cases.

Example 5.2

We generate $250$ samples of $n=120$ random numbers and $r=6$ groups. For each $i=1,\ldots,n$, we first choose a group $j\in\{1,\ldots,r\}$ randomly and then generate a random number $x_{i}$, normally distributed around $m_{j}$, where

[TABLE]

Figure 6 shows one sorted random sample $(x_{1},\ldots,x_{n})$.

In order to represent this data set by a graph, we set, following [18]

[TABLE]

and then use

[TABLE]

The tolerance value (here $\mathrm{tol}=10^{-4}$) is used to neglect small weights and avoid a complete graph.We denote by $g_{\mathrm{opt}}$ and $\delta_{\mathrm{opt}}$ the optimal number of clusters obtained by taking the spectral gaps and the SDA, respectively, as stability indicators in the criterion of maximal stability. We used $r_{\min}=4$ and $r_{\max}=8$ as bounds for the number of communities. Table 1 reports the frequency of occurrence of $k$ optimal clusters for $4\leq k\leq 8$ with both stability indicators for $\alpha=\frac{1}{2}$, and in Table 2 for $\alpha=\frac{1}{4}$. We conclude that the success of recognizing the number of communities strongly depends on the graph representation of the data set, but in both cases the SDA outperforms the spectral gap in finding the correct number of clusters.

6 Conclusion

Stability of spectral clustering for undirected weighted graphs is measured by the structured distance to ambiguity (SDA) that is introduced and studied here. We compare the SDA with the spectral gap as a commonly accepted heuristic stability indicator. The spectral gap is characterized as an unstructured distance to ambiguity (Theorem 3.1). In contrast to the spectral gap, the SDA takes into account the underlying graph structure, or more precisely, the $k$th SDA equals the distance of the given graph Laplacian to the set of Laplacian matrices of the same graph with perturbed weights and with coalescing $k$th and $(k+1)$-st eigenvalues. We derive a two-level algorithm to compute the SDA.

Stability indicators such as the structured or unstructured distance to ambiguity can be used for selecting the appropriate number of clusters, by maximizing the stability indicator (choosing the “most stable” clustering). Our numerical experiments show that the spectral gap, although it ignores the underlying graph structure, often yields clustering results similar to the SDA. Tested on the stochastic block model, the optimal number of clusters delivered by both stability indicators usually differs by at most 1. We also present a particular example where the SDA outperforms the spectral gap in finding the optimal number of clusters.

The at most slightly better results of the SDA for choosing the number of clusters come along with a much higher computational cost. The purpose of this paper is thus not to replace the spectral gap by the SDA for general clustering computations. The objective here was to provide a sound theoretical concept and an additional computational tool, against which clustering results obtained with computationally cheaper heuristic criteria can be tested in special situations.

Acknowledgments

We thank two anonymous Referees for their valuable comments and suggestions.

We thank Ulrike von Luxburg and Matthias Hein for a helpful and encouraging discussion.

Part of this work was developed during some visits to Gran Sasso Science Institute in L’Aquila and to the University of Tübingen. The authors thank both institutions for the very kind hospitality.

N. Guglielmi thanks the INdAM GNCS (Gruppo Nazionale di Calcolo Scientifico) and the Italian M.I.U.R. (PRIN 2017) for financial support.