A Unified Framework for Fair Spectral Clustering With Effective Graph Learning

Graph learning (GL) aims to infer the graph topology behind observed data, a prerequisite step for (fair) SC when similarity graphs are unavailable. Traditionally, graphs are constructed via some direct rules, such as $k$-nearest-neighborhood ($k$-NN), $\varepsilon$-nearest-neighborhood ($\varepsilon$-NN) [22], and sample correlation methods like Pearson correlation (PC). These methods may be limited in capturing similarity relationships between data pairs [23]. Thus, many works attempt to learn graphs from data adaptively, including the sparse representation (SR) method [24] and the low-rank representation method [25]. The emergence of adaptive neighbourhood graph learning (ANGL) [26] provides a new way that uses the probability of two samples being adjacent to measure the similarity between them. In [27], a possibilistic neighbourhood graph is proposed, an improved version of [26]. Recently, with the rise of graph signal processing (GSP) [28], many works attempt to learn graphs from the perspective of signal processing. One of the widely-used GSP-based GL methods postulates that signals are smooth over the corresponding graphs [29]. Intuitively, a smooth graph signal means the signal values of two connected nodes are similar[30], which is also a fundamental principle of SC [5]. Many methods are dedicated to learning graphs from smooth signals [31]. However, limited to our understanding, applying smoothness-based GL to SC has yet to be thoroughly explored, let alone FSC.

Many works focus on establishing a unified model for SC, which can be roughly divided into three categories. The first one integrates graph construction and spectral embedding [21, 32, 26]. They use an independent discretization step as post-processing. The second one is based on a given similarity graph. They integrate spectral embedding and discretization [33, 34, 35]. The third category integrates all three stages into a single objective function [20, 23, 36, 37, 38]. Our model differs from these models in two main ways. (i) Our framework utilizes a new graph construction method. (ii) We further consider fairness issues in clustering tasks.

Given an undirected graph $\mathcal{G}=\{\mathcal{V},\mathcal{E}\}$ of $D$ vertices, where $\mathcal{V}$ and $\mathcal{E}$ are the sets of vertices and edges of $\mathcal{G}$, respectively, its adjacency matrix $\mathbf{W}\in\mathbb{S}^{D\times D}$ is a symmetric matrix with zero diagonal entries and non-negative off-diagonal elements if the graph has non-negative edge weights and no self-loops. The Laplacian matrix of $\mathcal{G}$ is defined as $\mathbf{L}=\mathbf{D}-\mathbf{W}$, where $\mathbf{D}\in\mathbb{S}^{D\times D}$ is a diagonal matrix satisfying $\mathbf{D}_{[ii]}=\sum_{j=1}^{D}\mathbf{W}_{[ij]}$. Unnormalized SC aims to partition $D$ nodes into $K$ disjoint clusters $\mathcal{C}_{1},...,\mathcal{C}_{K}$, where $\mathcal{V}=\mathcal{C}_{1}\cup...\cup\mathcal{C}_{K}$, and $\mathcal{C}_{k}$ is the set containing nodes in the $k$–th cluster. The problem of unnormalized SC is equivalent to minimizing the $\mathrm{RatioCut}$ objective function [5], i.e.,

where $\mathcal{V}\setminus\mathcal{C}_{k}$ contains all nodes in $\mathcal{V}$ except those in $\mathcal{C}_{k}$, and

Let $\widetilde{\mathbf{U}}\in\mathbb{R}^{D\times K}$ be

Then, minimizing the $\mathrm{RatioCut}$ objective function (1) is equivalent to solving the following problem [5]

Due to the discrete constraint of (3), problem (4) is NP-hard. In practice, problem (4) is usually relaxed to

where $\mathbf{U}\in\mathbb{R}^{D\times K}$ is a relaxed continuous clustering label matrix, and $\mathbf{U}^{\top}\mathbf{U}=\mathbf{I}$ is adopted to avoid trivial solutions. The process of solving (5) is called spectral embedding. After obtaining $\mathbf{U}^{*}$, a common practice is to apply $k$means to the rows of $\mathbf{U}^{*}$ to yield final discrete clustering labels $\mathbf{Q}$, where $\mathbf{Q}\in\{0,1\}^{D\times K}$ is a binary cluster indicator matrix. The only non-zero element of the $i$-th row of $\mathbf{Q}$ indicates the cluster membership of the $i$-th node of $\mathcal{G}$.

