An MBO scheme for clustering and semi-supervised clustering of signed networks

TL;DR

This paper presents a graph-based diffuse interface model using the Ginzburg-Landau functional and an adapted MBO scheme for effective clustering of signed networks, including semi-supervised scenarios, with promising results on synthetic and real data.

Contribution

It introduces a novel MBO-based clustering method for signed networks that scales well and incorporates semi-supervised learning, improving upon existing approaches.

Findings

Effective clustering of signed networks demonstrated on synthetic and real data.

Method scales efficiently to large sparse networks.

Semi-supervised extension improves clustering accuracy.

Abstract

We introduce a principled method for the signed clustering problem, where the goal is to partition a graph whose edge weights take both positive and negative values, such that edges within the same cluster are mostly positive, while edges spanning across clusters are mostly negative. Our method relies on a graph-based diffuse interface model formulation utilizing the Ginzburg-Landau functional, based on an adaptation of the classic numerical Merriman-Bence-Osher (MBO) scheme for minimizing such graph-based functionals. The proposed objective function aims to minimize the total weight of inter-cluster positively-weighted edges, while maximizing the total weight of the inter-cluster negatively-weighted edges. Our method scales to large sparse networks, and can be easily adjusted to incorporate labelled data information, as is often the case in the context of semi-supervised learning. We tested our method on a number of both synthetic stochastic block models and real-world data sets (including financial correlation matrices), and obtained promising results that compare favourably against a number of state-of-the-art approaches from the recent literature.

Tables

Table 1. Table 1 : Summary of frequently used symbols in the paper from Section 3 onwards.

Symbol	Description
$\alpha$	percentage of known information
$A$	affinity matrix of the signed case
$A^{+}$	positive part of $A$
$A^{-}$	negative part of $A$
$B$	time evolution matrix in the MBO scheme with affinity matrix $N$
$\overline{B}$	time evolution matrix in the MBO scheme with affinity matrix $A$
$c_i$	size of clusters of the Signed Stochastic Block Model
$d\tau$	time discretization parameter
$D$	degree matrix of $N$
$\overline{D}$	degree matrix of $A$
$\underset{¯}{e}$	element of the vertices of the probability simplex
${\underset{¯}{e}}^{\pm }$	element of the vertices of ${\Sigma }_K^{\pm }$
$\epsilon$	parameter of the GL functional
$F_{\epsilon }$	GL functional
${\overline{F}}_{\epsilon }$	signed GL functional
$I_V$	Identity matrix with size $V$
$K$	number of clusters
$\Lambda$	eigenvalue matrix
$L$	Laplacian of $N$
$L^{+}$	Laplacian of $A^{+}$
$L_{\text{rw}}$	random walk Laplacian of $N$
$L_{\text{sym}}$	symmetric Laplacian of $N$
$\overline{L}$	signed Laplacian
${\overline{L}}_{\text{rw}}$	signed random walk Laplacian
${\overline{L}}_{\text{sym}}$	signed symmetric Laplacian
$\lambda$	sparsity probability
${\lambda }^{+}$, ${\lambda }^{-}$	trade-off parameters for must- and cannot-links, respectively
$m$	number of eigenvector considered in the expansion of $\overline{B}$
$n$	step number in the MBO algorithm
$N$	affinity matrix of the positive edges case
$N_{\tau }$	number of repetitions of the diffusion step
${\mathbb{N}}_{+}$	natural numbers excluding $0$.
$\eta$	noise probability

Table 2. Table 2 : Parameter values used to test the MBO scheme in the SSBM experiments.

Parameter	Value
Nodes	$V=1200$
Clusters	$K=2$, $3$, $4$, $5$, $8$, $10$
Size	Equal clusters size $V/K$
$d\tau$	0.1
$N_{\tau }$	3
$m$	Equal to no. of clusters $K$
$ϵ$	${10}^{-7}$
No. Experiments	$100$
Sparsity	$\lambda \in [0.02,0.1]$
Noise	$\eta \in [0.05,0.45]$
Known Information	$\alpha \in [0.05,0.9]$

Table 3. Figure 14 : Clusterings computed by Kmeans++, BNC, and MBO for the signed US Migration network, for varying number of clusters k 𝑘 k . Each color, extracted from a continuous color map, corresponds to a different cluster.

	Kmeans++	BNC	MBO
$k=5$
$k=10$
$k=15$
$k=20$
$k=30$
$k=40$
$k=50$

Equations

$$P_{\alpha }(C):=\alpha {\sum }_n+(1-\alpha ){\sum }_p,$$

$$\bar{L}:=\bar{D}-A,$$

$${\bar{L}}_{\text{rw}}:=I-{\bar{D}}^{-1}A,$$

$${\bar{L}}_{\text{sym}}:=I-{\bar{D}}^{-1/2}A{\bar{D}}^{-1/2},$$

$${\min }_{\{x_1,\dots ,x_k\}\in {\mathcal{I}}_k}\left(\sum _{c=1}^k\frac{x_c^T(D^{+}-A)x_c}{x_c^T{x}_c}\right),$$

$${\min }_{\{x_1,\dots ,x_k\}\in {\mathcal{I}}_k}\left(\sum _{c=1}^k\frac{x_c^T(D^{+}-A)x_c}{x_c^T\bar{D}{x}_c}\right).$$