Fair spectral clustering groups the vertices of $\mathcal{G}$ by considering fairness. If the nodes of $\mathcal{G}$ belong to $S$ sensitive groups $\mathcal{D}_{1},...,\mathcal{D}_{S}$, where $\mathcal{D}_{s}$ contains the nodes of the $s$-th sensitive group, we define the $\mathrm{Balance}$ of cluster $\mathcal{C}_{k}$ as [8]

where $[S]:=\{1,...,S\}$. The higher the $\mathrm{Balance}$ of each cluster, the fairer the clustering [8]. It is not difficult to check that $\underset{k\in[K]}{\min}\mathrm{Balance}(\mathcal{C}_{k})\leq\underset{s\neq s^{\prime}\in[S]}{\min}\;|\mathcal{D}_{s}|/|\mathcal{D}_{s^{\prime}}|$. Thus, this notion of fairness is asking for a clustering where the fraction of different sensitive groups in each cluster is approximately the same as that of the entire dataset $\mathcal{V}$ [14], which is also called group fairness. To incorporate this fairness notion into SC, a group-membership vector $\mathbf{f}_{s}\in\{0,1\}^{D}$ of $\mathcal{D}_{s}$ is defined, where $(\mathbf{f}_{s})_{[i]}=1$ if $i\in\mathcal{D}_{s}$ and $(\mathbf{f}_{s})_{[i]}=0$ otherwise, for $s\in[S]$ and $i\in[D]$. Then, we have the following lemma.

Lemma 1 states that the proportional representation of all sensitive attribute samples in each cluster can be interpreted as a linear constraint $\mathbf{F}^{\top}\widetilde{\mathbf{U}}=\mathbf{0}$. Under this fairness constraint, unnormalized SC is equivalent to the following problem

Similarly, we can relax (8) to

Following traditional SC, existing FSC models perform $k$means on rows of $\mathbf{U}$ to obtain discrete cluster labels $\mathbf{Q}$.

Spectral rotation [19] is an alternative to $k$means for obtaining discrete clustering results from continuous labels $\mathbf{U}$, which is formulated as

where the set $\mathcal{I}$ contains all discrete cluster indicator matrices, and $\mathbf{R}\in\mathbb{R}^{K\times K}$ is an orthonormal matrix. According to the spectral solution invariance property [19], if $\mathbf{U}$ is a solution of (5), $\mathbf{U}\mathbf{R}$ is another solution. A suitable $\mathbf{R}$ could facilitate $\mathbf{U}\mathbf{R}$ as close to $\mathbf{Q}$ as possible. In contrast, $k$means is performed directly on $\mathbf{U}$ obtained from spectral embedding, which may far deviate from the real discrete results. Thus, spectral rotation usually achieves superior performance than $k$means [19].

We first introduce a variant of the stochastic block model [39] (vSBM) to generate random graphs with cluster structures and sensitive attributes [14]. This model assumes that there are two or more meaningful ground-truth clusterings of the observed data, and only one of them is fair. Assume that $\mathcal{V}$ comprises $S$ sensitive groups and is partitioned into $K$ clusters such that $|\mathcal{D}_{s}\cap\mathcal{C}_{k}|/|\mathcal{C}_{k}|=\zeta_{s},s\in[S],k\in[K]$, for $\zeta_{s}\in(0,1)$ with $\sum_{s=1}^{S}\zeta_{s}=1$. Based on the clusters and sensitive groups, we construct a random graph by connecting two vertices $i$ and $j$ with a probability $\mathrm{Pr}(i,j)$ that depends on the clusters and sensitive groups of $i$ and $j$. We define

where $\pi_{C}:[D]\to[K]$ and $\pi_{S}:[D]\to[S]$ are two functions that assign a node $i\in\mathcal{V}$ to one of the clusters and sensitive groups, respectively. Let $\mathbf{L}^{*}$ be the real graph Laplacian matrix generated by the vSBM method and $\widehat{\mathbf{L}}$ be the Laplacian matrix estimated by any graph construction method. The matrix $\widehat{\mathbf{L}}$ is used as the input to fair spectral embedding in (9), and spectral rotation is utilized to obtain discrete cluster labels. Our goal is to derive a fair clustering error bound related to the estimation error between $\widehat{\mathbf{L}}$ and $\mathbf{L}^{*}$. Let us make some assumptions.

Assumption 1 is similar to the $(1+\epsilon)-$approximation of $k$means [40], which provides the estimation accuracy of spectral rotation. Assumption 2 is the same as that in Theorem 1 of [14], which is only made to facilitate theoretical analysis. In practice, Assumption 2 may be violated, which, however, does not affect the effectiveness of FSC algorithms [14]. Based on the two assumptions, we have the following proposition.

According to [14], the meaning of “the number of misclassified vertices is at most $D_{m}$” is that there exists a permutation of cluster indices such that the clustering results up to this permutation successfully predict all cluster labels but $D_{m}$ many vertices. Note that the error bound consists of two parts. The first one is caused by the difference between the expected and real graph produced by the vSBM method, which is similar to [14]. The second part is related to the estimation error of graph construction methods. The fair constraint affects clustering performance via $\mathbf{Z}$, which is a matrix determined by sensitive group-membership matrix $\mathbf{F}$. For convenience, $\mathbf{Z}^{\top}\mathbf{L}\mathbf{Z}$ is dubbed fair graph. Generally, the error bound in (15) depends on $K$, $D$ and $\epsilon$. If we divide (15) by $D$, we obtain the bound for the misclassification rate. The first part of the misclassification rate bound tends to zero as $D$ goes to infinity, meaning that if $\mathbf{L}^{*}$ is exactly estimated (the second part equals to zero), performing FSC via (9) and spectral rotation is weakly consistent [14]. However, $\mathbf{L}^{*}$ usually cannot be estimated exactly, causing an additional error for subsequent fair clustering results. If the fair graph estimation error does not increase quadratically as $D$, the second part of the misclassification rate bound will also decay to zero. Proposition 1 illustrates that a well-estimated graph $\widehat{\mathbf{L}}$, which is close to $\mathbf{L}^{*}$, brings a small misclassification error bound. Thus, it motivates us to seek a more effective method to construct accurate graphs from observed data.

Given $N$ observed data $\mathbf{X}_{o}\in\mathbb{R}^{D\times N}$, we need to infer the underlying similarity graph topology as the input to FSC algorithms. However, contaminated data may lead to poor graph estimation performance, as indicated in Proposition 1, which degrades subsequent fair clustering performance. Therefore, we propose a method to learn graphs from potentially noisy data $\mathbf{X}_{o}$, which is formulated as

where $\mathcal{L}:=\left\{\mathbf{L}:\mathbf{L}\in\mathbb{S}^{D\times D},\,\mathbf{L}\mathbf{1}=\mathbf{0},\,\mathbf{L}_{[ij]}\leq 0,\,\,i\neq j\right\}$ contains all Laplacian matrices. Moreover, $\xi$ and $\beta$ are parameters, and $\bm{\upsilon}\in\mathbb{R}^{D}$ is a vector of adaptive weights. We let $\mathbf{\Upsilon}:=\mathrm{diag}(\sqrt{\bm{\upsilon}})$, where $\sqrt{\bm{\upsilon}}=(\sqrt{\bm{\upsilon}_{[1]}},...,\sqrt{\bm{\upsilon}_{[D]}})^{\top}$. Eq.(16) is a joint model of denoising and smoothness-based GL [30].

1) Denoising: If $\mathbf{L}$ is fixed, the problem (16) becomes

The model is a node-adaptive graph filter, and $\bm{\upsilon}$ represents node weights. Specifically, given node weights $\bm{\upsilon}$, we have

Taking the derivative of (18) and setting it to zero, we obtain