$$x_{t}(i):=\left\{\begin{array}[]{rl}1,&\;\text{if node }i\in C_{t},\\ 0,&\;\text{otherwise },\\ \end{array}\right.$$

$$F_{\epsilon }\left[\rho \right]:=\frac{\epsilon }{2}{\int }_0^1{\left(\frac{d\rho }{ds}\right)}^{2}ds+\frac{1}{4\epsilon }{\int }_0^1{\left(1-{\rho }^2\right)}^{2}ds,$$

$$\frac{d\rho }{dt}\left(t,x\right)=-\frac{\delta F_{\epsilon }}{\delta \rho }\left[\rho \right]\left(t\right),$$

$$\frac{d\rho }{dx}\left(t,1\right)=\frac{d\rho }{dx}\left(t,0\right)=0\text{for}t>0.$$

$$\frac{d\rho }{dt}\left(t,x\right)=\epsilon \frac{d^2\rho }{d{x}^2}\left(t,x\right)-\frac{1}{\epsilon }\rho \left(t,x\right)\left({\rho }^2\left(t,x\right)-1\right),$$

$$\frac{d{F}_{\epsilon }\left[\rho \right]}{dt}\left(t\right)=-{\int }_0^1{\left(\frac{d\rho }{dt}\left(t,s\right)\right)}^{2}ds<0,$$

$${\rho }_{\epsilon }^{\ast }\left(x\right):=\tanh \left(\frac{x}{\epsilon \sqrt{2}}\right).$$

$$D_{ij}\coloneqq\left\{\begin{array}[]{cc}\begin{array}[]{c}\sum_{k=1}^{V}N_{ik}\in\mathbb{R}_{+0}\\ 0\end{array}&\begin{array}[]{c}\text{if }i=j\\ \text{if }i\neq j\end{array}\end{array}\right.,$$

$$\underline{\rho }\left(t\right):=\left\{{\rho }_0(t),{\rho }_1(t),\dots ,{\rho }_V(t)\right\}.$$

$$F_{\epsilon }(\underline{\rho }):=\frac{\epsilon }{2}⟨\underline{\rho },L\underline{\rho }⟩+\frac{1}{4\epsilon }\sum _{i=1}^V({\rho }_i^2-1{)}^2,$$

$$⟨\underline{\rho },L\underline{\rho }⟩=\sum _{i=1}^V{N}_{ij}({\rho }_i-{\rho }_j{)}^2\ge 0.$$

$$\frac{d{\rho }_i}{dt}\left(t\right)=-\epsilon {\left[L\underline{\rho }\right]}_i\left(t\right)-\frac{1}{\epsilon }{\rho }_i\left(t\right)\left({\rho }_i^2\left(t\right)-1\right),$$

$$\frac{{\underline{\rho }}^{n+\frac{1}{2}}-{\underline{\rho }}^n}{d\tau }=-\epsilon L{\underline{\rho }}^{n+\frac{1}{2}},$$

$${\underline{\rho }}^{n+\frac{1}{2}}={\left(I_V+\epsilon d\tau L\right)}^{-1}{\underline{\rho }}^n,$$

$${\underline{\rho }}^{n+\frac{1}{2}}={\left[X\left(I_V+d\tau \epsilon \Lambda \right)X^T\right]}^{-1}{\underline{\rho }}^n,$$

$$\rho_{i}^{n+1}=T\left(\rho_{i}^{n+\frac{1}{2}}\right)\coloneqq\left\{\begin{array}[]{c}\begin{array}[]{cc}+1&\text{if }\rho_{i}^{n+\frac{1}{2}}\geqslant 0\end{array}\\ \begin{array}[]{cc}-1&\text{if }\rho_{i}^{n+\frac{1}{2}}<0\end{array}\end{array}\right..$$

$$\Sigma_{K}\coloneqq\left\{\left(\begin{array}[]{cccc}\sigma_{1},&\sigma_{2},&\dotsc,&\sigma_{K}\end{array}\right)\in\left[0,1\right]^{K}:\sum_{k=1}^{K}\sigma_{k}=1\right\},$$

$$\underline{u}_{i}\coloneqq\left(\begin{array}[]{cccc}u_{i1},&u_{i2},&\dotsc,&u_{iK}\end{array}\right)^{T}\in\Sigma_{K},$$

$$U\coloneqq\left(\begin{array}[]{cccc}\underline{u}_{1},&\underline{u}_{2},&\dotsc,&\underline{u}_{V}\end{array}\right)^{T}\in\left(\Sigma_{K}\right)^{V}.$$

$$F_{\epsilon }(U):=\frac{\epsilon }{2}⟨U,LU⟩+\frac{1}{2\epsilon }\sum _{i=1}^V\left(\prod _{k=1}^K\frac{1}{4}{\Vert {\underline{u}}_i-{\underline{e}}_k\Vert }_1^2\right),$$

$$\underline{e}_{k}\coloneqq\left(\begin{array}[]{cccc}e_{1},&e_{2},&\dotsc,&e_{K}\end{array}\right)^{T}\in\left\{0,1\right\}^{K},$$

$$U^{n+\frac{1}{2}}=B{U}^n,\text{where}B:={\left[X\left(I_V+\frac{d\tau }{N_{\tau }}\Lambda \right)X^T\right]}^{-N_{\tau }}\in {\mathbb{R}}^{V\times V}$$

$${\underline{v}}_i:=\underset{\underline{y}\in {\Sigma }_K}{\mathrm{argmin}}{\Vert {\underline{u}}_i^{n+\frac{1}{2}}-\underline{y}\Vert }_2\in {\left[0,1\right]}^K.$$

$${\underline{u}}_i^{n+1}={\underline{e}}_k\text{with}k=\underset{j=1,\dots ,K}{\mathrm{argmin}}{\Vert {\underline{v}}_i-{\underline{e}}_j\Vert }_2.$$

Full text

1. 1 Introduction

2. 2 Related literature

3. 3 Ginzburg–Landau functional

4. 3.1 GLF as potential energy

5. 3.2 GLF and clustering

6. 3.3 Numerical solution of the ACE

7. 3.4 Generalization to multiple clusters

8. 4 Ginzburg–Landau functional on signed graphs

9. 4.1 Signed Laplacians

10. 4.2 MBO scheme for signed graphs

11. 5 GLF with constraints

12. 5.1 Must-link and cannot-link

13. 5.2 Fidelity and avoidance terms

14. 5.3 Anchors

15. 6 The algorithm

16. 7 Numerical experiments for a signed stochastic block model

17. 7.1 Best practice for MBO

18. 7.2 Comparison to other methods

19. 7.3 Must-link and cannot-link

20. 7.4 Fidelity and avoidance terms

21. 7.5 Anchors

22. 8 Numerical experiments for the Barabási–Albert model

23. 9 Image segmentation

24. 10 Time series correlation clustering

25. 10.1 S&P 1500 Index

26. 10.2 Foreign exchange

27. 11 US migration

28. 12 Summary and future directions

29. Appendix

30. A.1 Projection and Thresholding on $\Sigma_{K}$ and $\Sigma_{K}^{\pm}$

Keywords: MBO, clustering, signed networks, graph Laplacians, spectral methods, time series

1 Introduction

Clustering is one of the most widely used techniques in data analysis, and aims to identify groups of nodes in a network that have similar behaviour or features. It is arguably one of the first tools one uses when facing the task of understanding the structure of a network.

Spectral clustering methods have become a fundamental tool with a broad range of applications in applied mathematics and computer science, and tremendous impact on areas including scientific computing, network science, computer vision, machine learning and data mining, such as ranking and web search algorithms. On the theoretical side, understanding the properties of the spectrum (both eigenvalues and eigenvectors) of the adjacency matrix and its various corresponding Laplacians, is crucial for the development of efficient algorithms with performance guarantees, and leads to a very mathematically rich set of open problems.

Variational methods have a long history in mathematics and physics. In a variational model, the aim is to find minimizers of a given functional (which is often called an ‘energy’ in recognition of the field’s origins, which are rooted in physics, even in situations where this terminology is not physically justified). When formulated in a continuum setting, such models are closely related to partial differential equations (PDEs) which arise as the Euler–Lagrange equations or gradient flow equations in the model. These equations are often nonlinear and can have other interesting properties, such as nonlocality. In recent years it has been recognized that such properties can be useful for clustering problems as well [1, 2]. When attempting a hard clustering problem, i.e. a clustering problem in which each node is to be assigned to exactly one cluster, the approach is usually to relax the exact combinatorial problem. For example, in the case of spectral relaxations, one typically leverages the spectrum of certain graph-based matrix operators, in particular their eigenvalues and especially eigenvectors. This introduces the question of how to translate the real-valued outcomes of the relaxed problem into a discrete cluster label. To this end, non-linearities can prove useful as they can be used to drive the real-valued approximate solution away from the decision boundaries. Moreover, what is non-local in a continuum setting, can be made local in a graph setting, by adding edges between nodes that are located far away in continuum space. Furthermore, since the equations which result from graph-based variational methods tend to resemble discretized PDEs, the large existing literature on numerical methods for solving PDEs often provides fast algorithms for solving such equations. Considerations of similar nature have led to a growing interest in the use of variational methods in a graph setting [3, 4, 5, 6].

In certain applications, negative weights can be used to denote dissimilarity or distance between a pair of nodes in a network. The analysis of such signed networks has become an increasingly important research topic in the past decade, and examples include social networks that contain both friend and foe links, shopping bipartite networks that encode like and dislike relationships between users and products [7], and online news and review websites such as Epinions [8] and Slashdot [9] that allow users to approve or denounce others [10]. Such networks may be regarded as signed graphs, where a positive edge weight captures the degree of trust or similarity, while a negative edge weight captures the degree of distrust or dissimilarity. Further efforts for the analysis of signed graphs have lead to a wide range of applications including edge prediction [11, 10], node classification [12, 13], node embeddings [14, 15, 16, 17], and node ranking [18, 19]. We refer the reader to [20] for a recent survey on the topic.

Finally, another main motivation for our proposed line of work arises in the area of time series analysis, where clustering time series data is a task of particular importance [21]. It arises in a number of applications in a diverse set of fields ranging from biology (analysis of gene expression data [22]), to economics (time series that capture macroeconomic variables [23]). A particularly widely used application originates from finance, where one aims to cluster the large number of time series corresponding to financial instruments in the stock market [24, 25]. In all the above application scenarios, a popular approach in the literature is to consider the similarity measure given by the Pearson correlation coefficient (and related distances), that measures the linear dependence between pairs of variables and takes values in $[-1,1]$. Furthermore, the resulting empirical correlation matrix can be construed as a weighted network, where the signed edge weights capture pairwise correlations, which reduces the task of clustering the multivariate time series to that of clustering the underlying signed network. For scalability and robustness reasons, one can further consider thresholding on the p-value associated to each individual sample correlation entry [26], following tests of statistical significance, which renders the network sparse, and thus amenable to larger scale computations. We refer the reader to the seminal work of Smith et al. [27] for a detailed survey and comparison of various methodologies for casting time series data into networks, in the context of fMRI time series. Their overall conclusion that correlation-based approaches are typically quite successful at estimating the connectivity of brain networks from fMRI time series data, provides sound motivation for the development of algorithms for analysis of signed networks and correlation clustering.

Contribution. Our main contribution is to extend the Ginzburg–Landau functional to the setting of clustering signed networks, and show the applicability of the MBO scheme in this context, for different types of constraints that arise in the semi-supervised setting. As a side contribution, we provide a more computationally efficient way to determine node-cluster membership (as an alternative to the traditional projection/thresholding step). In addition, we show extensive numerical experiments on both synthetic and real data, showcasing the competitivity of our approach when compared to state-of-the-art methods, across various noise and sparsity regimes.

Paper outline. The remainder of this paper is organized as follows. Section 2 is a summary of related work from the signed clustering literature. Section 3 gives a brief overview of the Ginzburg–Landau functional. Section 4 extends the Ginzburg–Landau framework to the setting of signed graphs, by explicitly handling negative weights and modifying the MBO algorithm accordingly. Section 5 pertains to the semi-supervised learning setting, and incorporates various types of a priori information available to the user. Section 6 summarizes the structure of the algorithm we propose. Section 7 contains numerical experiments under a signed stochastic block model, and includes various scenarios in the semi-supervised setting. Section 8 details the outcome of numerical experiments under the Barabási–Albert model. Section 9 considers an application to image segmentation, Section 10 to clustering correlation matrices arising from financial time-series data, and Section 11 to computational social science using a network resulting from migration patterns between US counties. Section 12 summarizes our findings and details future research directions. Finally, we defer to the Appendix the technical details that pertain to the projection/thresholding step.

Notation. We summarize our notation in Table 1.

2 Related literature

The problem of clustering signed graphs traces its roots to the work of Cartwright and Harary from the 1950s from social balance theory [28, 29]. Consider, for example, three mutually connected nodes from a signed graph, and define the relationship between two nodes to be positive if the edge in-between is $+1$ (they are “friends”), and negative if the edge in-between is $-1$ (they are “enemies”). We say that such a triad is balanced if the product of its signs is positive, and that the signed graph is balanced if all of its cycles are positive [29], i.e., each cycle contains an even number of negative edges. Cartwright and Harary proved the so-called structure theorem that if a signed graph $G=(V,E)$ is balanced, then the set of vertices, $V$ can be partitioned into two subsets, called plus-sets, such that all of the positive edges are between vertices within a plus-set and all negative edges are between vertices in different plus-sets.

Davis [30] proposed a weaker notion of balance, weak balance theory, which relaxed the balanced relationship by allowing an enemy of one’s enemy to also be an enemy, and showed that a network with both positive and negative edges is k-weakly balanced if and only if its nodes can be clustered into $k$ groups such that edges within groups are positive and edges between groups are negative. For a $k$-balanced network, the problem of finding the partition of the node set can be viewed as an optimization problem: find the $k$-partition $C$ which minimizes $P_{\alpha}:\Phi\rightarrow\mathbb{R}$ given by

[TABLE]

over $\Phi$, the set of all possible $k$-partitions. Here $\sum_{p}$ is the number of positive arcs between plus-sets, $\sum_{n}$ is the number of negative arcs within plus-sets, and $\alpha\in\left(0,1\right)$ is a weighting constant [31]. If the graph is $k$-balanced the set of minimizers of $P$ will not change when $\alpha$ varies; if the graph is not $k$-balanced, the parameter $\alpha$ allows control over which kind of error is penalized more strongly, positive ties between groups or negative ties within groups.

Motivated by this characterization, the $k$-way clustering problem in signed networks amounts to finding a partition into $k$ clusters such that most edges within clusters are positive and most edges across clusters are negative. As an alternative formulation, one may seek a partition such that the number of violations is minimized, i.e., negative edges within the cluster and positive edges across clusters. In order to avoid partitions where clusters contain only a few nodes, one often prefers to incentivize clusters of similar large size or volume, leading to balancing factors such as the denominators $x_{c}^{T}x_{c}$ in (5) or $x_{c}^{T}\bar{D}x_{c}$ in (6), discussed further below.

Weak balance theory has motivated the development of algorithms for clustering signed networks. Doreian and Mrvar [32] proposed a local search approach in the spirit of the Kernighan–Lin algorithm [33]. Yang et al. [34] introduced an agent-based approach by considering a certain random walk on the graph. Anchuri et al. [35] proposed a spectral approach to optimize modularity and other objective functions in signed networks, and demonstrated its effectiveness in detecting communities from real world networks. In a different line of work, known as correlation clustering, Bansal et al. [36] considered the problem of clustering signed complete graphs, proved that it is NP-complete, and proposed two approximation algorithms with theoretical guarantees on their performance. On a related note, Demaine and Immorlica [37] studied the same problem but for arbitrary weighted graphs, and proposed an O($\log n$) approximation algorithm based on linear programming. For correlation clustering, in contrast to k-way clustering, the number of clusters is not given in advance, and there is no normalization with respect to size or volume.

Kunegis et al. [38] proposed new spectral tools for clustering, link prediction, and visualization of signed graphs, by solving a signed version of the 2-way ratio-cut problem via the signed graph Laplacian [39]

[TABLE]

where $\bar{D}$ is the diagonal matrix with $\bar{D}_{ii}=\sum_{i=1}^{n}|d_{ij}|$. The authors proposed a unified approach to handle graphs with both positive and negative edge weights, motivated by a number of tools and applications, including random walks, graph clustering, graph visualization and electrical networks. The signed extensions of the combinatorial Laplacian also exist for the random-walk normalized Laplacian and the symmetric graph Laplacian

[TABLE]

[TABLE]

the latter of which is particularly suitable for skewed degree distributions. Here $I$ denotes the identity matrix. Recent work, by a subset of the authors in our present paper, provided the first theoretical guarantees, under a signed stochastic block model, for the above Signed Laplacian $\bar{L}$ [40].

In related work, Dhillon et al. [41] proposed a formulation based on the Balanced Ratio Cut objective

[TABLE]

and the closely related Balanced Normalized Cut (BNC) objective for which the balance factors $x_{c}^{T}x_{c}$ in the denominators are replaced by $x_{c}^{T}\bar{D}x_{c}$,

[TABLE]

Here $\mathcal{I}_{k}$ denotes the collection of all sets containing the indicator vectors of a $k$-partition of the node set $V$, i.e. if $\{x_{1},\ldots,x_{k}\}\in\mathcal{I}_{k}$, then each of the $k$ vectors has entries

[TABLE]

for a collection of disjoint clusters $C_{t}$ which partitions the node set $V$. The same authors claimed that the approach of Kunegis et al. [38] via the Signed Laplacian faces a fundamental weakness when directly extending it to $k$-way clustering, and obtained better results on several real data sets via their proposed algorithm which optimizes the BNC objective function. On the other hand, [40] reported good numerical performance of the Signed Laplacian even for higher values of $k$, aside from providing a robustness analysis of it for $k=2$.

A rather different approach for signed graph clustering has been investigated by Hsieh et al. [42], who proposed a low-rank model, based on the observation that the problem of clustering a signed network can be interpreted as recovering the missing structure of the network. The missing relationships are recovered via state-of-the-art matrix completion algorithms (MC-SVP), which first complete the network and then recover the clusterings via the k-means algorithm using the top $k$ eigenvectors of the completed matrix. An alternative method (MC-MF), scalable to very large signed networks, is based on a low-rank matrix factorization approach, which first completes the network using matrix factorization and derives the two low-rank factors $U,H\in\mathbb{R}^{n\times k}$, runs k-means on both $U$ and $V$, and selects the clustering which achieves a lower BNC objective function (6).

Both strong and weak balance are defined in terms of local structure at the level of triangles. As observed early in the social balance theory literature, local patterns between the nodes of a signed graph induce certain global characteristics. Chiang et al. [43] exploited both local and global aspects of social balance theory, and proposed algorithms for sign prediction and clustering of signed networks. They defined more general measures of social imbalance based on cycles in the graph, and showed that the classic Katz measure, ubiquitous for unsigned link prediction, also has a balance theoretic interpretation in the context of signed networks. The global structure of balanced networks leads to a low-rank matrix completion approach via convex relaxation, which scales well to very large problems. This approach is made possible by the fact that the adjacency matrix of a complete k-weakly balanced network has rank 1 if $k\leq 2$, and rank $k$ for all $k>2$ [43].

The recent work of Mercado et al. [44] presented an extended spectral method based on the geometric mean of Laplacians. In another recent work, a subset of the authors showed in [45] that, for $k=2$, signed clustering can be formulated as an instance of the group synchronization problem over $\mathbb{Z}_{2}$, for which spectral, semidefinite programming relaxations, and message passing algorithms have been studied; extensions to the multiplex setting and incorporation of constraints have also been considered. Finally, we refer the reader to [46, 47] for recent surveys on clustering signed and unsigned graphs.

3 Ginzburg–Landau functional

The Ginzburg–Landau functional (GLF) was originally introduced in the materials science literature to study phase separation, and has since also found many uses in image processing and image analysis due to its connections with the total variation functional [48, 49]. In this section, we briefly review the main idea behind the GLF together with a description of the (local) minimizer(s) of the functional in terms of solutions of the Allen–Cahn equation (ACE). Then, we will show how this approach can be adapted to the discrete case of clustering over signed graphs.

3.1 GLF as potential energy

In this introductory presentation, which is based on [50, 51], we restrict to the one dimensional case for simplicity. Consider the function $\rho$ in the Sobolev space $\mathscr{W}^{1,2}\left(\left[0,1\right]\right)$. The GLF $F_{\varepsilon}:\mathscr{W}^{1,2}\left(\left[0,1\right]\right)\mapsto\left[0,\infty\right)$ is defined as

[TABLE]

for some $\varepsilon>0$. The second term has the shape of a double-well and is often labelled $W\left(\rho\right)\coloneqq\left(1-\rho^{2}\right)^{2}\in\mathbb{R}_{+}$. The GLF is interpreted as the potential energy for the physical observable $\rho$. At this stage, $\rho$ is also understood as time dependent and we want to find solutions of

[TABLE]

where $\frac{\delta F}{\delta\rho}$ is the functional derivative of $F$, which lead in their time evolution to the local minima of the GLF. A common choice for the boundary values of the problem (9) is the von Neumann condition

[TABLE]

In this way, calculating the functional derivative of the right-hand side in (9) gives rise to the ACE

[TABLE]

where the last term is often called $W^{\prime}\left(\rho\right)\coloneqq\rho\left(\rho^{2}-1\right)\in\mathbb{R}$. Moreover, it is possible to verify by direct calculation that the GLF is decreasing in time on the solution of the ACE, that is

[TABLE]

as is expected for the evolution of a physical observable along its energy potential. The meaning of the coefficient $\varepsilon$ can be understood as follows. A stationary solution for the ACE, such that $\frac{d\rho}{dt}\left(t,x\right)=0\,\forall t\in\mathbb{R}_{+}$ in (11), is given by

[TABLE]

By taking the first derivative with respect to the spatial coordinate, we can see that the parameter $\varepsilon$ is responsible for the sharpness of the solution (see Figure 1). Thus, given the values of the function $\rho^{*}_{\varepsilon}$ at the boundary of its domain, the ACE can be understood as the evolution of the layer profile connecting the two boundary values of the system.

3.2 GLF and clustering

Using tools from nonlocal calculus, see Section 2.3 in [1], it is possible to write a version of the GLF where the density function is defined on a network. To fix the notation, we denote as $G\left(V,E,N\right)$ a positive-weighted undirected graph with vertex set $V$, edge set $E$ and affinity matrix $N$ which is a symmetric matrix whose entries $N_{ij}$ are zero if $i=j$ or $i\not\sim j$222The symbol $i\sim j$, respectively $i\not\sim j$, denotes that nodes $i$ and $j$ are connected, respectively not connected, by an edge. and nonnegative otherwise. From the affinity matrix, ones defines the diagonal degree matrix as

[TABLE]

where $V$ is the cardinality of the vertex set. When the entries of the affinity matrix are positive, the associated network is often called an unsigned graph. In the same way, the discrete equivalent of the density function is now a column vector on the vertex set, that is $\underline{\rho}:\mathbb{N}\times\left[0,T\right]\mapsto\mathbb{R}^{V}$, written as

[TABLE]

Following [3], it turns out that the network version of (8) for unsigned graphs, also called graph GLF (GGLF), is given by

[TABLE]

where $L\coloneqq D-N\in\mathbb{R}^{V\times V}$ is the (combinatorial) graph Laplacian, and $\langle\cdot,\cdot\rangle$ denotes the usual inner product on ${\mathbb{R}}^{V}$. It is important to notice that the following expression holds true

[TABLE]

The inequality in (17) shows that $L$ is positive semidefinite (i.e., $L$ has non-negative, real-valued eigenvalues, including [math]) and this assures the non-negativity of the GGLF (16). Again, the minimization of the GGLF can shed light on the interpretation of the terms composing the potential. The second term in $F_{\varepsilon}$ will force minimizers $\rho\in{\mathbb{R}}^{V}$ to have entries which are close to $-1$ or $1$ in order to keep $(\rho_{i}^{2}-1)^{2}$ small. The first term in $F_{\varepsilon}$ penalizes the assignment to different clusters of nodes that are connected by highly weighted edges. This means that minimizers of $F_{\varepsilon}$ are approximate indicator vectors of clusters that have few highly weighted edges between them. In the end, we interpret nodes having different sign as belonging to different clusters.

As before, minimizing the GGLF using the method of the gradient descent, we obtain the ACE for networks which we write component-wise as

[TABLE]

where $\left[L\underline{\rho}\right]_{i}$ indicates the $i$-th component of the vector $L\underline{\rho}$. Instead of the combinatorial graph Laplacian $L$, it is possible to use the random walk Laplacian $L_{\text{rw}}:=I-D^{-1}N$ (which is positive semidefinite), or the symmetric graph Laplacian $L_{\text{sym}}:=I-D^{-1/2}ND^{-1/2}$ (symmetric positive semidefinite). In fact, [52] argues for the use of $L_{\text{sym}}$, but no rigorous results are known that confirm this is necessarily the best choice for each application.

3.3 Numerical solution of the ACE

Various methods have been employed to numerically compute (approximate) minimizers of $F_{\varepsilon}$. In [52] a gradient flow approach is used and (18), interpreted as a system of ordinary differential equations (ODEs), is solved. The hope is that the gradient flow dynamics will converge to a local minimizer of $F_{\varepsilon}$. Under some assumptions [53], this hope is theoretically justified [54], but numerically the nonconvexity of $F_{\varepsilon}$ can cause some issues. In [52] such issues were tackled with a convex splitting scheme.

Another strategy which has been developed to approximately solve the ACE for networks consists in using a variation of the Merriman–Bence–Osher (MBO) scheme (also known as threshold dynamics). This was first introduced in [55, 56] in a continuum setting as a proposed algorithm to approximate flow by mean curvature. Following [1], we first neglect the terms $W^{\prime}\left(\rho_{i}\right)$ and discretize (18) in the following way

[TABLE]

where $d\tau\in\mathbb{R}_{+}$ is the length of a small time interval to be tuned and the superscript $n$ indicates the step of the iterative procedure ($n+\frac{1}{2}$ simply specifies an intermediate step between the $n$-th step and the subsequent one which will be labelled as $n+1$). The above leads to the MBO equation (MBOE)

[TABLE]

where $I_{V}\in\left\{0,1\right\}^{V\times V}$ is the identity matrix. After performing the eigendecomposition of the Laplacian, (20) becomes

[TABLE]

where $X$ is a matrix whose columns are the normalized eigenvectors of $L$ and $\Lambda$ is the relevant eigenvalue diagonal matrix. Since the double-well potential pushes the components of $\underline{\rho}$ to be either $+1$ or $-1$, its role can be approximated via a thresholding step, in the following way

[TABLE]

This procedure can be iterated until a convergence criterion is met. In practice, the diffusion step (20) is repeated $N_{\tau}\in\mathbb{N}_{+}$ times before the thresholding step (22). This corresponds to multiplying $\underline{\rho}^{n}$ in (20) by $\left[X\left(I_{V}+\frac{d\tau}{N_{\tau}}\varepsilon\Lambda\right)X^{T}\right]^{-N_{\tau}}$ (note the division of $d\tau$ by $N_{\tau}$).

The graph version of MBO has been used successfully in recent years in various applications, both applied to $F_{\varepsilon}$ [1], its multi-class extension [57], and its signless Laplacian variant [58]. In general, we notice that the choice of $\varepsilon$ also poses a problem. Smaller values of $\varepsilon$ are expected to lead to minimizers $\underline{\rho}$ that are closer to being indicator vectors, however the discrete length scales of the problem impose a lower bound on the $\varepsilon$ values that can be resolved. The MBO scheme avoids this problem since the $W^{\prime}\left(\rho_{i}\right)$ are substituted with a thresholding step. In this way, we can fix $\varepsilon=1$. Using the MBO scheme instead of the gradient flow does not solve all problems: It eliminates the terms in the ODEs that originated with the non-convex term in $F_{\varepsilon}$ and it removes $\varepsilon$ from consideration, but instead a value for $d\tau$ has to be chosen. The discrete length scales in the problem disallow values that are too small, yet $d\tau$ has to be small enough to prevent oversmoothing by the diffusion term $\left[L\underline{\rho}\right]_{i}$.

3.4 Generalization to multiple clusters

The Ginzburg–Landau method has been adapted to allow for clustering or classification into $K>2$ clusters [59, 60, 61, 62]. We briefly review this line of work here, since our proposed approach for clustering signed graphs is for general $K$ as well. This construction relies upon an assignment or characteristic matrix whose entries correspond to the weight for a node $i$ to belong to a cluster $k$. More precisely, consider the simplex

[TABLE]

which corresponds to row vectors whose components sum up to $1$. Then, we can associate to each node a characteristic vector333We generally refer to the vector whose entries are the “intensity” of belongingness to a cluster as characteristic vector also if these coefficients are not only [math] or $1$., written as a column vector

[TABLE]

which specifies the weight for the node $i$ to belong to each of the clusters in the network. So, the characteristic matrix of all the nodes in the system is defined as

[TABLE]

In the end, as explained in [61], the generalization of the GLF reads as

[TABLE]

where $\left\|\cdot\right\|_{1}$ is the taxicab metric, $\langle U,LU\rangle=\text{Tr}\left(U^{T}LU\right)$ and $\underline{e}_{k}$ correspond to the canonical basis vector defined as

[TABLE]

whose components satisfy $e_{i=k}=1$ and $e_{i\neq k}=0$. In this context, the MBO scheme for the characteristic matrix reads as

[TABLE]

where we have fixed $\varepsilon=1$, as explained at the end of Section 3.3, and we have included the contribution from the diffusion iteration directly in $B$ as explained below equation (22). Then, following [61], we project each row of $U^{n+\frac{1}{2}}$ on the simplex $\Sigma_{K}$. We call the corresponding vector

[TABLE]

It is important to notice that in the intermediate step $\underline{u}_{i}^{n+\frac{1}{2}}\in\mathbb{R}^{K}$. The thresholding step consists in determining the closest vertex to $\underline{v}_{i}$ of $\Sigma_{K}$, that is444When $\underline{v}_{i}$ has the same distance to multiple vertices, one of those is picked at random with uniform probability.

[TABLE]

In the end, we have U^{n+1}=\left(\begin{array}[]{cccc}\underline{u}_{1}^{n+1},&\underline{u}_{2}^{n+1},&\dotsc,&\underline{u}_{V}^{n+1}\end{array}\right)^{T}. In Section 4.2, we will use a simplified version of the MBO scheme which reduces the projection (29) and thresholding (30) to a single step (the equivalence of these two versions will be proved in the Appendix A.1).

4 Ginzburg–Landau functional on signed graphs

In this section, we define a new affinity matrix that also incorporates negative weights, and investigate how the degree and Laplacian matrices can be suitably modified. In this way, we also show how negative weights can be interpreted in the context of the GL functional. Moreover, we discuss how the standard MBO algorithm can be simplified.

4.1 Signed Laplacians

In the following, we will consider connected weighted undirected simple graphs (hence with symmetric affinity matrices and no self-loops). The corresponding affinity matrix $A$ is a symmetric matrix whose entries $A_{ij}$ are zero if $i=j$ or $i\not\sim j$ and nonzero real otherwise. This can be decomposed in two matrices $A^{+}\in\mathbb{R}_{+0}^{V\times V}$ and $A^{-}\in\mathbb{R}_{+0}^{V\times V}$ whose components are defined as

[TABLE]

respectively, so that

[TABLE]

In this decomposition, note that $A^{+}$, respectively $A^{-}$, takes into account the affinity, respectively antipathy, between nodes. Along the same lines as [38], we define the signed degree matrix as

[TABLE]

where $\left|A_{ik}\right|$ is the absolute value of $A_{ik}$. As for the affinity matrix, $\bar{D}$ can be decomposed into $D^{+}\in\mathbb{R}_{+0}^{V\times V}$ and $D^{-}\in\mathbb{R}_{+0}^{V\times V}$, with components

[TABLE]

so that

[TABLE]

Following [38], the signed Laplacian is defined as

[TABLE]

together with its variants

[TABLE]

and

[TABLE]

Kunegis et al. [38] showed that the signed Laplacian $\bar{L}$ is positive semidefinite, and furthermore, unlike the ordinary Laplacian $L$, the signed Laplacian $\bar{L}$ is strictly positive-definite for certain graphs, including most real-world networks. It is known that (see for instance [63, 64]) if a matrix $Y_{1}$ is positive semidefinite then $Y_{2}\coloneqq XY_{1}X^{T}$ is also positive semidefinite provided that $X$ is invertible. We can apply this to show from the definition (38) that $\bar{L}_{\text{sym}}$ is also positive semidefinite since $\bar{D}$ is diagonal with rank $V$. Positive semidefiniteness of $\bar{L}_{\text{rw}}$ follows from the second equality in (38) since it shows that $\bar{L}_{\text{rw}}$ and $\bar{L}_{\text{sym}}$ are similar matrices.

In particular, for (36) we can see that a natural decomposition stemming from (32) and (35) occurs,

[TABLE]

where we stress that $L^{+}\in\mathbb{R}^{V\times V}$ and $Q^{-}\in\mathbb{R}_{+0}^{V\times V}$. The latter is known in the literature as signless Laplacian [65] and, to the best of our knowledge, first appeared, although without this nomenclature, in [66].

To understand the meaning of the different terms contributing to the signed Laplacian, first consider the following quadratic form

[TABLE]

for some $\underline{x}\in\mathbb{R}^{V}$. Using the definition of $D^{-}$ together with the fact that $A^{-}$ is symmetric with null diagonal, we have for the first term555Here, we use the following shorthand notation $\sum_{j<i}\equiv\sum_{\begin{subarray}{c}j=1\\ j<i\end{subarray}}^{V}$, together with $\sum_{j>i}\equiv\sum_{\begin{subarray}{c}j=1\\ j>i\end{subarray}}^{V}$ and $\sum_{j\neq i}\equiv\sum_{j<i}+\sum_{j>i}$.

[TABLE]

and for the second

[TABLE]

Putting all the contributions together gives

[TABLE]

The minimum of (43) is achieved when $x_{i}=-x_{j}$. In the simplified case in which the node vector is limited to dichotomic values such as $\underline{x}\in\left\{-1,1\right\}^{V}$, one can envision the system to be made up of $K=2$ clusters and the value of $x_{i}$ (being only $+1$ or $-1$) specifies the node-cluster association. In addition, with this particular choice, (43) counts the number of times, weighted for $A^{-}_{ij}$, two nodes were placed incorrectly in the same cluster (which we want to minimize). We note that choosing as dichotomic values $\left\{0,1\right\}$ would reach a different conclusion since the result of the minimization would be $x_{i}=x_{j}=0$, that is, the nodes would belong to the same cluster. In regard to $L^{+}$ (being the Laplacian of the affinity matrix $A^{+}$), using either $\left\{-1,1\right\}$ or $\left\{0,1\right\}$ leads to the same interpretation as for (17). A more rigorous point of view based on the TV measure can be found in [58] and references therein.

The elements of the canonical basis, which correspond to pure node-cluster associations, can be suitably transformed to take values $\left\{-1,1\right\}$, via the following expression

[TABLE]

where the associated simplex is

[TABLE]

In addition, since each row of $U$ from (25) is a linear combination of $\underline{e}_{k}$ whose coefficients sum up to $1$, the characteristic matrix in the new representation can be obtained as

[TABLE]

In the end, by plugging $\bar{L}$ in place of $L$ in (26), we have that the signed GL functional is given by

[TABLE]

where the different normalization of the terms comes from (44) and (46). In other words, through the signed GL functional we want to minimize the presence of positive links between clusters as well as negative links within each one. We can verify that the following equality holds

[TABLE]

where, on the right-hand side, the Laplacian of the GL functional is $L^{+}$, and we used the fact that $L^{+}1_{V}=0$ together with known properties of the trace of the product of two matrices.

4.2 MBO scheme for signed graphs

The MBO scheme arising from the signed GL has the same structure as in (28), that is

[TABLE]

together with $\bar{X}$ and $\bar{\Lambda}$ which are the eigenvector and eigenvalue matrices of $\bar{L}$, respectively666Here we have used another definition for $\bar{B}$ which is different from the definition of $B$ in (28). In fact, since $\bar{L}$ is real and symmetric its eigenvector matrix is orthogonal, that is $\bar{X}^{T}=\bar{X}^{-1}$. So, we have that $\left[\bar{X}\left(I_{V}+\frac{d\tau}{N_{\tau}}\bar{\Lambda}\right)\bar{X}^{T}\right]^{-N_{\tau}}=\left[\left(I_{V}+\frac{d\tau}{N_{\tau}}\bar{\Lambda}\right)\bar{X}^{T}\right]^{-N_{\tau}}\bar{X}^{-N_{\tau}}=\left(\bar{X}^{T}\right)^{-N_{\tau}}\left(I_{V}+\frac{d\tau}{N_{\tau}}\bar{\Lambda}\right)^{-N_{\tau}}\bar{X}^{-N_{\tau}}=\left[\bar{X}\left(I_{V}+\frac{d\tau}{N_{\tau}}\bar{\Lambda}\right)^{-1}\bar{X}^{T}\right]^{N_{\tau}}$.. As explained in the Appendix A.1, the succession of projection (which in [61] is performed using the algorithm developed in [67]) and thresholding steps can be reduced to determine the entry of the largest component of each row of $\left(U^{\pm}\right)^{n+\frac{1}{2}}$, so that777As before, when $\left(\underline{u}_{i}^{\pm}\right)^{n+1}$ has the same distance to multiple vertices, one of those is picked at random with uniform probability.

[TABLE]

However, we note that the calculation of the argmax in (51) is independent of the representation chosen for the characteristic matrix. In fact, using (46) we can see that

[TABLE]

Since $\bar{B}1_{V}$ has constant rows, the location of the largest entry in a row of $\left(U^{\pm}\right)^{n+\frac{1}{2}}$ is the same as in $\bar{B}U^{n}$, that is from (51)

[TABLE]

since $\sum_{h=1}^{V}\bar{b}_{ih}$ is independent of $j$ by definition. The factor $2$ can also be neglected since it does not change the actual node-cluster association. We indicate this equivalence with the symbol

[TABLE]

In other words, after applying $\bar{B}$, the position of the largest entry of $\left(\underline{u}_{i}^{\pm}\right)^{n+1}$, resulting from the argmax in (51), is the same as the one of $\left(\underline{u}_{i}\right)^{n+1}$. In this way, we can keep using a sparse node-cluster association matrix by picking $U$ instead of $U^{\pm}$.

5 GLF with constraints

Now, we extend the GLF to three possible types of constraints. Here, we distinguish between two different types of constraints: soft, when the final result of the algorithm does not necessarily satisfy the original constraint, and hard, when the output matches the constraint.

5.1 Must-link and cannot-link

In practice, the affinity matrix usually results from measurements of a system’s properties (see Section 9 for example). However, in the same spirit as in [68], additional information in the form of both positive (must-link) and negative (cannot-link) edge weights, can also be known and could have a different relative strength with respect to measurements. In quantitative terms, we will embed must-links in $A^{+}$ and cannot-links in $A^{-}$ and generate the new degree matrices $D^{+}$ and $D^{-}$ accordingly. To this end, we define

[TABLE]

where $M\in\mathbb{R}_{+0}^{V\times V}$ is the must-link matrix, and $\lambda^{+}\in\mathbb{R}_{+}$ is a trade-off parameter. The role of $M$ consists in adding links to nodes which should be belong to the same cluster or to increase the weight of a link to reinforce node-cluster associations. In general, it is a symmetric sparse matrix with positive weights (whose values can also be different among each couple of nodes). We indicate the resulting degree matrix for $\mathcal{A}^{+}$, defined in the same fashion as in (33), with $\mathcal{D}^{+}\in\mathbb{R}_{+0}^{V\times V}$. Similarly

[TABLE]

where $C\in\mathbb{R}_{+0}^{V\times V}$ is the cannot-link matrix, and $\lambda^{-}\in\mathbb{R}_{+}$ the trade-off parameter. In the same fashion as $M$, the matrix $C$ can include links to nodes which should not belong to the same cluster or to decrease even more the weight of a link pushing node couples to belong to different clusters. In this case as well, $C$ is a symmetric sparse matrix with positive weights, where the greater the weight the more the two nodes are incentivised to belong to different clusters. Also here we have $\mathcal{D}^{-}\in\mathbb{R}_{+0}^{V\times V}$ as the degree matrix of $\mathcal{A}^{-}$. In the end, we have that the GLF in the presence of must- and cannot-links is given as

[TABLE]

where $\mathcal{L}^{+}\coloneqq\mathcal{D}^{+}-\mathcal{A}^{+}$, $\mathcal{Q}^{-}\coloneqq\mathcal{D}^{-}+\mathcal{A}^{-}$ so that $\mathcal{L}\coloneqq\mathcal{L}^{+}+\mathcal{Q}^{-}$. As in (50), we now need to determine the matrix $\bar{\mathcal{B}}$ resulting from the eigendecomposition of $\mathcal{L}$.

5.2 Fidelity and avoidance terms

Soft constraints can also be implemented at the level of the characteristic matrix mainly in two ways. First, as a fidelity term (see for instance [61]) which attracts the dynamics of the ACE towards a particular characteristic matrix $\hat{U}$. At the level of the GLF this translates into a term which increases the overall energy when the node-cluster association in $\hat{U}$ is not matched. Second, as an avoidance term which pushes the ACE dynamics to avoid particular node-cluster combinations expressed in a characteristic matrix $\tilde{U}$. Again, for the GLF this means an energy increase when the wrong node-cluster association is met. We notice that while the first constraint points the system to a particular combination and increasing the energy every time this condition is not met, the second does the opposite, that is only one combination increases the energy and leaves all the others as equally possible. The GLF for an affinity matrix with positive edges only reads as

[TABLE]

where $\circ$ is the Hadamard product and $R^{fi}$ (and $R^{av}$) is a diagonal matrix whose element $R_{ii}^{fi}$ ($R_{ii}^{av}$) is positive if and only if there exists a $j\in\{1,\ldots,V\}$ such that $\hat{U}_{ij}\neq 0$ ($\tilde{U}_{ij}\neq 0$) and is zero otherwise. At the level of the MBOE, following [61] the first step in the procedure is given by

[TABLE]

We can obtain an equivalent expression for the signed GLF using (46). The result is

[TABLE]

Now the MBOE (50) becomes

[TABLE]

and after the transformation $\eqref{eq:Utrasf}$ becomes

[TABLE]

The factor $\frac{1}{4}$ can be embedded in the definitions of $R^{fi}$ and $R^{av}$. As explained in Section 4.2, since $1_{V}$ has constant rows, the term $\bar{B}1_{V}$ does not modify the result of the argmax problem. Moreover, since the global factor is positive, the presence of the coefficient $2$ in front of the brackets in (62) does not affect which entry is the largest and it can be neglected. So, we can consider only the following terms

[TABLE]

5.3 Anchors

In general, anchors correspond to a subset of the network nodes whose values are fixed once and for all. For example, in the context of sensor network localization (that aims to find coordinates of sensors in $\mathbb{R}^{d}$ from given pairwise distance measurements; typically $d=\{2,3\}$), anchors are sensors that are aware of their location, and anchor-based algorithms make use of their existence when computing the coordinates of the remaining sensors [69]. In the context of the molecule problem from structural biology applications [70] (motivated by NMR), there are molecular fragments (i.e., embeddings of certain subgraphs) whose local configuration is known in advance up to an element of the special Euclidean group SE(3), rather than the Euclidean group E(3). In other words, in the terminology of signed clustering for $k=2$, there exists a subset of nodes (anchors) whose corresponding ground truth assigned value $\{-1,1\}\in\mathbb{Z}_{2}$ (i.e., cluster membership) is known a priori, and can be leveraged to find the values of the remaining non-anchor nodes.

In the context of multiple clusters, anchors are understood as those nodes whose cluster membership is known. From the point of view of the MBO, we can consider anchors as fixed points in the dynamics, that is from (19)888For this presentation we use the notation for the unsigned case, but the arguments apply to $\bar{L}$ as well.

[TABLE]

Anchors are implemented in the following way. Suppose we know that node $z$ always belongs to cluster $k$, that is

[TABLE]

including $n=0$ and the intermediate steps $n+\frac{1}{2}$. To include this a priori information in the MBO scheme, it is sufficient to use at the level of (28) an evolution matrix (to be used in place of the original $B$) which leaves the anchor node unaffected by the dynamics. The latter can be defined as

[TABLE]

In general, to take into account the presence of a second anchor node it is sufficient to modify all the relevant entries of $B^{a}_{ij}$ in the same fashion as in equation (66). This procedure can be repeated for all the anchor nodes.

6 The algorithm

This section is devoted to describing the general structure of the MBO scheme. First, we will show the basic implementation and then we will explain how to modify the relevant structures to incorporate different types of constraints.

The MBO scheme for signed Laplacians we propose is given in Algorithm 1. In addition to the graph variables (i.e., the number of nodes, clusters and the signed Laplacian), we also need to specify

• •

$d\tau\in\mathbb{R}_{+}$ which is the length of the interval of the time discretization procedure as in (19),

• •

$N_{\tau}\in\mathbb{N}$ which is the number of times we repeat the diffusion step in (20),

• •

$m\in\mathbb{N}$ which is the number of eigenvectors and eigenvalues considered in the eigendecomposition of $\bar{B}$ in (50) (we refer the reader to the next paragraph for a theoretical motivation of this choice).

The corresponding eigenvector and eigenvalue matrices (thus having $m$ eigenvectors and eigenvalues) are indicated as $\bar{X}_{m}$ and $\bar{\Lambda}_{m}$, respectively. The parameter $\epsilon$ is defined later in (67). In regard to $d\tau$, in Section 7.1 we show, empirically on a synthetic data set, that the outcomes of algorithm 1 are independent of the chosen value as long as it remains small. In [61], it is reported that the optimal value for the number of diffusion step iterations is $N_{\tau}=3$. Finally, the number of eigenvectors considered is chosen to be equal to the number of clusters set a priori since, at a heuristic level, the eigenvalues empirically exhibit a big jump after the $K$-th term (although this may not always hold true, see [71]). On a related note, we refer the reader to the celebrated result of Lee, Gharan, and Trevisan on higher-order Cheeger inequalities for undirected graphs with $k\geq 2$ clusters [72], which motivates the use of the top $k$ eigenvectors and generalizes initial results on $k=2$ [73]. Moreover, the version of the algorithm we propose does takes into account the projection/thresholding step in a different fashion, as explained in Section 4.2, which allows us to shorten the running time. To include various constraints, it suffices to modify only a few relevant parts of the code in algorithm 1. In particular, we consider

• •

Must-link and cannot-link soft constraints - here, we just need to use $\mathcal{B}$ in place of $\bar{B}$ in (54),

• •

Fidelity and avoidance terms - the evolution equation given by $\bar{B}$ now includes the additional terms as in (59),

• •

Anchor nodes - the matrix $\bar{B}$ is modified, as explained in Section 5.3, to include the fixed nodes in the MBO dynamics.

In all cases, as a stop criterion we use, as in [61],

[TABLE]

All the constraints can also be combined together since they act at different levels in the scheme. However, we notice that including anchors automatically rules out possible soft constraints acting on the same node.

7 Numerical experiments for a signed stochastic block model

This section is devoted to present the performance of the MBO scheme for signed Laplacians on synthetic data sets. This approach is useful since it allows a direct comparison of the results with the ground truth. Here, we consider a signed stochastic block model (SSBM) followed by the definition of the cost function used to compare the output of the MBO algorithm with the ground truth.

Given a set of numbers $\left\{c_{1},c_{2},\dotsc,c_{K}\right\}\in\mathbb{N}$ representing the size of $K$ clusters of $V$ nodes, the corresponding affinity matrix for the fully connected case is given in block form as

[TABLE]

where $1_{c_{i}}$ is a square matrix of size $c_{i}$ whose entries are all $1$ and $-1_{c_{i}\times c_{j}}\in\left\{-1\right\}^{c_{i}\times c_{j}}$ is a matrix with $c_{i}$ rows and $c_{j}$ columns whose elements are all $-1$ (for consistency, we have that $\sum_{k=1}^{K}c_{k}=V$). Here, we see that all the nodes belonging to the same (respectively different) cluster have a $+1$ (respectively $-1$) link. The corresponding output vector (see Algorithm 1) for the ground truth is given as

[TABLE]

where each entry corresponds to a node and its value to the cluster (indexed from $1$ to $K$) the node belongs to, and its characteristic matrix is

[TABLE]

An instance of the SSBM is obtained from (68) by introducing sparsity (understood as null entries) and noise (which corresponds to a change of the sign of the elements $S_{ij}$) with probabilities $\lambda$ and $\eta$, respectively. More specifically, an edge is present independently with probability $\lambda$, while the sign of an existing edge is flipped with probability $\eta$. In quantitative terms, an SSBM matrix $A$, which we fix to be symmetric and null-diagonal, is generated from $S$ via the following probabilistic mixture model. Consider a row index $i=1,\dotsc,V$ and a column index $j=i,\dotsc,V$ (the $j$ index starts from $i$ so to consider only the upper triangular part in the definition of the elements of $A$), we define the affinity matrix as

[TABLE]

Then, the corresponding Laplacian can be computed as detailed in Section 4.1.

As explained in Section 6, the output of the MBO algorithm is a vector $\underline{v}^{*}$ which holds the information related to the node-cluster association in the same spirit as (69). Therefore, it is sensible to compare $\underline{v}^{*}$ with $\underline{s}$ to measure the ability of MBO to retrieve the ground truth clustering. To quantify success in terms of robustness to noise and sampling sparsity, we use the popular adjusted rand index (ARI) defined as follows [74]. Consider the vector $\underline{v}^{*}$ and its associated nodes partition $\mathscr{P}\left(\underline{v}^{*}\right)=\left\{\mathscr{P}_{1}\left(\underline{v}^{*}\right),\mathscr{P}_{2}\left(\underline{v}^{*}\right),\dotsc,\mathscr{P}_{K}\left(\underline{v}^{*}\right)\right\}$ into $K$ clusters, where $\mathscr{P}_{i}\left(\underline{v}^{*}\right)$ is a set whose elements are those nodes belonging to cluster $i$. Also for the ground truth vector, it is possible to construct a partition $\mathscr{P}\left(\underline{s}\right)$ which is in general different from $\mathscr{P}\left(\underline{v}^{*}\right)$. Then, the ARI is defined as

[TABLE]

where $a_{i}$, $b_{j}$ and $p_{ij}$ are extracted from the contingency table

[TABLE]

given

[TABLE]

the number of nodes in common between $\mathscr{P}_{i}\left(\underline{v}^{*}\right)$ and $\mathscr{P}_{j}\left(\underline{s}\right)$.

For the numerical tests, the different types of constraints described in Section 5 can be included by extracting, at random, information from the ground truth. More specifically, $M$ and $C$ can be obtained from $S$ for the must/cannot-link case, and $\hat{U}$ and $\tilde{U}$ from $U_{s}$ for the fidelity/avoidance term etc. We indicate the percentage of information overlapping with the ground truth using the parameter $\alpha\in\left[0,1\right]$. Since we are incorporating constraints using ground truth information, it is clear that the results will be closer to the ground truth for increasing $\alpha$. The purpose of this set of experiments is simply to provide a proof of concept for our approach.

In summary, the numerical experiments we consider are generated using the following set of parameters summarized in Table 2, some of which are also explained in Section 6.

Unless stated otherwise, the initial condition we use is the clustering induced by a segmentation based on a bottom eigenvector, that is an eigenvector associated with the smallest nonzero positive eigenvalue999Since the (signed) Laplacians we use are all positive semidefinite (and not identically zero) [38], a nonzero positive eigenvalue can always be found. When the geometric multiplicity of this eigenvalue is greater than one, we pick at random one of the elements of its eigenspace., of the Laplacian matrix (each experiment is specified by a different Laplacian matrix which will be used for its initial condition as well as the MBO evolution). In particular, we consider the interval given by the smallest and largest entries in the aforementioned eigenvector, and create $K$ equally-sized bins spanning this interval. Finally, each node will be allocated in its corresponding bin according to its value in the eigenvector. The resulting characteristic matrix is used as the initial condition.

7.1 Best practice for MBO

As explained in Section 4.1, one may use different versions of the Laplacians given by (36), (37) and (38). In particular, in the context of spectral clustering of (unsigned) graphs, $L_{\text{rw}}$ and $L_{\text{sym}}$ typically provide better results compared to $L$, especially in the context of skewed degree distributions [71, 75]. We observe a similar behaviour for their respective signed counterparts $\bar{L}_{\text{rw}}$, $\bar{L}_{\text{sym}}$ and $L$, as can be seen in Figure 2. We remark that in subsequent experiments, we will only use $\bar{L}_{\text{sym}}$ as it achieves the highest accuracy. Moreover, we have tested the additional degree of freedom resulting from the time discretization procedure by running MBO with different orders of magnitude for the parameter $d\tau$. Results are detailed in Figure 3, and the overlapping curves illustrate that the results are not sensitive to the choice of the parameter $d\tau$.

7.2 Comparison to other methods

The performance of our MBO algorithm is compared to that of the Balanced Normalized Cut (BNC) [41], and with the results achieved by standard Kmeans$++$ clustering [76, 77]. In this respect, as can be seen in Figure 4, MBO always provides better performance than Kmeans$++$, and improves upon BNC when sparsity and noise are high.

7.3 Must-link and cannot-link

In this section, we detail how MBO can be applied in the context of constrained clustering. We have run the MBO algorithm for $\mathcal{L}=\mathcal{L}^{+}$, corresponding to the must-link scenario, and $\mathcal{L}=\bar{L}+\mathcal{Q}^{-}$, corresponding to the case in which only cannot-links are included as additional information in the definition of the Laplacian. In practice, we have considered the symmetric version of $\mathcal{L}$ in the same fashion as in (38). In the plots shown in Figure 5, we see that, for high sparsity and noise, the constrained clustering improves upon the standard signed Laplacian in both cases. In addition, the two approaches can be combined to obtain better results. In this case, we have included the known information up to $\alpha=0.1$ to not modify considerably the sparsity of the affinity matrix.

7.4 Fidelity and avoidance terms

The second soft constraint we include in the MBO scheme is at the level of the characteristic matrix. The fidelity term pushes the characteristic matrix to be similar to $\hat{U}$ and increases the energy for every row $U$ which is different from $\hat{U}$. Instead, for the avoidance term the energy of the GL functional is increased only when there is an exact match between two rows of $U$ and $\tilde{U}$. This means that, in general, the effect of the avoidance term will be less influential since it provides an energy increase fewer times. This can be seen in Figure 6, where the black curve is always below the red one. We add that, in order to make the effect of the fidelity and avoidance terms more evident, we fix the nonzero diagonal elements of the matrices $R$ and $Q$ in (59) equal to $30$.

7.5 Anchors

In this setting, we are introducing in the MBO scheme a hard constraint, as we need to a priori modify the matrix responsible for the evolution of $U$, as detailed in Section 5.3. As shown in Figure 7, anchors always improve upon the unconstrained signed Laplacian, as expected.

8 Numerical experiments for the Barabási–Albert model

In this section, we apply the MBO scheme to signed networks for which the underlying measurement graph (indicating the presence or absence of an edge) is generated from the Barabási–Albert (BA) model [78]. In a nutshell, the algorithm to generate the BA network begins with a connected network with $V_{0}$ nodes. At each time step, a new node is added to the graph and is connected to $\nu\leqslant V_{0}$ pre-existing nodes. The probability that a link from the new node connects to a generic node $i$ is given by

[TABLE]

where the sum indexed by $j$ runs over all pre-existing nodes. We call $H\in\left\{0,1\right\}^{V\times V}$ the resulting BA network which can be used to build the relevant signed affinity matrix in the following way. As before, consider a row index $i=1,\dotsc,V$ and a column index $j=i,\dotsc,V$ (the $j$ index starts from $i$ so to consider only the upper triangular part in the definition of the elements of $A$), we define the signed BA network as

[TABLE]

where the sparsity of the BA network is included via $H_{ij}$ (the noise probability $\eta$ is understood as in (71)) and $S_{ij}$ are the elements of $S$ in equation (68). We compared the numerical performance of MBO with that of BNC and Kmeans$++$, in the same spirit as in Section 7.1. The results are given in Figure 8, and are consistent with those obtained for the SSBM model.

9 Image segmentation

This section considers a proof-of-concept application of the Ginzburg–Landau functional under different types of constraints, as introduced in Section 5, in the context of image segmentation. As an example, we consider a picture of Cape Circeo in Sabaudia (LT, Italy) [79], Figure 9 (a), in which it is possible to distinguish four main clusters which are the sand, sea, mountain, and sky.

The network associated to the picture has been obtained in the following way. First, each pixel has a one-to-one correspondence with a node in the resulting graph. Second, each pixel is linked in the network only to a subset of its nearest-neighbouring pixels (in particular, we have considered each pixel to be connected to all the others which are present within a radius of $5$, as specified in the MATLAB Package The graph analysis toolbox: Image processing on arbitrary graphs, see [80, 81]). Third, the weights of the affinity matrix corresponding to the picture can be obtained by considering the weight of the link between two pixels as a function of their color difference. In fact, considering the colors red, green, and blue (RGB) as the basis of a 3-dimensional vector space, each color can be described in terms of its RGB components in the same way as a vector can be characterized as the vector sum of the canonical basis elements weighted by some coefficients. In this way, we can naturally define the distance of two colors as a function of the Euclidean distance of their vectors in the RGB space. We used the distance function

[TABLE]

where $b$ is an external parameter (which we fix to $b=14$) and $col\left(i\right)$ is the corresponding color vector of pixel $i$ (note that the norm in (77) is not squared). The set of all the differences $\left\|col\left(i\right)-col\left(j\right)\right\|$ in a picture is normalized according to their z-score using the MATLAB command normalize [82]. The MBO schemes used for this image segmentation problem run using $K=4$ clusters as an input parameter. For constrained clustering, we have considered both must-links and cannot-links, as highlighted in Figure 9. One iteration (multiple iterations do not change the result) of the MBO scheme yields the clustering highlighted in Figure 10 (a). Afterwards, the resulting $U$ has been used as the initial condition for two different constrained clustering optimizations, whose results are shown in Figure 10 (b) and (c) for the must-links and cannot-links, respectively. It can be inferred that the cluster recovery works well, except for a small area in the middle of the picture; however, we remark that even in the original picture (that is Figure 9 (a) without blue and red square) it is difficult to the naked eye to understand whether that patch (that is the bit of sea to the bottom right of the mountain) corresponds to sea or mountain, just by looking at its color.

In Figure 11, we have taken into account fidelity and avoidance terms as highlighted in plot (a). The same information has been used for both the fidelity and the avoidance parts. For the former, the rows of the matrix $\hat{U}$ associated to the nodes in the colored patches have entries equal to $1$ at the corresponding cluster column (we have picked column $1$ for red, $2$ for blue, $3$ for cyan, and $4$ for green). All the remaining entries of $\hat{U}$ are equal to [math]. For the latter, we have considered an energy increase every time nodes are attributed to clusters to which they should not belong. This results in a constraint of the form $\tilde{U}=\hat{U}\left(1_{K}-I_{K}\right)$. Starting from a random configuration of $U$, a single application of the MBO scheme using $\hat{U}$ yields the segmentation shown in Figure 11 (b) and a separate run of the algorithm using $\tilde{U}$ yields to (c). We note that the node-cluster association recovery results agree with the expected boundaries present within the image, with a minor uncertainty in the central area (which is a feature present also in the constrained clustering case).

We have also considered the clusters highlighted in Figure 11 (a) as anchors, and the results are identical to those shown in Figure 11 (b).

10 Time series correlation clustering

This section details the outcomes of our numerical experiments on two time series data sets. Since ground truth data is not available, we compare the output of the algorithms by plotting the network adjacency matrix sorted by cluster membership.

10.1 S&P 1500 Index

We use time series daily prices for $n=1500$ financial instruments, constituents of the S&P 1500 Composite Index, during 2003-2015, containing approximately $n_{d}=3000$ trading days [83]. We consider daily log returns of the prices,

[TABLE]

where $P_{i,t}$ denotes the market close price of instrument $i$ on day $t$. In order to promote both positive and negative weights in the final correlation matrix we aim to cluster, we also compute the daily market excess return for each instrument,

[TABLE]

where $R_{\text{SPY},t}$ denotes the daily log return of SPY, the S&P 500 index ETF which is used as a proxy (weighted average) for the market. Next, we calculate the Pearson correlation between historical returns for each pair of instruments $(i,j)$, which ultimately provides the edge weight $A_{ij}$ in our signed graph. Figure 12 shows a segmentation of the US market into $K\in\{10,30\}$ clusters. The cluster structure is more obvious to the naked eye for the BNC and MBO methods, as opposed to the plain Kmeans$++$ clustering. The large cluster recovered by MBO could potentially correspond to the 1000 instruments which are within the S&P 1500 Index and outside of the S&P 500, but we defer this investigation to future work.

10.2 Foreign exchange

This experiment pertains to clustering a foreign exchange FX currency correlation matrix, based on daily special drawing rights (SDR) exchange value rates [84]. As shown in Figure 13, only BNC and MBO are able to recover four clusters, associated respectively to the EURO €, US Dollar $, UK Pound Sterling £, and Japanese Yen Y

. These are precisely the four currencies that, in certain percentage weights, define the value of the SDR reserve. To the naked eye, the adjacency matrix (sorted by cluster membership) recovered by MBO exhibits sharper transitions between clusters, when compared to BNC and especially Kmeans$++$.

11 US migration

The final real world data set we consider is a migration network based on the 2000 US Census report, which captures migration patterns between all pairs of counties in the US during the 1995-2000 time frame [85, 86]. The network has $V=3107$ nodes, where each node corresponds to a county in mainland US. The adjacency matrix of the given network is asymmetric, with $M_{ij}\neq M_{ji},\;1\leq i,j\leq V$. Similar to the approach in [87], we first symmetrize the input matrix, and let $\tilde{M}_{ij}=M_{ij}+M_{ji}$ denote the total number of people that migrated between county $i$ and county $j$ (thus $\tilde{M}_{ij}=\tilde{M}_{ji}$) during the five-year period. Finally, we turn the network into a signed one by approximately centering the distribution of edge weights at 0. More specifically, from each nonzero entry we remove the median value $\tilde{m}$ of the distribution of nonzero edge weights in $\tilde{M}$, and denote the resulting symmetric matrix by $A$, with $A_{ij}=\tilde{M}_{ij}-\tilde{m}$. In Figure 14, we compare the outputs of plain Kmeans++, BNC, and our MBO approach, for a varying number of clusters $K\in\{5,10,15,20,30,40,50\}$. It can be seen that the plain Kmeans++ approach is not able to return any meaningful results across all values of $k$, and a similar statement holds true for BNC for low values of $K\in\{5,10,15,20\}$: in both cases, a single cluster captures almost all the nodes in the network. BNC returns meaningful results —in the sense that there are no clusters which contain only a few nodes— for larger numbers of clusters $K\in\{30,40,50\}$, aligning well with state boundaries. On the other hand, our proposed MBO algorithm returns meaningful results across all values of $K$, and the clusters align well with the true boundaries between the states, as previously observed in [87].

12 Summary and future directions

We introduced graph-based diffuse interface models utilizing the Ginzburg–Landau functionals, using an adaptation of the classic numerical Merriman–Bence–Osher (MBO) scheme for minimizing such graph-based functionals. Our approach has the advantage that it can easily incorporate various types of labelled data available to the user, in the context of semi-supervised clustering.

We have provided extensive numerical experiments that demonstrate the robustness to noise and sampling sparsity of the MBO scheme, under a signed stochastic block model. In particular, for sparse graphs and larger number of clusters, we often obtain more accurate results compared to state-of-the-art methods such as the BNC and the signed Laplacian.

Our work opens up a number of interesting directions for future investigation. It would be insightful to investigate the performance of the MBO scheme in the setting of very sparse graphs in the stochastic block model with $p=\Theta(\frac{1}{n})$, thus below the connectivity threshold, and potentially incorporate regularization terms considered in the recent literature [88, 89].

As a by-product of our research, we expect that the low-dimensional embeddings of signed networks (e.g., correlation matrices or related similarity/dissimilarity measures) arising from high-dimensional time series data obtained via the MBO scheme, could be of independent interest in the context of robust dimensionality reduction in multivariate time series analysis. Such embeddings could be used to learn nonlinear mappings from multivariate time series data, which has broad applications ranging from financial data analysis to biomedical research. Such mappings could prove useful for exploiting weak temporal correlations inherent in sequential data, and can be subsequently augmented with a regression or classification step, with the end goal of improving out-of-sample prediction. In a financial setting, this corresponds to building nonlinear risk/factor models for financial instruments, a task of paramount importance when it comes to understanding the risks that govern a financial system [90, 91].

An interesting potential future direction is the development of an MBO scheme for detecting structures in directed graphs. The vast majority of methods for structural analysis of graphs, for tasks such as clustering and community detection, have been developed for undirected graphs, and a number of challenges arise when considering directed graphs, complicating the immediate extension of existing algorithms for the undirected case. Various directed extensions typically consider a symmetrized version of the adjacency matrix [92, 93, 94, 95]. In ongoing work by a subset of the authors [96], we consider the problem of clustering directed graphs, where the clustering structure arises due to the existence of certain pairs of imbalanced cuts in the graph. More specifically, an imbalanced cut between two sets of vertices X and Y is one where the majority of the edges are directed from X to Y, with few edges being directed in the opposite direction from Y to X. The goal is to find $k$ clusters such that the sum of all (or many of) the ${k\choose 2}$ pairwise imbalanced cuts is maximized, while also penalizing the existence of small clusters in order to obtain a more balanced partition. In [96], we consider a Hermitian matrix derived directly from the skew-symmetric adjacency matrix, and recover an embedding from the spectrum of its corresponding (Hermitian) graph connection Laplacian, potentially after normalization [69, 97]. Extending the Ginzburg–Landau functional and the MBO scheme to the setting of directed graphs and Hermitian Laplacians is an appealing future research direction.

Acknowledgements

MC and AP acknowledge support from The Alan Turing Institute EPSRC grant EP/N510129/1 and seed funding project SF029 “Predictive graph analytics and propagation of information in networks". AP also acknowledges support from the National Group of Mathematical Physics (GNFM-INdAM), by Imperial College together with the Data Science Institute and Thomson-Reuters Grant No. 4500902397-3408 and EPSRC grant EP/P002625/1. YvG did most of the work which has contributed to this paper when he was a lecturer at the University of Nottingham. YvG acknowledges that this project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 777826. Francesco Cosentino is acknowledged for his contribution in the derivation of the proof of Lemma 2.

Appendix

A.1 Projection and Thresholding on $\Sigma_{K}$ and $\Sigma_{K}^{\pm}$

Let $K\in\mathbb{N}_{+}$. In this section, we show that computation of the projection on $\Sigma_{K}$ (or $\Sigma_{K}^{\pm}$) followed by the thresholding step in the MBO scheme in (29)–(30) is equivalent to determining a maximal component of every vector $\underline{u}_{i}^{n+\frac{1}{2}}$. To fix notation, the vertices of $\Sigma_{K}$ are the elements of the canonical basis in $\mathbb{R^{K}}$ (for instance, for $K=3$ we have $\underline{e}_{1}=\left(1,0,0\right)^{T}$, $\underline{e}_{2}=\left(0,1,0\right)^{T}$ and $\underline{e}_{3}=\left(0,0,1\right)^{T}$), while those of $\Sigma_{K}^{\pm}$ are defined as $\underline{e}_{k}^{\pm}\coloneqq 2\underline{e}_{k}-\underline{1}_{K}\in\left\{-1;1\right\}^{V}$ (that is, for the $K=3$ case, we have $\underline{e}^{\pm}_{1}=\left(1,-1,-1\right)^{T}$, $\underline{e}^{\pm}_{2}=\left(-1,1,-1\right)^{T}$ and $\underline{e}^{\pm}_{3}=\left(-1,-1,1\right)^{T}$).

We first remind the reader of a useful equality for vectors $\underline{a}$, $\underline{b}\in\mathbb{R}^{K}$:

[TABLE]

where we have used the standard inner product in ${\mathbb{R}}^{K}$: $\langle\underline{a},\underline{b}\rangle_{2}:=\sum_{k=1}^{K}a_{k}b_{k}$. The following lemma lists three direct consequences of (80).

Lemma 1**.**

Given $\underline{a}$, $\underline{b}\in\mathbb{R}^{K}$ and $i,j\in\{1,\ldots,K\}$, define $\underline{b}^{*}\in\mathbb{R}^{K}$ by

[TABLE]

Then

[TABLE] 2. 2.

Let $\underline{x}\in\mathbb{R}^{K}$ and let $\underline{e}_{i},\underline{e}_{j}\in\mathbb{R}^{K}$ be two vertices of the simplex $\Sigma_{K}$. Then $\left\|\underline{x}-\underline{e}_{i}\right\|_{2}\geqslant\left\|\underline{x}-\underline{e}_{j}\right\|_{2}$ iff $x_{j}\geqslant x_{i}$, with equality between the norms iff $x_{i}=x_{j}$. 3. 3.

Let $\underline{x}\in\mathbb{R}^{K}$ and let $\underline{e}_{i}^{\pm},\underline{e}_{j}^{\pm}\in\mathbb{R}^{K}$ be two vertices of the simplex $\Sigma_{K}^{\pm}$. Then $\left\|\underline{x}-\underline{e}_{i}^{\pm}\right\|_{2}\geqslant\left\|\underline{x}-\underline{e}_{j}^{\pm}\right\|_{2}$ iff $x_{j}\geqslant x_{i}$, with equality between the norms iff $x_{i}=x_{j}$.

Proof.

We apply (80) with $\underline{c}=\underline{b}^{*}$. Since $\|\underline{b}\|_{2}=\|\underline{b}^{*}\|_{2}$ and

[TABLE]

the result follows. 2. 2.

We apply (80) with $\underline{a}=\underline{x}$, $\underline{b}=\underline{e}_{i}$, and $\underline{c}=\underline{e}_{j}$. Since $\|\underline{e}_{i}\|_{2}=1=\|\underline{e}_{j}\|_{2}$ and $\langle\underline{x},\underline{e}_{j}-\underline{e}_{i}\rangle_{2}=x_{j}-x_{i}$, it follows that

[TABLE]

and thus the result follows. 3. 3.

This follows from a proof completely analogous to the previous one, once we note that $\|\underline{e}^{\pm}_{i}\|_{2}=\sqrt{K}=\|\underline{e}^{\pm}_{j}\|_{2}$ and $\langle\underline{x},\underline{e}^{\pm}_{j}-\underline{e}^{\pm}_{i}\rangle_{2}=2(x_{j}-x_{i})$ where $\underline{x}\in\mathbb{R}^{K}$.

∎

Lemma 2**.**

Let $\underline{x}\in\mathbb{R}^{K}$ and

[TABLE]

be its projection on $\Sigma_{K}$101010Existence and uniqueness of $\underline{x}^{\perp}$ is a well-known result, see for instance Lemma 2.2 in [98]. Given two components $x_{j}$ and $x_{i}$ of $\underline{x}$, if $x_{j}\geqslant x_{i}$ then $x_{j}^{\perp}\geqslant x_{i}^{\perp}$. The analogous result holds for the projection of $\underline{x}$ on $\Sigma_{K}^{\pm}$.

Proof.

We will prove the result for the projection on $\Sigma_{K}$. A completely analogous proof gives the result for $\Sigma_{K}^{\pm}$.

First note that the result is trivially true if $i=j$, hence we assume that $i\neq j$. For a proof by contradiction, assume that $x_{j}\geq x_{i}$ and $x_{j}^{\perp}<x_{i}^{\perp}$. Define $\underline{x}^{\perp*}$ by

[TABLE]

Then, by the first result of Lemma 1 (with $\underline{a}=\underline{x}$ and $\underline{b}=\underline{x}^{\perp}$), we have that

[TABLE]

The inequality follows from our assumptions on $x_{i}$, $x_{j}$, $x_{i}^{\perp}$, and $x_{j}^{\perp}$. Since $\underline{x}^{\perp}\in\Sigma_{K}$, we have $\sum_{k=1}^{K}x^{\perp*}_{k}=\sum_{k=1}^{K}x^{\perp}_{k}=1$, hence $\underline{x}^{\perp*}\in\Sigma_{K}$. Thus, per definition (85), we have $\|\underline{x}-\underline{x}^{\perp*}\|_{2}^{2}-\|\underline{x}-\underline{x}^{\perp}\|_{2}^{2}\geq 0$. It follows that $\|\underline{x}-\underline{x}^{\perp*}\|_{2}=\|\underline{x}-\underline{x}^{\perp}\|_{2}$. By uniqueness of the minimizer in (85) we deduce that $\underline{x}^{\perp*}=\underline{x}^{\perp}$. In particular that means that $x^{\perp}_{j}=x^{\perp*}_{j}=x^{\perp}_{i}$, which is a contradiction. ∎

Applying Lemma 1 (result 2) to the thresholding step in (30), we find that $k$ in (30) is such that the $k^{\text{th}}$ entry of $\underline{v}_{i}$ is not smaller than all its other entries. Furthermore, Lemma 2 applied to (29) tells us that $\underline{v}_{i}$ and $\left(\underline{u}_{i}^{n+\frac{1}{2}}\right)$ have their largest entries at the same positions. Thus (29) and (30) are equivalent to $\underline{u}_{i}^{n+1}=\underline{e}_{k}$, where $k$ is such that the $k^{\text{th}}$ entry of $\left(\underline{u}_{i}^{n+\frac{1}{2}}\right)$ is not smaller than all its other entries.