We let $\mathbf{K}:=\left(\mathbf{I}+{\xi}(\mathbf{\Upsilon}^{\top}\mathbf{\Upsilon})^{-1}\mathbf{L}\right)^{-1}$, which is positive definite and has eigen-decomposition $\mathbf{K}=\mathbf{\Theta}\mathbf{\Lambda}\mathbf{\Theta}^{\top}$ with eigenvalues matrix $\mathbf{\Lambda}$ and eigenvectors matrix $\mathbf{\Theta}$. Moreover, $\mathbf{\Lambda}=\mathrm{diag}\left(\frac{1}{1+\xi\lambda_{1}},....,\frac{1}{1+\xi\lambda_{D}}\right)$, where $0=\lambda_{1}\leq,....,\leq\lambda_{D}$ are the eigenvalues of $(\mathbf{\Upsilon}^{\top}\mathbf{\Upsilon})^{-1}\mathbf{L}$. From the perspective of graph spectral filtering (GFT) [28], $\mathbf{K}\mathbf{X}_{o}=\mathbf{\Theta}\mathbf{\Lambda}\mathbf{\Theta}^{\top}\mathbf{X}_{o}$ can be interpreted as that the observed graph signals (columns of $\mathbf{X}_{o}$) are first transformed to the graph frequency domain via $\mathbf{\Theta}^{\top}$, attenuated GFT coefficients according to $\mathbf{\Lambda}$, and transformed back to the nodal domain via $\mathbf{\Theta}$. It is observed from $\mathbf{\Lambda}$ that the graph filter $\mathbf{K}$ is low-pass since the attenuation is stronger for larger eigenvalues. Thus, the graph filter can suppress the high-frequency component of raw data $\mathbf{X}_{o}$, alleviating the noises on the graph to output the “noiseless” signal $\mathbf{X}$.

Our graph filter $\mathbf{K}$ differs from the Auto-Regressive graph filter $\left(\mathbf{I}+{\xi}\mathbf{L}\right)^{-1}$ [41] in that we assign each node an individual weight $\bm{\upsilon}_{[i]},i=1,...,D$. The reason for using $\bm{\upsilon}$ is that the measurement noise of different nodes may be heterogeneous. If the $i$-th node signal (the $i$-th row of $\mathbf{X}_{o}$) has a small noise scale, a large $\bm{\upsilon}_{[i]}$ should be assigned to the fidelity term of node $i$ in (17) to ensure $\mathbf{X}_{[i,:]}$ is are close to the corresponding observation $(\mathbf{X}_{o})_{[i,:]}$ [42]. When we cannot know the noise scale a priori, we can adaptively learn $\bm{\upsilon}$ from the data. Specifically, given $\mathbf{X}$, the problem (17) becomes

Intuitively, solving (20) will assign a large $\bm{\upsilon}_{[i]}$ to node $i$ if $\mathbf{X}_{[i,:]}$ is close to $(\mathbf{X}_{o})_{[i,:]}$, as expected.

2) Graph learning: If we have obtained the “noiseless” signals $\mathbf{X}$ via the graph filter $\mathbf{K}$, the problem (16) becomes

The first Laplacian quadratic term of (21) is equivalent to

which measures the average smoothness of data $\mathbf{X}$ over the graph $\mathbf{L}$ [30]. The second term of (21) contains regularizers that endow the learned graphs with desired properties. The $\log$-degree term is to control node degree, and the Frobenius norm term is to control graph sparsity. Our model (21) can learn a graph suitable for graph-based clustering tasks for the following reasons. (i) It is observed from (22) that minimizing the smoothness is to seek a graph whose similar vertices (node signals) are closely connected, which is consistent with the fundamental principle of SC. (ii) The $\log$-degree term can avoid isolated nodes, which is crucial for SC, especially for normalized SC [5]. (iii) The Frobenius norm term can lead to a sparse graph, which may remove redundant and noisy edges.

The model (21) is similar to the ANGL method [26] since both construct graphs by minimizing the smoothness. The main differences lie in three aspects. (i) Our model removes the sum-to-one constraint in the ANGL method—the degree of each node is forced to be one—since the constraint makes the output graphs sensitive to noisy points [27]. Removing this constraint allows our model to capture more complex similarity relationships. (ii) We add a $\log$-degree term to ensure the learned graph has no isolated nodes. (iii) The input data of (21) are those produced by a low-pass graph filter.

3) Discussion: We try to explain why our method is effective in constructing graphs from observed data. If data $\mathbf{X}_{o}$ have a clustering structure, they should follow the assumption of cluster and manifold, i.e., the data in the same cluster are close to each other. According to [43], smooth signals containing low-frequency parts tend to follow the cluster and manifold assumption. Thus, if $\mathbf{L}$ accurately represents the graph behind observed data, the denoising part of our model has two functions. First, it filters out the high-frequency components of the observed graph signals that correspond to noises. Second, it produces smooth signals that have a clearer clustering structure, which could facilitate subsequent clustering. To better illustrate the effectiveness of the node-adaptive filter, Fig.2 depicts the t-SNE [44] results of our methods on the MNIST dataset, where four clusters correspond to four randomly selected digits. We can see that raw data are entangled together. In contrast, the denoised data $\mathbf{X}$ by the graph filter are clearly separated, meaning that the denoising part of our model can produce cluster-friendly signals. From the perspective of GL, our model (21) learns a graph minimizing the smoothness of data, i.e., the nodes corresponding to similar signals are closely connected. Thus, the learned graph can effectively capture similarity relationships between data and preserve clustering structures. Consequently, the denoising operation and the smoothness-based GL can enhance each other collaboratively to bring a high-quality graph for subsequent fair clustering tasks.

In this subsection, we build an end-to-end FSC framework that inputs observed data $\mathbf{X}_{o}$ and node attributes $\mathbf{F}$ and directly outputs discrete cluster labels. As shown in Fig.3, our model consists of four modules, i.e., denoising, graph learning, fair spectral embedding, and discretization. First, we construct graphs from the observed data by the proposed method (16). Once we obtain $\mathbf{L}$, the Laplacian matrix together with $\mathbf{F}$ can be directly used to perform fair spectral embedding (9) to obtain continuous clustering label matrix $\mathbf{U}$. Finally, we leverage spectral rotation (10) instead of $k$means to obtain discrete cluster labels. In addition to the reason for superior performance as stated in Section III-B, we utilize spectral rotation since it can be flexibly integrated into an end-to-end framework. Integrating all the above subtasks into a single objective function, we obtain

where $\mu$ and $\gamma$ are two parameters. All modules are not simply put together, but are bridged by two Laplacian quadratic terms, i.e., $\frac{\xi}{N}\mathrm{Tr}(\mathbf{X}^{\top}\mathbf{L}\mathbf{X})$ and $\mu\mathrm{Tr}(\mathbf{U}^{\top}\mathbf{L}\mathbf{U})$. First, $\frac{\xi}{N}\mathrm{Tr}(\mathbf{X}^{\top}\mathbf{L}\mathbf{X})$ can be viewed as the graph Tikhonov regularizer of the denoising task (18) to output smooth signals [41]. On the other hand, it measures smoothness in the GL task to capture the similarity relationships between data. Second, $\mu\mathrm{Tr}(\mathbf{U}^{\top}\mathbf{L}\mathbf{U})$ together with $\mathbf{F}^{\top}\mathbf{U}=\mathbf{0}$ performs fair spectral embedding as input to the discretization. It is also used to impose structural constraints on the constructed graph, which is discussed in the next subsection. The four modules are coupled with each other to achieve the overall optimal results for all subtasks.

To better understand how the fairness constraint works, we introduce a new variable matrix $\mathbf{Y}\in\mathbb{R}^{(D-S+1)\times K}$ and let $\mathbf{U}=\mathbf{Z}\mathbf{Y}$, where $\mathbf{Z}$ is the matrix defined in Proposition 1. The matrix $\mathbf{F}$ encodes sensitive information, as does $\mathbf{Z}$. Then, problem (23) can be rephrased in term of $\mathbf{Y}$ as

In (24), the fairness constraint $\mathbf{F}^{\top}\mathbf{{U}}=\mathbf{0}$ is removed. We conduct spectral embedding on the fair graph $\mathbf{Z}^{\top}\mathbf{L}\mathbf{Z}$, which encodes graph topology and sensitive information simultaneously, instead of $\mathbf{L}$ to ensure fair clustering. The impact of fairness constraints is discussed in the next subsection.

The basic formulation (23) is flexible and has many possible extensions. Here are some examples. (i) We can replace spectral rotation with improved spectral rotation [45] to further improve discretization performance. (ii) We can introduce self-weighted features into (23) to determine the importance of different features in assigning cluster labels [46]. (iii) We can extend (23) from unnormalized SC (5) to normalized SC [47]. (iv) We can also incorporate individual fairness [17] into our unified model. We place the details of these extensions in the supplementary material for completeness.

1) Connections to community-based GL models: If we only focus on GL and fair spectral embedding, Eq.(24) becomes

where $\mathbf{X}$ is regarded as the “noiseless” data here. According to Ky Fan’s theorem [48], we have $\underset{\mathbf{Y}^{\top}\mathbf{Y}=\mathbf{I}}{\min}\;\mathrm{Tr}(\mathbf{Y}^{\top}\mathbf{Z}^{\top}\mathbf{L}\mathbf{Z}\mathbf{Y})=\sum_{k=1}^{K}\widetilde{\lambda}_{k}$, where $\widetilde{\lambda}_{k}$ is the $k$ smallest eigenvalue of $\mathbf{Z}^{\top}\mathbf{L}\mathbf{Z}$. Thus, the problem (25) can be rephrased as

Note that $\mathbf{Z}^{\top}\mathbf{L}\mathbf{Z}$ is a semi-positive definite matrix, i.e., $\widetilde{\lambda}_{k}\geq 0$. Minimizing (26) is equivalent to forcing $\sum_{k=1}^{K}\widetilde{\lambda}_{k}\to 0$ if $\mu$ is large enough. That is, (26) encourages the fair graph to have $K$ connected components. Therefore, (26) can be viewed as the community-based GL model, which has been widely studied. For example, [49] lets the first $K$ smallest eigenvalues of Laplacian matrices be zero to obtain the community structures, which can be relaxed to the last term in (26). The works [50, 26] force the rank of Laplacian matrices to $D-K$, which can also be interpreted as minimizing the sum of the first $K$ eigenvalues. Furthermore, [51] adds a term $\mathrm{Tr}(\mathbf{\Xi}^{\top}\mathbf{L}\mathbf{\Xi})$ to impose community constraints, where $\mathbf{\Xi}\mathbf{\Xi}^{\top}$ contains the value 1 for within-community edges only and 0 everywhere else. Although closely related, (25) differs from existing community-based GL models in two key aspects. First, the basic GL models are different. Our model is based on the smoothness-based GL, while [49, 51] are based on statistical GL models like Graphical Lasso. Besides, [50, 26] are based on the ANGL method. Second, our model imposes the community constraint on the fair graph $\mathbf{Z}^{\top}\mathbf{L}\mathbf{Z}$ instead of on $\mathbf{L}$ like existing works. Thus, the fairness constraint may affect the topology of the learned graph to obtain fair clustering. We will test the impact of the fairness constraint in the experimental section.

2) Connections to unified SC models: If we remove the denoising module and fairness constraint, our model becomes

Again, $\mathbf{X}$ is treated as the observed data. The model (27) is an end-to-end SC model. Here, we discuss the connections between our model and those unified SC models integrating graph construction, spectral embedding, and discretization. As stated in Remark 1, we focus on basic formulations without additional extensions. The first model we compare is [20]

where $\alpha_{U},\mu_{U}$, and $\gamma_{U}$ are constant parameters. Moreover, $\mathcal{W}=\left\{\mathbf{W}:\mathbf{W}\in\mathbb{S}^{D\times D},\mathbf{W}\geq 0,\mathrm{diag}(\mathbf{W})=\mathbf{0}\right\}$ is the set containing all adjacency matrices. This is a unified SC model that leverages the sparse representation method [24] to construct graphs, which is different from our GL method.

Another unified SC model [23, 36] is concluded as

where $\beta_{J},\mu_{J}$, and $\gamma_{J}$ are constant parameters. The graph construction method in (29) is the ANGL method [26]. We have discussed the difference between our graph construction method and the ANGL in the previous subsection.

In summary, our model (27) differs from the existing unified SC models (28)-(29) mainly in the graph construction method. As Proposition 1 states, an accurate GL method can boost fair clustering performance. In the experimental section, we develop fair versions of (28)-(29) and compare them with our model (23) to illustrate the superiority of our model.

Our algorithm alternately updates $\mathbf{L},\mathbf{U},\mathbf{R}$, $\mathbf{Q}$, $\mathbf{X}$, and $\bm{\upsilon}$ in (23), i.e., updating one with the others fixed. For clarity, we omit the iteration index here. The following derivations are the updates in one iteration.

1) Update $\mathbf{L}$: The sub-problem of updating $\mathbf{L}$ is

The problem can be rewritten in terms of $\mathbf{W}$

where

and $Reg_{W}(\mathbf{W})=-\mathbf{1}^{\top}\mathrm{log}(\mathbf{W}\mathbf{1})+{\beta}\lVert\mathbf{W}\rVert_{\mathrm{F}}^{2}$. By the definition of $\mathcal{W}$, the free variables of $\mathbf{W}$ are the upper triangle elements. Thus, we define a vector $\mathbf{w}\in\mathbb{R}^{P},P:=\frac{D(D-1)}{2}$, satisfying that $\mathbf{w}=\mathrm{Triu}(\mathbf{W})$, where $\mathrm{Triu}(\cdot):\mathbb{R}^{D\times D}\to\mathbb{R}^{P}$ is a function that converts the upper triangular elements of a matrix into a vector. Then, the problem (31) is equivalent to

where $\mathbf{p}=\mathrm{Triu}(\mathbf{P}$), $\mathbf{S}\in\mathbb{R}^{D\times P}$ is a linear operator satisfying $\mathbf{S}\mathbf{w}=\mathbf{W}\mathbf{1}$ [30]. The problem (33) is convex, and we employ the algorithm in [52] to solve the problem. The complete algorithm flow is presented in the supplementary materials. After obtaining the estimated ${\mathbf{w}}$, we let ${\mathbf{W}}=\mathrm{iTriu}({\mathbf{w}})$, where $\mathrm{iTriu}(\cdot):\mathbb{R}^{P}\to\mathbb{R}^{D\times D}$ is the inverse $\mathrm{Triu}$ operation. The operation $\mathrm{iTriu}({\mathbf{w}})$ converts ${\mathbf{w}}$ into an adjacency matrix, where ${\mathbf{w}}$ corresponds to the upper triangle elements of ${\mathbf{W}}$. Finally, we calculate the Laplacian matrix from ${\mathbf{W}}$ and feed it into subsequent updates of other variables.

2) Update $\mathbf{U}$: The sub-problem of updating $\mathbf{U}$ is

Like (24), (34) can be cast into a problem of variable $\mathbf{Y}$

This is a typical quadratic optimization problem with orthogonal constraints. Let $\phi(\mathbf{Y})$ be the objective function of (35). We have that $\phi(\mathbf{Y})$ is differential, and $\nabla_{\mathbf{Y}}\,\phi(\mathbf{Y})=2\mu\mathbf{Z}^{\top}\mathbf{L}\mathbf{Z}\mathbf{Y}-2\gamma\mathbf{Z}^{\top}\mathbf{Q}\mathbf{R}^{\top}$. Thus, the problem can be efficiently solved via the algorithm in [53]. After obtaining $\mathbf{Y}$, we let $\textbf{U}=\textbf{Z}\textbf{Y}$.

3) Update $\mathbf{R}$: The sub-problem of updating $\mathbf{R}$ is

It is the orthogonal Procrustes problem with a closed-form solution [54]. Assuming that $\mathbf{\Theta}_{L}$ and $\mathbf{\Theta}_{R}$ are the left and right matrices of SVD of $\mathbf{Q}^{\top}\mathbf{U}$, the solution to (36) is [54]

4) Update $\mathbf{Q}$: The sub-problem of updating $\mathbf{Q}$ is

The optimal solution to (38) is as follows:

5) Update $\mathbf{X}$: The sub-problem of updating $\mathbf{X}$ is (17), which has a closed solution (19). However, matrix inversion is computationally expensive with complexity $\mathcal{O}(D^{3})$. Fortunately, $\left(\mathbf{\Upsilon}^{\top}\mathbf{\Upsilon}+{\xi}\mathbf{L}\right)^{-1}$ is symmetric, sparse, and positive definite. We can hence solve (17) efficiently using conjugate gradient (CG) algorithm without matrix inverse [55].

6) Update $\bm{\upsilon}$: The sub-problem of updating $\bm{\upsilon}$ is (20). Taking the derivative of (20) and setting it to zero, we have

It is observed that the updates of $\mathbf{L}$, $\mathbf{U}$, $\mathbf{R}$, $\mathbf{Q}$, $\mathbf{X}$, and $\bm{\upsilon}$ are coupled with each other. Updating one variable depends on the other variables, leading to an overall optimal solution. The complete procedure is presented in Algorithm 1.

1) Convergence analysis: It is challenging to obtain a globally optimal solution to (23) since it is not jointly convex for all variables. However, our algorithm for solving each sub-problem can reach its optimal solution. Specifically, when we update $\mathbf{L}$, the problem (32) is convex, and the corresponding algorithm is guaranteed to converge to the global optimum [52]. When updating $\mathbf{U}$, we use the algorithm in [53] to solve the problem (35), which can converge to the global optimum [53]. The updates of $\mathbf{Q}$, $\mathbf{R}$ and $\bm{\upsilon}$ have closed-form solutions. Despite updating $\mathbf{X}$ via (18) has a closed solution, we update $\mathbf{X}$ using CG, which is guaranteed to converge [55]. In summary, the update of each variable converges in our algorithm. In reality, the whole algorithm converges well, which is verified experimentally in Section VI.

2) Complexity analysis: In one iteration, our algorithm consists of six parts, which we analyze one by one below. As stated in [52], the update of $\mathbf{L}$ requires $\mathcal{O}(T_{1}D^{2})$ costs, where $T_{1}$ is the average number of iterations of updating $\mathbf{w}$. The computational cost can be further reduced if the average number of neighbors per node is fixed; see [56] and analysis therein. The computational complexity of our algorithm for updating $\mathbf{U}$ is $\mathcal{O}(T_{2}(DK^{2}+K^{3}))$ according to [53], where $T_{2}$ is the average number of iterations of the algorithm in [53]. When updating $\mathbf{R}$, we perform SVD on $\mathbf{Q}^{\top}\mathbf{U}\in\mathbb{R}^{K\times K}$, which costs $\mathcal{O}(K^{3})$. The updates of $\mathbf{Q}$ and $\bm{\upsilon}$ require $\mathcal{O}(DK^{2})$ and $\mathcal{O}(DN)$, respectively. Finally, the complexity of using CG to update $\mathbf{X}$ is $\mathcal{O}(T_{3}DN)$, where $T_{3}$ is the average number of iterations of the CG algorithm.

1) Graph generation: For synthetic data, we leverage the vSBM method to generate random graphs with sensitive attributes. Specifically, we let $\zeta_{s}=\frac{1}{S},a=0.8,b=0.2,c=0.15$, and $d=0.05$. After obtaining connections among nodes, we assign each edge a random weight between $[0.1,2]$. Finally, we normalize the edge weights to satisfy $\mathrm{Tr}(\mathbf{L}^{*})=D$.

2) Signal generation: We generate $N$ observed signals from the following Gaussian distribution [29]

where $\mathbf{\Sigma}_{e}=\mathrm{diag}(\sigma_{1}^{2},...,\sigma_{D}^{2})$ and $\sigma_{i}$ is the noise scale of the $i$-th node. As stated in [29], signals generated in this way are smooth over the corresponding graph.

3) Evaluation metrics: In topology inference, determining whether two vertices are connected can be regarded as a binary classification problem. Thus, we employ the F1-score ($\mathrm{FS}$) to evaluate classification results

where $\mathrm{TP}$ is true positive rate, $\mathrm{TN}$ is true negative rate, $\mathrm{FP}$ is false positive rate, and $\mathrm{FN}$ is false negative rate. We also use the estimation error (${\mathrm{EE}}$) to evaluate the learned graph

For a fair comparison of $\mathrm{EE}$, we normalize the learned graphs to $\mathrm{Tr}(\widehat{\mathbf{L}})=D$. For fair clustering, we use the same two metrics as in [14]: clustering error ($\mathrm{CE}$) and Balance ($\mathrm{Bal}$)

where $\widehat{\pi}_{C}(i)$ is the estimated cluster label of node $i$ (after proper permutation), and ${\pi}_{C}(i)$ is the ground-truth. The metric $\mathrm{Balance}$ measures the average balance of all clusters.

4) Baselines: The comparison baselines are list in Table.I. The model Fairlets is the fair version of $k$median [8]. Models 3-5 are the implementations of [14] using different graph construction methods. FGLASSO [57] is the only model that jointly performs graph construction and fair spectral embedding. FSRSC and FJGSED are the fair versions of unified SC models (28) and (29). The formulations and algorithms for the two models are placed in the supplementary material.