Bounds on Perfect Node Classification: A Convex Graph Clustering Perspective

TL;DR

This paper introduces a new convex optimization approach for perfect node classification in graphs, leveraging community structure and node-specific data to improve accuracy under milder conditions.

Contribution

It proposes a novel spectral graph clustering framework that integrates node labels and features, providing theoretical guarantees for perfect community recovery.

Findings

Guarantees perfect community recovery under milder conditions

Demonstrates synergy between graph structure and node-specific information

Validates approach with numerical experiments

Abstract

We present an analysis of the transductive node classification problem, where the underlying graph consists of communities that agree with the node labels and node features. For node classification, we propose a novel optimization problem that incorporates the node-specific information (labels and features) in a spectral graph clustering framework. Studying this problem, we demonstrate a synergy between the graph structure and node-specific information. In particular, we show that suitable node-specific information guarantees the solution of our optimization problem perfectly recovering the communities, under milder conditions than the bounds on graph clustering alone. We present algorithmic solutions to our optimization problem and numerical experiments that confirm such a synergy.

Equations

$$\underset{\mathbf{L}\in \mathcal{B}}{\min }\Vert \mathbf{L}-\mathbf{A}{\Vert }_1+{\mu }_0\Vert \mathbf{L}{\Vert }_{\ast },$$

$$\underset{\mathbf{L}\in \mathcal{B}}{\min }\Vert \mathbf{L}-\mathbf{A}{\Vert }_1+{\mu }_0\Vert \mathbf{L}{\Vert }_{\mathcal{A}}.$$

$$\underset{\begin{array}{c}\mathbf{L}\in \mathcal{B},\{{\lambda }_i\ge 0,{\mathbf{e}}_i\in {\mathbb{S}}^{n-1}\end{array}{\}}_i}{\min }$$

$$\mathbf{L}=\sum _i{\lambda }_i{\mathbf{e}}_i{\mathbf{e}}_i^{\top },$$

$$\underset{\begin{array}{c}\mathbf{L}\in \mathcal{B},\{{\mathbf{\theta }}_v\in \Theta {\}}_{v\in \mathcal{V}},\\ \{{\lambda }_i\ge 0,{\mathbf{e}}_i\in {\mathbb{S}}^{n-1},{\mathbf{\theta }}_i\in \Theta {\}}_i\end{array}}{\min }$$

$$\mathbf{L}=\sum _i{\lambda }_i{\mathbf{e}}_i{\mathbf{e}}_i^{\top },{\mathbf{\theta }}_v=\sum _i{\lambda }_i{ϵ}_{i,v}^2{\mathbf{\theta }}_i,\forall v\in \mathcal{V}$$

$$\mathcal{A}:=\Big{\{}\big{(}\mathbf{e}\mathbf{e}^{\top},\{\epsilon_{v}^{2}\bm{\theta}\}_{v\in\mathcal{V}}\big{)}\;\Big{|}\;\mathbf{e}=(\epsilon_{v})_{v},\;\|\mathbf{e}\|_{2}=1,\;\bm{\theta}\in\Theta\Big{\}}.$$

$$\underset{\begin{array}{c}\mathbf{U}\in \mathcal{U},\{{\lambda }_i\ge 0,{\mathbf{a}}_i{\}}_i\end{array}}{\min }$$

$$\mathbf{U}=\sum _i{\lambda }_i{\mathbf{a}}_i,$$

$$\underset{\begin{array}{c}\mathbf{U}\in \mathcal{U},\{{\lambda }_i\ge 0,{\mathbf{a}}_i{\}}_i\end{array}}{\min }$$

$$\mathbf{U}=\sum _i{\lambda }_i{\mathbf{a}}_i$$

$$n_{v,j}:=\sum _{u\in {\mathcal{C}}_j}{A}_{uv}.$$

$$N_{i,j}^{+}=\sum _{v\in {\mathcal{C}}_i}{n}_{v,j}^{+},{\rho }_{i,j}^{+}=\frac{N_{i,j}^{+}}{n_i{n}_j}.$$

$$\left|\frac{n_{v,j}^{+}}{n_j}-{\rho }_{ij}^{+}\right|\le \delta .$$

$$\underset{{\hat{\mathbf{\theta }}}_i\in \Theta }{\min }\mu \sum _{i=1}^K{n}_i{F}_i({\hat{\theta }}_i)+{\mu }_1\sum _{i<j}{n}_i{n}_j\Vert {\hat{\mathbf{\theta }}}_i-{\hat{\mathbf{\theta }}}_j\Vert ,$$

$$\nabla F_i({\hat{\mathbf{\theta }}}_i)+\gamma \sum _{j\ne i}\frac{n_j}{n}\frac{{\hat{\mathbf{\theta }}}_i-{\hat{\mathbf{\theta }}}_j}{\Vert {\hat{\mathbf{\theta }}}_i-{\hat{\mathbf{\theta }}}_j\Vert }:={\mathbf{n}}_i\in -\partial I_{\Theta }({\hat{\mathbf{\theta }}}_i),$$

$$\mu \sim \text{nicefrac}n{n}_i\text{if }{n}_i\le \sqrt{n},\mu \sim \text{nicefrac}n\sqrt{n_i}\text{if }{n}_i\ge \sqrt{n},$$

$$\rho \le \text{nicefrac}{n}_i^{2}n\text{if }{n}_i\le \sqrt{n},\rho \le c\gamma \cdot \text{nicefrac}{n}_{i}n\text{if }{n}_i\ge \sqrt{n}.$$

$${\mathbf{n}}^{\ast }\in \partial I_{\mathcal{U}}({\mathbf{U}}^{\ast }),\underset{\mathbf{a}\in \mathcal{A}}{\sup }{\mathbf{a}}^{\top }{\mathbf{z}}^{\ast }\le {\mu }_0,{\mathbf{g}}^{\ast }\in \partial \varphi ({\mathbf{U}}^{\ast }),{\mathbf{g}}^{\ast }+{\mathbf{z}}^{\ast }+{\mathbf{n}}^{\ast }=0.$$

$${\tilde{\mathbf{\theta }}}_i=\arg \underset{{\bar{\mathbf{\theta }}}_i\in \Theta }{\min }$$

$$\tilde{\mathbf{E}}=\arg \underset{\{{\bar{\mathbf{e}}}_i{\}}_{i=1}^r}{\min }$$

$$\mathbf{L}=\sum _{i=1}^r{\bar{\mathbf{e}}}_i{\bar{\mathbf{e}}}_i^{\top }\in \mathcal{B},{\mathbf{\theta }}_v=\sum _{i=1}^r{\overline{ϵ}}_{i,v}^2{\bar{\mathbf{\theta }}}_i\in \Theta ,$$

$$f_{\mathrm{feature}}(\mathbf{x};\mathbf{R}):=\frac{1}{m}\left({\mathbf{x}}^{\top }{\mathbf{R}}^{-1}\mathbf{x}+\text{Tr}(\mathbf{R})\right),$$

$$f_{\mathrm{label}}(\tilde{y}=j;\mathbf{\pi }):=-{\pi }_j.$$

$$f_{v}(\mathbf{z};\bm{\theta})=\left\{\begin{array}[]{cc}f_{\mathrm{feature}}(\mathbf{x};\mathbf{R})+\beta f_{\mathrm{label}}(\tilde{y};\bm{\pi})&v\in\mathcal{V}_{\mathrm{train}}\\ f_{\mathrm{feature}}(\mathbf{x};\mathbf{R})&v\in\mathcal{V}_{\mathrm{test}}\end{array}\right.$$

$${\mathbb{S}}_{{\rho }_{-},{\rho }_{+}}=\left\{\mathbf{R}\mid \mathbf{R}={\mathbf{R}}^{\top },\forall \mathbf{x}\in {\mathbb{R}}^n:{\rho }_{-}\Vert \mathbf{x}{\Vert }_2^2\le {\mathbf{x}}^{\top }\mathbf{Rx}\le {\rho }_{+}\Vert \mathbf{x}{\Vert }_2^2\right\}.$$

$$\Vert \mathbf{x}{\Vert }_{\mathcal{A}}=\inf \left\{\sum _i{c}_i:\mathbf{x}=\sum _i{c}_i{}_i,{}_i\in \mathcal{A},c_i\ge 0\right\}.$$

$${\mathcal{A}}_{{\ell }_1}=\{\pm {\mathbf{e}}_i|i=1,\dots ,n\},$$

$$\Vert \mathbf{x}{\Vert }_{{\mathcal{A}}_{{\ell }_1}}=\sum _{i=1}^n|x_i|=\Vert \mathbf{x}{\Vert }_1,$$

$${\mathcal{A}}_{\ast }=\{\mathbf{u}{\mathbf{v}}^{\top }|\mathbf{u}\in {\mathbb{R}}^m,\mathbf{v}\in {\mathbb{R}}^n,\Vert \mathbf{u}{\Vert }_2=\Vert \mathbf{v}{\Vert }_2=1\}.$$

$$\Vert \mathbf{X}{\Vert }_{{\mathcal{A}}_{\ast }}=\sum _{i=1}^{\min (m,n)}{\sigma }_i(\mathbf{X})=\Vert \mathbf{X}{\Vert }_{\ast },$$

$${\mathcal{A}}_{\mathrm{sym},+}=\{\mathbf{u}{\mathbf{u}}^{\top }:\Vert \mathbf{u}{\Vert }_2=1\},{\mathcal{A}}_{\mathrm{sym},\pm }=\{\pm \mathbf{u}{\mathbf{u}}^{\top }:\Vert \mathbf{u}{\Vert }_2=1\}.$$

$$\Vert \mathbf{X}{\Vert }_{{\mathcal{A}}_{\mathrm{sym},+}}=\sum _i{\lambda }_i=\mathrm{tr}(\mathbf{X})\text{for }\mathbf{X}⪰0,$$

$$\Vert \mathbf{X}{\Vert }_{{\mathcal{A}}_{\mathrm{sym},\pm }}=\sum _i|{\lambda }_i|=\Vert \mathbf{X}{\Vert }_{\ast }.$$

Full text

1 Introduction

Consider a population of individuals $v\in\mathcal{V}$ and their given feature vectors $\{\mathbf{x}_{v}\}$. The goal of transductive classification is to assign class labels $\hat{y}_{v}$ to a group $\mathcal{V}_{\mathrm{test}}$ of test individuals in $\mathcal{V}$ by observing the noisy labels $\tilde{y}_{v}$ of the remaining group $\mathcal{V}_{\mathrm{train}}$. Node classification (Bhagat et al., 2011; Xiao et al., 2022) is a variation on this theme, where the population is endowed with an additional unweighted graph $\mathcal{G}=(\mathcal{V},E)$, with the edges $E$ representing some relationship within the population. One may naturally expect that under suitable conditions, incorporating $\mathcal{G}$ enhances classification. In this paper, we provide conditions, under which such a gain provably exists.

Our study assumes the homophily property (McPherson et al., 2001), which means that the individuals (nodes) belonging to the same class are more likely to be connected in the graph. Many real-world graph-structured datasets, such as citation networks (Cora (Sen et al., 2008), Citeseer (Sen et al., 2008), Pubmed (Namata et al., 2012)), social networks (Reddit Social Networks (Hamilton et al., 2017)), and commercial networks (Amazon Co-Purchase (McAuley et al., 2015)) exhibit homophily. In these datasets, the graph $\mathcal{G}$ possesses partitions that tend to align with node-specific information $\{\mathbf{x}_{v},\tilde{y}_{v}\}$. As a result, we observe a natural relation between node classification and graph clustering (Englert & Glantz, 1999; Schaeffer, 2007) methods, which identify such partitions, typically under the assumption that well-connected groups of nodes belong to the same class or cluster.

Under homophily, graph clustering and node classification are complementary and can benefit from each other when approached together. This motivates us to develop a framework that integrates graph clustering and classification into a unified optimization problem. Our approach is inspired by spectral graph clustering (Belkin & Niyogi, 2001; Ng et al., 2001; Hajek et al., 2016) and convex optimization techniques for graph clustering (Korlakai Vinayak et al., 2014; Chen et al., 2014; Li et al., 2021), all of which share a common approach: seeking a low-rank, positive semi-definite (PSD) approximation of the adjacency matrix $\mathbf{A}$ of $\mathcal{G}$ (Korlakai Vinayak et al., 2014):

[TABLE]

where $\|\cdot\|_{1}$ is the $\ell_{1}$-norm (sum of absolute values of the entries of the matrix), $\|\cdot\|_{*}$ is the nuclear norm (sum of the singular values of the matrix), and $\mathcal{B}$ denotes the set of symmetric matrices with entries in the interval $[0,1]$, and with ones on the diagonal.

Considering that the nuclear norm is a specific instance of the atomic norm (Bhaskar et al., 2013) (see Appendix A for definitions and details), we may replace $\|\mathbf{L}\|_{*}$ with $\|\mathbf{L}\|_{\mathcal{A}}$, where the atomic set $\mathcal{A}$ contains rank-1 matrices $\mathbf{a}_{i}=\mathbf{e}_{i}\mathbf{e}_{i}^{\top}$ as atoms, i.e., $\mathcal{A}\coloneqq\{\mathbf{e}\mathbf{e}^{\top}:\|\mathbf{e}\|_{2}=1\}$. With this reformulation, Problem (1) becomes:

[TABLE]

Considering the atomic norm definition (Appendix A, equation (19)) and noting that, for optimization purposes, $\|\mathbf{L}-\mathbf{A}\|_{1}$ can be equivalently expressed as $-\langle\bar{\mathbf{A}},\mathbf{L}\rangle$,111This follows from the fact that for $a_{ij}\in\{0,1\}$ and $l_{ij}\in[0,1]$, we have $|l_{ij}-a_{ij}|=-l_{ij}\bar{a}_{ij}+a_{ij}$. one can equivalently express (1) as:

[TABLE]

where the elements of the polarized adjacency matrix $\bar{\mathbf{A}}$ are $1$ if $\{u,v\}\in E$ and $-1$ otherwise. The summations are over the infinite elements as the atomic set contains infinite number of atoms, and this hold for the rest of paper if the upper bound of the sum is not defined.

This optimization serves as the foundation of our framework, providing an opportunity to build upon it by incorporating node-specific information. It also highlights that if the graph has $r$ well-separated clusters, the optimal solution of (3) has rank $r$ and can be expressed as $\mathbf{L}=\mathbf{E}\bm{\Lambda}\mathbf{E}^{\top}$, where $\bm{\Lambda}=\mathrm{diag}(\lambda_{1},\ldots,\lambda_{r})$ and $\mathbf{E}=[\,\mathbf{e}_{1},\ldots,\mathbf{e}_{r}\,]$ is the matrix whose columns are $\{\mathbf{e}_{i}\}_{i}$. Each row of $\mathbf{E}$, denoted by $\bm{\epsilon}_{v}$, represents the embedding vector of node $v\in\mathcal{V}$ in the low-rank space. Under some conditions on the graph, these embeddings become one-hot vectors that perfectly indicate the cluster membership of each node, enabling exact cluster recovery.

Theoretical Contributions.

Firstly, we propose a novel optimization framework for transductive node classification by generalizing the convex graph clustering formulation in (3). Our framework integrates the graph structure with node-specific information, including features and partial labels, and is inspired by principles of multimodal learning. Secondly, we establish that our framework guarantees perfect label recovery, i.e., the exact recovery of the true labels for all nodes, under some structural assumptions. To the best of our knowledge, our results provide the first rigorous proof that node features and graph structure can have a synergistic effect in node classification.

Algorithmic Contributions.

We propose CADO, an efficient scalable alternating conditional gradient algorithm for solving our optimization framework with a fixed number of atoms. Each step of this algorithm reduces to tractable subproblems with closed-form solutions under our studied model. We validate our framework through experimental studies showing that CADO effectively solves our proposed problem. We also validate our theoretical findings through our experiments.222The code is available at this link.

Paper Organization.

Section 2 introduces our generalized optimization framework. Section 3 establishes the perfect recovery guarantees of our framework. In Section 4, we analyze the computational complexity of solving our framework, followed by an algorithmic solution. Section 5 reports our experimental results, corroborating our theoretical findings, and Section 6 discusses related work and concludes the paper.

2 Node Classification via Atomic Norms

In this section, we extend the convex graph clustering formulation in (3) to incorporate node-specific information. Let $\mathbf{z}_{v}$ denote the information associated with node $v\in\mathcal{V}$, where $\mathbf{z}_{v}=\mathbf{x}_{v}$ for test nodes ($v\in\mathcal{V}_{\mathrm{test}}$) and $\mathbf{z}_{v}=(\mathbf{x}_{v},\tilde{y}_{v})$ for training nodes ($v\in\mathcal{V}_{\mathrm{train}}$). We assume a generic set of models $\bm{\theta}\in\Theta$ and loss functions $f_{v}(\mathbf{z}_{v};\bm{\theta})$ that compute how well $\bm{\theta}$ fits $\mathbf{z}_{v}$, with a concrete example provided in Section 5.1. Based on these components, we formulate the following optimization problem:

[TABLE]

where $\epsilon_{i,v}$ is the $v^{\text{th}}$ element of $\mathbf{e}_{i}$, and the variable $\bm{\theta}_{i}\in\Theta$ is a specific model for each cluster. The additional term in the objective function of (4) is similar to an empirical risk function, but each data point $\mathbf{z}_{v}$ is instead examined against an individual model $\bm{\theta}_{v}$. However, the constraint $\bm{\theta}_{v}=\sum_{i}\lambda_{i}\epsilon^{2}_{i,v}\bm{\theta}_{i}$ does not allow these models to be completely independent. In particular, when the embedding vectors $\bm{\epsilon}_{v}=(\epsilon_{i,v})_{i}$ coincide with the one-hot vectors of the corresponding clusters, we have $\bm{\theta}_{v}=\bm{\theta}_{i}$ for all nodes $v$ in the $i^{\text{th}}$ class that we denote by $\mathcal{C}_{i}$. In this way, the classifier benefits from the low-rank clustering mechanism.

To express the optimization problem in (4) more compactly as a regularized atomic norm problem, similar to (2)–(3), we define the joint variable \mathbf{U}:=\big{(}\mathbf{L},\{\bm{\theta}_{v}\}_{v\in\mathcal{V}}\big{)}, the atoms \mathbf{a}_{i}:=\big{(}\mathbf{e}_{i}\mathbf{e}_{i}^{\top},\{\epsilon_{i,v}^{2}\bm{\theta}_{i}\}_{v\in\mathcal{V}}\big{)}, and the atomic set

[TABLE]

With this notation, (4) can be rewritten in the following compact form:

[TABLE]

where $\phi(\mathbf{U})$ denotes the objective function in (4), and $\mathcal{U}:=\mathcal{B}\times\Theta^{|\mathcal{V}|}$ represents the feasible set of all variables \big{(}\mathbf{L},\{\bm{\theta}_{v}\}_{v\in\mathcal{V}}\big{)} with $\mathbf{L}\in\mathcal{B}$ and $\bm{\theta}_{v}\in\Theta$.

2.1 Regularization by Sum of Norms

The formulation in (5) allows for an arbitrary number of atoms, each potentially corresponding to a distinct cluster. This flexibility can lead to over-parameterization, where the term in $\phi$ corresponding to node-specific information has an inherent tendency to over-parametrize and assign an individual cluster to each node. While the optimization terms corresponding to the graph may prevent over-parametrization, such a scenario hardly reflects a synergy between them and the node-specific information. In response, we utilize the well-known sum-of-norms (SON) (Hocking et al., 2011; Lindsten et al., 2011; Panahi et al., 2017; Tan & Witten, 2015; Sun et al., 2021) regularization of the node-specific information, given by the regularization function $\mathbf{R}(\mathbf{U}):=\sum_{u<v}\|\bm{\theta}_{u}-\bm{\theta}_{v}\|$. Accordingly, adding this regularization leads to the following variant of (5):

[TABLE]

This regularization enforces that if two nodes belong to the same cluster, then their corresponding models $\bm{\theta}_{u}$ and $\bm{\theta}_{v}$ are encouraged to be equal.

Equation (6) represents our final optimization framework, which builds upon the convex graph clustering formulation while incorporating node-specific information. In the following sections, we examine this framework from two perspectives. First, we analyze the conditions under which it achieves perfect recovery of the underlying clusters. Then, we study the computational complexity of solving the problem and introduce the CADO algorithm, which iteratively computes the solution.

3 Perfect Recovery

In this section, we provide conditions under which the global solution of (6) perfectly recovers the clusters and hence the class predictions. Our analysis considers a large population and adopts the standard planted partition model for the graph $\mathcal{G}$, which is well-known to be easily resolvable when the number of classes is fixed. In such a setting, the graph structure alone suffices for perfect recovery, and the node-specific information offers no additional benefit. To address this, we focus on a more challenging regime where the number of classes grows with the problem size, making the graph information alone insufficient and motivating the integration of node-specific information.

3.1 Setup

Under the planted partition model, the population $\mathcal{V}$ is assumed to be partitioned into clusters $\mathcal{C}_{1},\mathcal{C}_{2},\ldots,\mathcal{C}_{K}$, with $n_{i}$ individuals in cluster $\mathcal{C}_{i}$. Consistent with the homophily assumption, these clusters are reflected in both the graph structure $\mathcal{G}$ and the node-specific information $\mathbf{z}_{v}$, as discussed separately and in detail below.

3.1.1 Graph

Let $N_{ij}:=\sum_{u\in\mathcal{C}_{i},v\in\mathcal{C}_{j}}A_{uv}$ be the connectivity of different partitions. Naturally, we are interested in a case where $N_{ii}$ is sufficiently larger than $N_{ij}$ for $i\neq j$. Next, we let $n_{v,j}$ be the number of connections of node $v$ to nodes in cluster $\mathcal{C}_{j}$, i.e.

[TABLE]

We take $n^{+}_{v,j}=n_{j}-n_{v,j}$ if $v\in\mathcal{C}_{j}$, and $n^{+}_{v,j}=n_{v,j}$ if $v\notin\mathcal{C}_{j}$. We also define

[TABLE]

These quantities capture the level of misconnections in the graph and vanish in the ideal case of perfectly separated clusters. For a fixed real number $\delta>0$, we formalize and bound the level of cluster separability using the following assumptions.

Assumption 3.1** ($\delta-$Homogeneity).**

For any node $v\in\mathcal{C}_{i}$ and any cluster $\mathcal{C}_{j}$, it holds that

[TABLE]

Assumption 3.2** ($\delta-$Visibility).**

For any two clusters $\mathcal{C}_{i},\mathcal{C}_{j}$, it holds that $\rho^{+}_{ij}<\frac{1}{2}-\delta$.

Moreover, we define $\tilde{A}_{uv}:=A_{uv}-\nicefrac{{N_{ij}}}{{n_{i}n_{j}}}$ and refer to the matrix $\tilde{\mathbf{A}}=(\tilde{A}_{uv})$ as the centered adjacency matrix. For each cluster $\mathcal{C}_{i}$, the matrix $\tilde{\mathbf{A}}=(\tilde{A}_{uv})$ can be divided into three blocks: the intra-cluster block $\tilde{\mathbf{A}}_{i}$ corresponds to the edges with both ends in $\mathcal{C}_{i}$, the extra-cluster block $\tilde{\mathbf{A}}^{i}$ contains edges with no ends in $\mathcal{C}_{i}$, and the inter-cluster $\tilde{\mathbf{B}}_{i}\in\mathbb{R}^{n_{i}\times(n-n_{i})}$ block includes edges with exactly one end in $\mathcal{C}_{i}$. Accordingly, we consider the following assumption, where $\|\cdot\|_{\mathrm{op}}$ denotes the operator norm:

Assumption 3.3**.**

We consider a family of problems with a growing size $n$ and growing clusters $n_{i}$ such that for every $i$, we have $\|\tilde{\mathbf{A}}_{i}\|_{\mathrm{op}}=O(\sqrt{n_{i}})$, $\|\tilde{\mathbf{A}}^{i}\|_{\mathrm{op}}=O(\sqrt{n})$ and $\|\tilde{\mathbf{B}}_{i}\|_{\mathrm{op}}=O(\sqrt{n})$.

The centered adjacency matrix purely reflects the wrong edges and hence it implies a regime where the number of miss-connections is not too large. For instance, any random flip of edges with a constant probability satisfies this condition, no matter how many correct edges (i.e. edges for perfect clusters) may exist.

3.1.2 Node-Specific Vectors

In an ideal scenario, the optimization in (6) produces perfectly separated partitions, assigning a model $\hat{\bm{\theta}}_{i}$ to each cluster $\mathcal{C}_{i}$. These cluster-level models are associated with the following characteristic optimization problem:

[TABLE]

where $F_{i}(\bm{\theta}):=\frac{1}{n_{i}}\sum_{v\in\mathcal{C}_{i}}f_{v}(\mathbf{z};\bm{\theta})$ is the aggregate loss function of each model. The solutions of $\hat{\bm{\theta}}_{i}$ are called biased centroids. For simplicity, we define $\gamma:=\nicefrac{{n\mu_{1}}}{{\mu}}$ and then introduce the following assumption for a given constant $R$:

Assumption 3.4** ($R-$separability).**

The biased centroids $\hat{\theta}_{i}$ are distinct and the following optimality condition holds for each $i$:

[TABLE]

where $\partial I_{\mathcal{U}}(\cdot)$ denotes the subdifferential of the indicator function $I_{\mathcal{U}}$ of $\mathcal{U}$. Moreover, for every $\bm{\theta}\in\Theta$ the relation $\braket{\mathbf{n}_{i},\bm{\theta}-\hat{\bm{\theta}}_{i}}\leq R$ holds for at most one index $i$.

This separability assumption ensures that node-specific information is cluster-separable. Specifically, it guarantees that centroids corresponding to different clusters are well-separated. Another key feature of the node-specific vectors is the variability in their gradients. This feature is reflected in our next assumption which is based on another constant $\rho$:

Assumption 3.5** (Gradient Variability).**

For all $u,v\in\mathcal{C}_{i}$, the gradient variability is bounded: $\|\nabla f_{v}(\hat{\bm{\theta}}_{i})-\nabla f_{u}(\hat{\bm{\theta}}_{i})\|\leq\rho$.

This assumption bounds the variability of node-specific information within clusters. It ensures that within each cluster, the local losses behave similarly around the corresponding centroid $\hat{\bm{\theta}}_{i}$, each node aligns closely with its associated centroid. Combined with Assumption 3.4, this ensures that feature-based classification is reliable and interpretable.

3.2 Main result

Here, we present our main theoretical results. The proofs are presented in Appendix C.

Theorem 3.6**.**

Assume a sequence of problems where $n,K$, and all $n_{i}$s grow to infinity such that $n_{i}\sim n_{j}$ and all $\rho_{ij}^{+}$ converge to fixed nonzero values. Fix $\gamma=\nicefrac{{n\mu_{1}}}{{\mu}}$ and suppose that Assumptions 3.1–3.4 hold with fixed $\delta$ and $R$, and that the feasible set $\Theta$ is bounded. Then, the followings hold true:

Depending on $\gamma,\delta$, there exists a constant $c>1$ such that for

[TABLE]

the ideal solution with $\mathbf{e}_{i}$, $\lambda_{i}=n_{i}$, and $\bm{\theta}=\hat{\bm{\theta}}_{i}$ is optimal in (6), provided Assumption 3.5 holds and

[TABLE] 2. 2.

Without the graph term $-\braket{\bar{\mathbf{A}},\mathbf{L}}$, the ideal solution is optimal if $c=1$, i.e., Assumption 3.5 holds with $\rho\leq\gamma\cdot\nicefrac{{n_{i}}}{{n}}.$ 3. 3.

Without the node-specific terms, graph clustering recovers clusters if $n=O(n_{i}^{2})$.

This theorem clearly shows an improvement by the combination of the graph and the node-specific information. First, in the few-large-clusters regime (i.e., $n<n_{i}^{2}$), graph information alone is sufficient for exact recovery. However, by combining graph and node-specific information, our model can also recover clusters in the many-small-clusters regime (i.e., $n_{i}^{2}<n$), even when node features alone are insufficient. Moreover, Compared to part (2), part (1) allows for at least a factor-$c$ larger variation within the gradients in the few-large-clusters regime, and yields even greater improvements in the opposite regime.

4 Complexity Analysis and Algorithmic Solutions

In this section, we analyze the computational aspects of solving the optimization problem in (6). We first establish that there exists an iterative algorithm that achieves polynomial-time convergence to an $\epsilon$-optimal solution. Although this result provides theoretical guarantees, the corresponding algorithm is computationally expensive and does not scale well in practice. Motivated by this limitation, we next propose an efficient and scalable algorithm, referred to as CADO, which solves the non-convex formulation in (4) by assuming a fixed number of atoms. Intuitively, this fixed-rank approximation is well justified by the SON regularization, which naturally promotes sparsity in the number of active atoms, and in practice we observe that CADO reliably converges to the optimal solution.

4.1 Polynomial-Time Convergence

Our analysis builds on the concept of a linear minimization oracle (LMO), a standard tool in conditional gradient methods. The LMO for a dictionary $\mathcal{A}$ is defined as follows:

Definition 4.1** (LMO).**

The linear minimization oracle (LMO) of a dictionary $\mathcal{A}$ is a map $\mathcal{O}_{\mathcal{A}}(\mathbf{w})$ that, for any vector $\mathbf{w}$, returns an atom $\mathbf{a}\in\mathcal{A}$ minimizing the inner product $\mathbf{a}^{\top}\mathbf{w}$.

To establish our convergence result, we make the following assumption:

Assumption 4.2**.**

The feasible set $\mathcal{U}$ is closed and convex, the atomic set $\mathcal{A}$ is bounded, and the function $\phi(\mathbf{U})$ is convex and lower semi-continuous. Moreover, there exists an optimal solution $(\mathbf{U}^{*},\{\lambda_{i}^{*}\}_{i})$ of (5) and vectors $\mathbf{z}^{*},\mathbf{n}^{*},\mathbf{g}^{*}$ satisfying

[TABLE]

Under this assumption, we obtain the following result:

Theorem 4.3**.**

Under Assumption 4.2, there exists an algorithm that returns an $\epsilon$-approximate solution of (5) in \mathrm{poly}\big{(}\nicefrac{{1}}{{\epsilon}}\big{)} number of LMOs, proximal evaluations of $\phi$, and orthogonal projections onto $\mathcal{U}$.

The proof of Theorem 4.3 is provided in Appendix B. This result also extends to (6) by replacing $\phi$ with its composite form $\phi^{\prime}:=\phi+\mu_{1}R$. While this establishes the polynomial-time solvability of the problem, the construction is primarily of theoretical interest and is not computationally efficient for large-scale instances. In particular, solving (5) exactly in practice is challenging due to the additional feasibility constraint $\mathbf{U}\in\mathcal{U}$. Although algorithms such as CoGEnT (Rao et al., 2015) efficiently handle the unconstrained version of the problem, extending them to the constrained setting remains an open question. To address these practical challenges, we turn to a scalable approach that operates on the non-convex formulation (4) and provides reliable performance in practice.

4.2 Efficient Algorithmic Solution: CADO

To solve (6) with a fixed number of atoms, we propose an alternating conditional gradient algorithm, termed CADO (Constrained Atomic Decomposition Optimization). The algorithm alternates between updating the embedding vectors $\{\bar{\mathbf{e}}_{i}\}$, which define the low-rank matrix $\mathbf{L}$, and updating the cluster-level models $\{\bar{\bm{\theta}}_{i}\}$, which determine the node-specific models $\{\bm{\theta}_{v}\}$. This alternating structure allows CADO to efficiently minimize the joint objective without explicitly enumerating the (infinite) set of atoms.

Each update step uses a conditional gradient (Frank–Wolfe) approach. At iteration $t$, a linear minimization oracle is solved for both the model parameters and the embedding vectors to find descent directions. Specifically, the LMOs take the following forms:

[TABLE]

where $\bar{\epsilon}_{i,v}$ denotes the $v$-th entry of the embedding vector $\bar{\mathbf{e}}_{i}$. The exact solutions of these LMOs depend on the specific form of the loss functions $f_{v}$ and the domain $\Theta$. In the application considered in Section 5.1, we show that both LMOs admit efficient closed-form solutions; these derivations and the detailed steps of the CADO algorithm are provided in Appendix D.

After obtaining the LMOs, the variables are updated via a convex combination with the previous iterates, ensuring feasibility throughout the iterations. While the alternating conditional gradient approach leads to a non-convex problem, in practice CADO converges quickly to stable solutions, as supported by our experimental results. The procedure is summarized in Algorithm 1.

5 Experimental Studies

In this section, we first present the specific case under study, for which the closed-form solutions of (14) and (13) are derived, and the specialized version of CADO is summarized in Algorithm 2. We then evaluate the proposed framework and the CADO algorithm on the node classification task. First, we demonstrate the synergy between the graph structure, node features, and labels within our unified optimization framework. Second, we empirically validate the effectiveness of the CADO algorithm. Additional experimental results are provided in Appendix E to further support our theoretical and algorithmic findings.

5.1 Case Study

The frameworks in (6) can cover a vast variety of scenarios by different choices of $f_{v}$. To be more concrete, we focus on one example. We assume $\mathbf{z}_{v}$ vectors are statistically independent and in the same class, they follow an identical distribution. Moreover, in the training group, given the true class label $y_{v}$, the features and the noisy labels are independent, i.e. $p(\mathbf{x}_{v},\tilde{y}_{v}\mid y_{v}=i)=p(\mathbf{x}_{v}\mid\mathbf{R}_{i})p(\tilde{y}_{v}\mid\bm{\pi}_{i})$.

We take $m-$dimensional Gaussian, centered feature vectors $\mathbf{x}_{v}\sim\mathcal{N}(0,\bar{\mathbf{R}}_{i})$ with the covariance matrix $\bar{\mathbf{R}}_{i}$, in each class. We assume that the features of each class are close to a different linear subspace. Specifically, the eigenvalues of $\bar{\mathbf{R}}_{i}$ are divided into two groups with large and small values. The large ones, which are greater than a fixed positive value $\rho_{+}$ correspond to a signal subspace, while the small ones, being less than another fixed value $\rho_{-}$ correspond to a noise subspace. In the training set of class $i$, we assume that $\tilde{y}_{v}=j$ occurs with a probability $\bar{\pi}_{ji}$ and denote $\bar{\bm{\pi}}_{i}=(\bar{\pi}_{ji})_{j}$. Naturally, we require $\rho_{+}>\rho_{-}$ and $\bar{\pi}_{ii}$ to be significantly larger than $\bar{\pi}_{ji}$ for $j\neq i$.

To incorporate these models in (6), we choose $\Theta$ and $f_{v}$ as follows. We take two type of functions that aim to estimate $\bar{\mathbf{R}}_{i}$ and $\bar{\bm{\pi}}_{i}$:

[TABLE]

Then, the loss function of each node is taken to be:

[TABLE]

We search for $\mathbf{R}$ over the set $\mathbb{S}_{\rho_{-},\rho_{+}}$ of all symmetric matrices whose eigenvalues are in the interval $[\,\rho_{-},\ \rho_{+}]$. In other words,

[TABLE]

Similarly, we search for $\bm{\pi}$ over $\Delta=\left\{\bm{\pi}=(\pi_{k})\mid\pi_{k}\geq 0,\ \sum_{k}\pi_{k}=1\right\},$ i.e., the standard simplex. Accordingly, we may take $\Theta=\mathbb{S}_{\rho_{-},\rho_{+}}\times\Delta$.

5.2 Experimental Setup

Data generation model.

We construct a synthetic dataset with $K$ clusters of equal size ($n=Kn_{0}$), using a stochastic block model (SBM) for the graph and Gaussian mixtures for node features.

Graph structure: Edges are generated such that nodes within the same cluster connect with probability $p$, and between clusters with probability $q$.

Node features: Each cluster generates centered Gaussian features with a covariance matrix whose eigenvalues encode signal and noise power. Specifically, $m_{\omega}$ eigenvalues are set to $\omega^{2}$ (noise), and the remaining eigenvalues are set equally to $\sigma^{2}$. The values of $\omega$ and $m_{\omega}$ control the noise level and govern the separability between feature subspaces.

node labels: Each cluster contributes $n_{\mathrm{t}}$ training labels. A training node is correctly labeled with probability $\pi$, and mislabeled uniformly at random otherwise. This allows us to explore the effect of label noise and training set size.

Hyperparameters.

We summarize the key hyperparameter here; a comprehensive list is provided in Appendix E. Unless explicitly varied in the experiments, we use the following default settings: we consider $K=3$ clusters with $n_{0}=300$ nodes per cluster, and set the number of atoms to $r=K$. The intra-cluster and inter-cluster edge probabilities are $p=0.1$ and $q=0.05$, respectively. Node features are 6-dimensional, with $m_{\omega}=4$ noisy eigenvalues of magnitude $\omega=0.04$ in their covariance matrix. This ensures reasonable separability of feature subspaces. For labels, we set the training label ratio to $\nicefrac{{n_{\mathrm{t}}}}{{n_{0}}}=0.2$ and the correct label probability to $\pi=1$.

For simplicity in experiments, we reparameterize the original scaling setup by introducing effective weights for each term. While our theoretical framework uses a scaling of $1$ for the graph term, $\mu$ for the feature term, and $\mu\beta$ for the label term, in practice we fix the weights to $\beta_{g}=1.0$, $\beta_{f}=2.5$, and $\beta_{l}=13.0$, corresponding to the graph, feature, and label terms, respectively.

Evaluation settings.

To evaluate the contribution of each information source—graph structure, node features, and node labels—we conduct an ablation study by selectively enabling or disabling each term in our objective. Specifically, we evaluate configurations using only the graph term, only the feature term, all pairwise combinations, and the full model integrating all three. These ablations are implemented by setting the corresponding regularization weights to zero and solving the resulting problem using the CADO algorithm, except for the graph-only setting, where we apply spectral clustering instead of the ablated CADO variant.

5.3 Our Results

We present three sets of experiments in Figure 1, each analyzing the impact of one component while keeping others fixed. For each configuration, we report accuracy on test nodes.

Graph structure.

We vary intra-community probability $p$ while keeping features and labels fixed. Our framework outperforms graph-only and pairwise models, especially as the graph becomes less informative (i.e., $p\approx q$).

Feature quality.

We vary the noise level $\omega$ in the feature covariance, degrading the signal-to-noise ratio. The results show that our method remains robust and leverages the graph structure when features become unreliable.

Training label availability.

We vary the training label ratio. Our framework effectively exploits a small number of noisy labels, when utilizing the graph and feature information.

Detailed discussions of these experiments together with extra simulations to support our theoretical and algorithmic contributions are provided in Appendix E. Overall, our results confirm that the proposed framework effectively combines all available information sources. CADO consistently recovers the true structure when the theoretical conditions hold and remains robust under increasing difficulty ($p\downarrow$, or $\omega\uparrow$). The first plot shows how features and labels boost spectral clustering when the graph is weak, while the last two demonstrate how graph structure enhances classification when node-specific information is degraded.

6 Related Works

We are not aware of prior work that formulates node classification through the lens of spectral clustering and atomic decomposition, nor that proposes a convex relaxation via atomic norm and SON regularization. Our approach bridges this gap by combining tools from spectral graph theory and convex optimization to analyze and solve node classification tasks.

Spectral Graph Clustering.

Classical spectral clustering methods (Belkin & Niyogi, 2001; Ng et al., 2001) uncover graph community structures by embedding nodes into low-dimensional spaces using eigenvectors of the graph Laplacian. Hajek et al. (2016) has explored convex relaxations, such as semidefinite programming, to provide exact recovery guarantees under stochastic block models. Our work builds on this perspective by interpreting clustering as an atomic decomposition problem and extending it to incorporate node-specific information (features and labels).

Transductive Node Classification.

Transductive node classification aims to infer labels for unlabeled nodes using a graph and partially observed labels. Classical approaches (Zhu et al., 2003) propagate label information under smoothness assumptions. More recently, Graph Neural Networks (GNNs) such as GCN (Kipf & Welling, 2017), GraphSAGE (Hamilton et al., 2017), and GAT (Veličković et al., 2018) have achieved strong empirical results by aggregating information from node neighborhoods. While GNNs implicitly encode smoothness and low-rank priors, they lack convexity and are difficult to analyze theoretically. In contrast, our formulation offers a principled convex framework that enables theoretical recovery guarantees under homophily.

Concluding Remarks.

In this paper, we proposed a novel optimization framework for node classification from a spectral clustering perspective. We introduced an algorithmic solution for the non-convex formulation with a fixed number of atoms and established theoretical guarantees for perfect recovery in the convex, regularized version of the framework. Our theoretical results, supported by experimental studies, demonstrate that leveraging both the graph structure and node-specific information reveals a strong synergy that leads to improved classification accuracy.

Mathematical Notation.

To simplify notation and enhance readability, we use the following conventions throughout the paper. Collections such as $\{\lambda_{a}\}_{a\in\mathcal{A}}$, $\{L_{vv}\}_{v\in\mathcal{V}}$, $\{L_{uv}\}_{u,v\in\mathcal{V}}$, and $\{\bm{\theta}_{v}\}_{v\in\mathcal{V}}$ are often abbreviated to $\{\lambda_{a}\}$, $\{L_{vv}\}$, $\{L_{uv}\}$, and $\{\bm{\theta}_{v}\}$ when the indexing set is clear from context. When two nodes $u$ and $v$ belong to the same cluster, we write $u,v\in\mathcal{C}_{i}$; otherwise, we denote $u\in\mathcal{C}_{i}$ and $v\in\mathcal{C}_{j}$ to indicate distinct clusters.

We denote vectors by bold lowercase letters (e.g., $\mathbf{x}$) and matrices by bold uppercase letters (e.g., $\mathbf{W}$). The entry in the $v^{\text{th}}$ row and $u^{\text{th}}$ column of a matrix $\mathbf{W}$ is denoted by $W_{vu}$, and its transpose by $\mathbf{W}^{\top}$. The Euclidean inner product between vectors is written as $\langle\cdot,\cdot\rangle$, and the matrix inner product is defined by $\langle\mathbf{A},\mathbf{C}\rangle:=\mathrm{tr}(\mathbf{A}\mathbf{C}^{\top})$, where $\mathrm{tr}(\cdot)$ denotes the trace operator.

The indicator function $I_{\Theta}(\bm{\theta})$ refers to the indicator of the feasible set $\Theta$, and $\partial I_{\Theta}(\bm{\theta})$ denotes the normal cone of $\Theta$ at point $\bm{\theta}$. The cardinality of a set $T$ is denoted by $|T|$. We use $\mathbf{1}_{n}$ and $\mathbf{0}_{n}$ to denote the $n$-dimensional all-one and all-zero vectors, respectively. An $m\times n$ all-zero matrix is denoted by $\mathbf{O}_{m\times n}$, with subscripts omitted when dimensions are clear from context.

Appendix A Atomic Norm Review

In many signal processing and machine learning applications, the objective is to reconstruct a signal that admits a simple or structured representation. The concept of atomic norms provides a unifying framework for this purpose by promoting simplicity in the signal representation. Specifically, an atomic norm can be used as a convex regularizer that induces a desired structure, such as sparsity or low-rankness, by relying on a predefined set of fundamental building blocks called atoms.

A.1 Atoms and Atomic Norms

The basic elements used to represent a signal are referred to as atoms. For a signal $\mathbf{x}$, the atomic set $\mathcal{A}$ is a collection of these basic elements, allowing the signal to be expressed as a nonnegative linear combination of a small number of atoms from $\mathcal{A}$. The atomic norm, denoted by $\|\mathbf{x}\|_{\mathcal{A}}$, quantifies the complexity of the signal in terms of how economically it can be represented using these atoms.

Formally, the atomic norm is defined as:

[TABLE]

This norm is often employed as a regularizer in optimization problems to encourage solutions that are structurally simple, such as sparse vectors or low-rank matrices. Different choices of the atomic set $\mathcal{A}$ yield different atomic norms. Two important examples are as follows:

• •

Sparsity ($\ell_{1}$-norm): For sparse signals, the atomic set consists of the signed canonical basis vectors:

[TABLE]

where $\mathbf{e}_{i}$ denotes the $i$-th canonical basis vector in $\mathbb{R}^{n}$. The induced atomic norm is:

[TABLE]

which is the familiar $\ell_{1}$-norm, widely used to promote sparsity.

• •

Low-rank matrices (nuclear norm): For low-rank matrix recovery, such as in matrix completion or graph clustering, the atomic set consists of rank-one matrices with unit norm:

[TABLE]

The induced atomic norm is:

[TABLE]

which is the nuclear norm, equal to the sum of the singular values of $\mathbf{X}$. This norm is the convex surrogate of the rank function and promotes low-rank structure.

For symmetric matrices, the above definitions can be simplified. If $\mathbf{X}=\mathbf{X}^{\top}$, we can define the symmetric atomic sets

[TABLE]

Then, with the eigenvalue decomposition $\mathbf{X}=\sum_{i}\lambda_{i}\,\mathbf{u}_{i}\mathbf{u}_{i}^{\top}$ (like graph Laplacian matrix):

[TABLE]

and, for general symmetric $\mathbf{X}$,

[TABLE]

In particular, for symmetric matrices the nuclear norm equals the sum of absolute eigenvalues, and for symmetric PSD matrices it reduces to the trace.

Appendix B Proof of Theorem 4.3

For convenience, we first provide a brief overview of the proof techniques employed.

B.1 Proof Overview

We cast the problem in (5) as a constrained atomic norm minimization, which we solve using a tailored ADMM-based algorithm. This method decouples the atomic norm regularization from additional structural constraints via a variable splitting strategy. To establish the complexity result, we demonstrate that our problem meets the conditions required for standard ADMM convergence, yielding an overall rate of $O(1/T)$.

B.2 Constrained Atomic Norm Optimization (CANO)

We consider the following general constrained atomic norm optimization problem, referred to as CANO:

[TABLE]

where $\|\cdot\|_{\mathcal{A}}$ is the atomic norm, and $\mathcal{X}\subseteq\mathbb{R}^{d}$ is a closed and convex set encoding additional constraints (e.g., feasibility constraints in (5)). This formulation introduces computational challenges due to the coupling of a non-polyhedral norm and convex constraint sets.

B.3 CANO-ADMM Algorithm

To solve (20), we apply ADMM by introducing an auxiliary variable $\mathbf{z}$ and rewriting the problem as:

[TABLE]

The augmented Lagrangian is defined as:

[TABLE]

where $\bm{\lambda}$ is the dual variable, and $\beta>0$ is the penalty parameter. The ADMM updates proceed as follows:

• •

Update $\mathbf{x}$:

[TABLE]

This step is solved using the CoGEnT algorithm Rao et al. (2015), which is designed for atomic norm-constrained problems.

• •

Update $\mathbf{z}$:

[TABLE]

This corresponds to projecting onto the set $\mathcal{X}$, which is assumed to be tractable via a gradient-projection update:

[TABLE]

where $\alpha$ is a step size and $\mathsf{P}_{\mathcal{X}}$ denotes Euclidean projection onto $\mathcal{X}$.

• •

Update $\bm{\lambda}$:

[TABLE]

B.4 Convergence Guarantees of CANO-ADMM

The convergence of ADMM with a rate of $O(1/T)$ in objective residuals and constraint violations is guaranteed under the following conditions:

$f(\mathbf{x})$ is convex and has a Lipschitz-continuous gradient; 2. 2.

The constraint sets $\|\mathbf{x}\|_{\mathcal{A}}\leq\tau$ and $\mathcal{X}$ are both closed and convex; 3. 3.

The subproblems are solvable to sufficient accuracy at each iteration.

Our setup satisfies all these assumptions:

• •

The function $f(\mathbf{x})$ is convex and differentiable (e.g., it is the sum of a linear term and convex node-wise losses);

• •

The atomic norm ball $\|\mathbf{x}\|_{\mathcal{A}}\leq\tau$ is convex by definition, and $\mathcal{X}$ is assumed to be convex and compact;

• •

The projection and LMO steps in CoGEnT are computationally tractable and converge efficiently.

Hence, CANO$-$ADMM converges at a rate $O(1/T)$ in the number of iterations. As each iteration requires a linear minimization oracle and projection operation, the overall runtime is polynomial in $1/\epsilon$, completing the proof of Theorem 4.3.

Appendix C Proof of Theorem 3.6

We begin with an overview of the proof methodology, followed by the problem setup and optimization formulation. We then outline the key elements, focusing on optimality conditions for exact recovery. The main proof is divided into three parts, each covering a different information regime: (i) node-specific, (ii) graph-based, and (iii) combined. For each, we construct a dual certificate via a “guess-and-golfing” approach and establish recovery conditions. We conclude by stating the theorem.

C.1 Proof Overview

Our proof follows a standard dual certificate approach. We aim to verify that the ideal solution—comprising one atom per cluster—satisfies the optimality (KKT) conditions of the convex problem (6). The core idea is to construct a dual certificate matrix $\mathbf{Z}$ (and associated subgradients $\{\mathbf{g}_{vu}\}$) such that the ideal solution is locally optimal.

To achieve this, we employ a guess-and-golfing strategy. Starting from a natural or trivial guess for $\mathbf{Z}$ (e.g., diagonal or block-structured), we iteratively refine it to meet the necessary conditions. Each refinement step aims to adjust the certificate to satisfy one or more of the following:

• •

sign constraints ensuring dual feasibility;

• •

inner product constraints capturing primal-dual consistency;

• •

and positive semidefiniteness of the gap matrix enforcing inequality optimality.

C.2 Problem Setup

Optimization Problem Restatement.

The optimization problem we analyze, with the goal of establishing conditions for perfect recovery, is our convex framework augmented with a sum-of-norms (SON) regularization term. For clarity and completeness, we restate it below, as in equation (6).

[TABLE]

In this formulation, $\mathbf{U}:=\left(\mathbf{L},\{\bm{\theta}_{v}\}_{v\in\mathcal{V}}\right)$. The feasible set is $\mathcal{U}:=\mathcal{B}\times\Theta^{|\mathcal{V}|}$, where $\mathcal{B}$ denotes the set of symmetric matrices with entries in $[0,1]$ and unit diagonal, and $\Theta$ is the feasible domain for node-specific models. The objective $\phi(\mathbf{U})$ is defined as: $\phi(\mathbf{U}):=-\braket{\bar{\mathbf{A}},\mathbf{L}}+\mu\sum_{v\in\mathcal{V}}f_{v}(\mathbf{z}_{v};\bm{\theta}_{v}),$ and $R(\mathbf{U})$ is the SON regularization term. Each atom $\in\mathcal{A}$ is of the form $\left(\mathbf{e}\mathbf{e}^{\top},\{\epsilon_{v}^{2}\bm{\theta}\}_{v\in\mathcal{V}}\right)$, and comes from the atomic dictionary $\mathcal{A}$ defined as:

[TABLE]

Using the definitions above, the optimization problem (26) can be explicitly given as:

[TABLE]

Equivalent Optimization Problem.

We simplify the problem in (28) by enforcing $L_{vv}=1$ for all $v$, which implies that the sum of coefficients $\lambda_{a}$ is fixed, and the conditions $L_{uv}\leq 1$ and $\bm{\theta}_{v}\in\Theta$ are automatically satisfied (see Lemma C.2). We also eliminate the auxiliary variables $\bm{\theta}_{v}$ by expanding their definition in the objective, leading to the following simplified problem:

[TABLE]

C.3 Proof Elements

Ideal Solution.

We define the ideal (perfect recovery) solution as one consisting of $K$ atoms—one per cluster—denoted by $\{(\hat{\lambda}_{i},\hat{\mathbf{e}}_{i},\hat{\bm{\theta}}_{i})\}_{i=1}^{K}$. These ideal atoms take the form:

[TABLE]

This yields an ideally partitioned matrix $\mathbf{L}$ with:

[TABLE]

Characteristic Optimization.

The optimal models $\hat{\bm{\theta}}_{i}$, which we refer to as biased centroids, are obtained by solving the following characteristic optimization problem, which jointly minimizes the empirical loss and the sum-of-norms penalty across clusters:

[TABLE]

where $F_{i}(\bm{\theta}):=\frac{1}{n_{i}}\sum_{v\in\mathcal{C}_{i}}f_{v}(\mathbf{z}_{v};\bm{\theta})$ is the aggregate loss over cluster $\mathcal{C}_{i}$. Note that the characteristic optimization in (32) corresponds to solving the full problem (29), given the ideal atoms $\{\hat{\lambda}_{i},\hat{\mathbf{e}}_{i}\}_{i=1}^{K}$ specified in (30). If we define $\gamma:=\nicefrac{{n\mu_{1}}}{{\mu}}$, the biased centroids satisfy the following condition, which derives from the first-order optimality condition of (32):

[TABLE]

Optimality Conditions.

To establish that the ideal solution is optimal for the simplified problem (29), it suffices to verify the first-order optimality (KKT) conditions at that point. Specifically, we aim to construct a dual certificate: a symmetric matrix $\mathbf{Z}\in\mathbb{R}^{n\times n}$ and a collection of subgradients $\{\mathbf{g}_{vu}\}$ such that the following optimality conditions are satisfied:

[TABLE]

[TABLE]

Condition (34) ensures complementary slackness with respect to the matrix variable $\mathbf{L}$. Also, importantly, the diagonal elements of $\mathbf{Z}$ are unconstrained and need not satisfy the sign conditions in (34), as they are not involved in the objective due to the fixed diagonal constraint $L_{vv}=1$. Inequality (35) must hold for all atoms $(\mathbf{e},\bm{\theta})\in\mathcal{A}$, while equality (36) holds for the $K$ ideal atoms $\{(\hat{\mathbf{e}}_{i},\hat{\bm{\theta}}_{i})\}_{i=1}^{K}$.

The terms $\mathbf{g}_{vu}$ in (35)–(36) correspond to subgradients of the SON regularization term, and are defined as:

[TABLE]

By symmetry of the SON term, we additionally require that $\mathbf{g}_{vu}=-\mathbf{g}_{uv}$. Hence, the KKT conditions reduce to constructing a dual matrix $\mathbf{Z}$ and a consistent set of subgradients $\{\mathbf{g}_{vu}\}$ that jointly satisfy equations (34)–(37).

Simplified Optimality Conditions.

We simplify the optimality conditions (35)–(36) by introducing a diagonal matrix $\mathbf{D}(\bm{\theta}):=\mathrm{diag}(d_{v}(\bm{\theta}))$, where each diagonal entry is defined as

[TABLE]

In the final line above, we used the optimality condition of the characteristic optimization (33) and substituted $\gamma=\nicefrac{{n\mu_{1}}}{{\mu}}$. Here, $y_{v}$ is the index of the cluster containing node $v$. Using this definition, we can equivalently rewrite the original optimality conditions as:

[TABLE]

Note that $\mathbf{g}_{vu}$ is undefined when $u,v\in\mathcal{C}_{i}$ (i.e., both nodes belong to the same cluster). In this case, $\mathbf{g}_{vu}$ acts as a free design parameter and must satisfy the subgradient constraint $\|\mathbf{g}_{vu}\|\leq 1$. We jointly construct the collection $\{\mathbf{g}_{vu}\}_{u,v\in\mathcal{C}_{i}}$ and the dual certificate $\mathbf{Z}$ such that the full set of optimality conditions—(34), (39)-(40), along with the subgradient condition (37)—are all satisfied.

C.4 Perfect Recovery with Node-Specific Information Only

We begin by constructing a simple candidate for the dual certificate $\mathbf{Z}$ and the subgradients $\{\mathbf{g}_{vu}\}$ that satisfy the optimality conditions in the absence of graph information. Specifically, we consider a diagonal matrix $\mathbf{Z}$ defined as:

[TABLE]

where $d_{v}(\bm{\theta})$ is given in equation (38). Due to the structure of the ideal atoms $\{\hat{\mathbf{e}}_{i}\}$, this choice guarantees that condition (40) is satisfied.

To enforce condition (39), we define:

[TABLE]

for which we have the following result that ensures(39).

[TABLE]

By Assumption 3.5, we have $\|\nabla f_{v}(\hat{\bm{\theta}}_{i})-\nabla f_{u}(\hat{\bm{\theta}}_{i})\|\leq\rho$, which ensures $\|\mathbf{g}_{vu}\|\leq 1$ if and only if:

[TABLE]

This defines the condition under which perfect recovery is guaranteed using node-specific information alone. Note that this construction does not rely on the graph, and condition (34) holds trivially if $\bar{\mathbf{A}}$ corresponds to an ideal graph.

C.5 Perfect Recovery with Graph Information Only

In the absence of node-specific information (i.e., when $\mu=0$), the optimization reduces to pure graph clustering. In this case, the objective becomes minimizing $-\langle\bar{\mathbf{A}},\mathbf{L}\rangle$ over convex combinations of rank-one outer products, and the sum-of-norms regularizer becomes inactive. Consequently, the dual certificate $\mathbf{Z}$ must satisfy the following simplified optimality conditions:

[TABLE]

Our goal is to construct a certificate $\mathbf{Z}$ that satisfies these conditions, certifying that the ideal solution is indeed optimal. We proceed via a step-by-step design using a golfing strategy, where we iteratively refine a candidate $\mathbf{Z}$ to satisfy each condition.

Step 1: Structured Guess.

We begin with a matrix $\mathbf{S}\in\mathbb{R}^{n\times n}$ defined blockwise according to:

[TABLE]

where $a_{ii},a_{ij}>0$ are scalar parameters. This form penalizes incorrect intra-cluster disconnections and inter-cluster connections. For brevity, we will denote this matrix as $\mathbf{S}=(0,-a_{ii},0,a_{ij})$ in the remainder of the proof. This step helps to satisfy the condition (43).

Step 2: Projection for Orthogonality.

To ensure that $\mathbf{Z}$ is orthogonal to the ideal atoms (i.e., satisfies (45)), we project $\mathbf{S}$ onto the orthogonal complement of the span of ideal atoms:

[TABLE]

and $\mathbf{E}=[\hat{\mathbf{e}}_{1},\dots,\hat{\mathbf{e}}_{K}]$ is the matrix of ideal cluster indicators. This construction ensures $\hat{\mathbf{e}}_{i}^{\top}\mathbf{Z}\hat{\mathbf{e}}_{i}=0$ for all $i$, as required by (45).

Step 3: Sign Condition.

Next, we enforce the sign constraints (43) by analyzing the entrywise form of $\mathbf{Z}-\bar{\mathbf{A}}$. Using Assumptions 3.1–3.2, we compute worst-case deviations after projection and derive bounds:

[TABLE]

Step 4: Negative Definiteness.

Finally, to ensure $\mathbf{Z}\preceq 0$ as required by (44), we must choose $\lambda$ large enough so that $\mathbf{P}^{\perp}(\mathbf{S}-\lambda\mathbf{I})\mathbf{P}^{\perp}$ has all non-positive eigenvalues. A sufficient condition is:

[TABLE]

where $a:=\max\{a_{ii},a_{ij}\}$ and $\tilde{\mathbf{A}}$ is the centered adjacency matrix introduced in the Section 3.1.1. Random matrix theory gives $\lambda_{\max}(\tilde{\mathbf{A}})=\Omega(\sqrt{n})$, so $\lambda=\Omega(\sqrt{n})$ suffices. To ensure that the sign conditions remain valid, we must also enforce $\lambda/n_{i}\to 0$, which yields the requirement:

[TABLE]

This confirms that perfect recovery is achievable using graph information alone, provided that inter- and intra-cluster connectivities are well-separated (Assumptions 3.1–3.2), and the number of nodes grows at most quadratically with respect to the cluster sizes (i.e., the number of clusters does not grow too quickly).

C.6 Perfect Recovery with Combined Graph and Node-specific Information

The goal of this section is to demonstrate that by combining graph structure and node-specific information, we obtain looser requirements on the parameters $\lambda$, $\mu$, and $\rho$ for perfect recovery. We construct a dual certificate $\mathbf{Z}$ that incorporates both sources of information and verify when it satisfies the optimality conditions (34), (39), and (40).

We define the dual certificate as:

[TABLE]

and the subgradients $\mathbf{g}_{vu}$ as

[TABLE]

Verifying Optimality Conditions.

We begin with (40), which holds due to $\bar{\mathbf{D}}$’s diagonal structure and the projection $\mathbf{P}^{\perp}$.

Next, we focus on verifying (39). Substituting the dual certificate $\mathbf{Z}$, we get:

[TABLE]

Bounding Node-specific Term.

Under Assumption 3.4, there exists at most one index $i$ for which

[TABLE]

Without loss of generality, assume this is index $i=0$. We define the diagonal matrix $\mathbf{F}\in\mathbb{R}^{n\times n}$ by:

[TABLE]

where only the first entry (corresponding to $i=0$) is [math] and others are $R$. Thus,

[TABLE]

Bounding Gradient Deviation Term.

Assuming the feasible set $\Theta$ is bounded as:

[TABLE]

and using the Cauchy-Schwarz inequality, we get:

[TABLE]

Combining Bounds.

Putting the above together into (48), we obtain:

[TABLE]

We want the matrix

[TABLE]

to be positive semidefinite. This matrix has the same block structure as $\tilde{\mathbf{A}}$, except that all diagonal entries are shifted by $\mu R+\lambda-\zeta$, except in block $i=0$, which is shifted by $\lambda-\zeta$.

By invoking Lemma C.3, we can determine sufficient conditions on $\lambda,\mu,\rho$ that ensure positive semidefiniteness and hence perfect recovery.

We now present the improved conditions for perfect recovery by leveraging both graph and node-specific information. The result is split into two regimes depending on the relative size of cluster size $n_{i}$ and the total number of nodes $n$.

Case 1: Many small clusters ($n_{i}\leq\sqrt{n}$).

In this case, we assume the number of clusters is large. We choose:

[TABLE]

to satisfy the assumptions of Lemma C.3. Then the condition: $\rho\mu-\mu_{1}n_{i}\leq n_{i}$ implies:

[TABLE]

This setting yields a substantial gain over the graph-only regime, particularly in high-cluster settings.

Case 2: Few large clusters ($n_{i}\geq\sqrt{n}$).

Here, the number of clusters is small and each cluster contains many nodes. We choose:

[TABLE]

Using similar logic, we obtain the recovery condition:

[TABLE]

Both cases demonstrate clear improvement; case 1 provides a strong gain due to tighter scaling with $\nicefrac{{n_{i}^{2}}}{{n}}$, and case 2 achieves a meaningful gain for larger clusters, improving over the baseline graph-only condition.

C.7 Final Theorem

Theorem C.1**.**

Assume a sequence of problems where $n,K$, and all $n_{i}$s grow to infinity such that $n_{i}\sim n_{j}$ and all $\rho_{ij}$ converge to fixed nonzero values. Fix $\gamma=\nicefrac{{n\mu_{1}}}{{\mu}}$ and suppose that Assumptions 3.1–3.4 hold with fixed $\delta$ and $R$, and that the feasible set $\Theta$ is bounded. Then, the followings hold true:

Depending on $\gamma,\delta$, there exists a constant $c>1$ such that for

[TABLE]

the ideal solution with $\mathbf{e}_{i}$, $\lambda_{i}=n_{i}$, and $\bm{\theta}=\hat{\bm{\theta}}_{i}$ is optimal in (6), provided Assumption 3.5 holds and

[TABLE] 2. 2.

Without the graph term $-\braket{\bar{\mathbf{A}},\mathbf{L}}$, the ideal solution is optimal if $c=1$, i.e., Assumption 3.5 holds with $\rho\leq\gamma\cdot\nicefrac{{n_{i}}}{{n}}.$ 3. 3.

Without the node-specific terms, graph clustering recovers clusters if $n=O(n_{i}^{2})$.

C.8 Lemmas

Lemma C.2**.**

Let $\bm{\theta}_{v}:=\sum_{a\in\mathcal{A}}\lambda_{a}\epsilon_{a,v}^{2}\bm{\theta}_{a}$ and $\mathbf{L}:=\sum_{a\in\mathcal{A}}\lambda_{a}\mathbf{e}\mathbf{e}^{\top}$, where each $\bm{\theta}_{a}\in\Theta$ and $\|\mathbf{e}\|=1$. Assume that for all $v\in\mathcal{V}$, $\sum_{a\in\mathcal{A}}\lambda_{a}\epsilon_{a,v}^{2}=1.$ Then, it follows that $\bm{\theta}_{v}\in\Theta$, $L_{uv}\leq 1$ for all $u,v\in\mathcal{V}$, and $\sum_{a\in\mathcal{A}}\lambda_{a}=n$.

Proof.

If $\sum_{a\in\mathcal{A}}\lambda_{a}\epsilon_{a,v}^{2}=1$ holds, we have:

• •

Since each $\bm{\theta}_{a}\in\Theta$ and $\sum_{a\in\mathcal{A}}\lambda_{a}\epsilon_{a,v}^{2}=1$, it follows that $\bm{\theta}_{v}$ is a convex combination of elements in $\Theta$, and hence $\bm{\theta}_{v}\in\Theta$ due to the convexity of $\Theta$.

• •

For $\mathbf{L}=\sum_{a\in\mathcal{A}}\lambda_{a}\mathbf{e}\mathbf{e}^{\top}$, the Cauchy–Schwarz inequality implies

[TABLE]

• •

Summing the constraint over all $v$, considering $\|\mathbf{e}\|^{2}=1$, gives:

[TABLE]

∎

Lemma C.3**.**

Consider the symmetric block matrix

[TABLE]

where $\mathbf{A},\mathbf{C}$ are symmetric matrices satisfying $\mathbf{A}\succeq\alpha\mathbf{I}$ and $\mathbf{C}\succeq\beta\mathbf{I}$ for some $\alpha,\beta\geq 0$. Then $\mathbf{D}$ is positive semidefinite if $\sigma^{2}(\mathbf{B})\leq\alpha\beta,$ where $\sigma(\mathbf{B})$ denotes the largest singular value of $\mathbf{B}$.

Proof.

Take any vector $\mathbf{x}=[\mathbf{x}_{1}^{\top},\mathbf{x}_{2}^{\top}]^{\top}$. Then, the quadratic form associated with $\mathbf{D}$ can be written as:

[TABLE]

Using the given conditions $\mathbf{A}\succeq\alpha\mathbf{I}$ and $\mathbf{C}\succeq\beta\mathbf{I}$, we have: $\mathbf{x}_{1}^{\top}\mathbf{A}\mathbf{x}_{1}\geq\alpha\|\mathbf{x}_{1}\|^{2},\mathbf{x}_{2}^{\top}\mathbf{C}\mathbf{x}_{2}\geq\beta\|\mathbf{x}_{2}\|^{2}.$ Applying these lower bounds and the definition of singular value ($\|\mathbf{B}\mathbf{x}_{2}\|\leq\sigma(\mathbf{B})\|\mathbf{x}_{2}\|$), it follows by the Cauchy–Schwarz inequality that:

[TABLE]

Now, complete the square explicitly to factorize the expression clearly:

[TABLE]

Since $\sigma^{2}(\mathbf{B})\leq\alpha\beta$ (implying $\sigma(\mathbf{B})\leq\sqrt{\alpha\beta}$), both terms in this expression are nonnegative. Thus:

[TABLE]

Hence, $\mathbf{D}$ is positive semi-definite. ∎

Appendix D Details of the CADO Algorithm

This section provides the full implementation details of the CADO algorithm. We describe how the embedding vectors and model parameters are updated using conditional gradient steps, relying on problem-specific linear minimization oracles (LMOs). Each component is addressed in a separate subsection, followed by the complete algorithm pseudo-code.

D.1 Embedding Update via Conditional Gradient

In this step, we update the embedding vectors $\{\bar{\mathbf{e}}_{i}\}_{i=1}^{r}$, which define the low-rank matrix $\mathbf{L}=\sum_{i}\bar{\mathbf{e}}_{i}\bar{\mathbf{e}}_{i}^{\top}$, while keeping the class-wise models $\{\bar{\bm{\theta}}_{i}\}$ fixed. The update is performed by solving the following LMO:

[TABLE]

where $\phi$ is the objective function from (4):

[TABLE]

Constraint simplification.

We now simplify the constraints using the structure of $\mathbf{L}$ and $\bm{\theta}_{v}$. First, the constraint $\mathbf{L}\in\mathcal{B}$ requires: $L_{vv}=1$ for all $v$, $0\leq L_{uv}\leq 1$ for all $u,v$, and symmetry.

Given $\mathbf{L}=\sum_{i}\bar{\mathbf{e}}_{i}\bar{\mathbf{e}}_{i}^{\top}$, the diagonal condition $L_{vv}=1$ becomes:

[TABLE]

This implies that the vector $\bar{\bm{\epsilon}}_{v}=(\bar{\epsilon}_{i,v})_{i}$ lies on the unit sphere, and that:

[TABLE]

is a convex combination of class-wise models, hence satisfying $\bm{\theta}_{v}\in\Theta$ automatically. Also, under this condition:

[TABLE]

and the upper bound constraint $L_{uv}\leq 1$ is also satisfied. The only remaining constraint is the nonnegativity of $L_{uv}$, which is equivalent to:

[TABLE]

Reparameterization.

To simplify the optimization, we define: $W_{v,i}\coloneqq\bar{\epsilon}_{i,v}^{2}.$ Let $\mathbf{W}\in\mathbb{R}^{n\times r}$ collect these squared embedding values. Under this reparameterization, the low-rank matrix becomes $\mathbf{L}=\sum_{i=1}^{r}\sqrt{W_{:,i}}\sqrt{W_{:,i}}^{\top},$ and the node-specific models satisfy $\bm{\theta}_{v}=\sum_{i=1}^{r}W_{v,i}\bar{\bm{\theta}}_{i}.$ The constraint $L_{vv}=\sum_{i}W_{v,i}=1$ together with $W_{v,i}\geq 0$ implies that each row of $\mathbf{W}$ lies in the probability simplex.

Thus, the LMO reduces to:

[TABLE]

We omit the non-negativity constraint $\mathbf{L}\geq 0$, as it is empirically satisfied due to the structure of $\bar{\mathbf{A}}$. If needed, it can be enforced via ADMM.

Gradient Computation.

Let $\mathbf{R}_{v}=\sum_{i}W_{v,i}\bar{\mathbf{R}}_{i}$, $\bm{\pi}_{v}=\sum_{i}W_{v,i}\bar{\bm{\pi}}_{i}$. The gradient of the objective in (53) with respect to $W_{v,i}$ consists of:

• •

Structural term:

[TABLE]

• •

Feature loss:

[TABLE]

• •

Label loss (training nodes only):

[TABLE]

Solution.

The optimization in (53) is linear in each row $\mathbf{W}_{v,:}$ and constrained to the simplex. Therefore, the optimal solution is a one-hot vector with 1 at the coordinate corresponding to the minimum gradient value:

[TABLE]

where

[TABLE]

Once $\mathbf{W}$ is computed, we recover embeddings via the following equations. Taking positive roots preserves the symmetry and ensures non-negativity of $\mathbf{L}$.

[TABLE]

D.2 Model Update via Conditional Gradient

Given the updated embedding vectors $\{\bar{\mathbf{e}}_{i}\}_{i=1}^{r}$, we update the class-wise models $\{\bar{\bm{\theta}}_{i}=(\bar{\mathbf{R}}_{i},\bar{\bm{\pi}}_{i})\}_{i=1}^{r}$ by solving one LMO over each atom. Specifically, each model parameter is updated via:

[TABLE]

where $\phi$ is the global objective from (4), and the embedding matrix $\mathbf{W}$ (with $W_{v,i}=\bar{\epsilon}_{i,v}^{2}$) is fixed.

Since each $\bar{\bm{\theta}}_{i}$ consists of a feature model $\bar{\mathbf{R}}_{i}\in\mathbb{S}_{\rho_{-},\rho_{+}}$ and a label distribution $\bar{\bm{\pi}}_{i}\in\Delta$, the LMO separates into two independent problems:

[TABLE]

Gradient computation.

Let $\mathbf{R}_{v}=\sum_{i=1}^{r}W_{v,i}\bar{\mathbf{R}}_{i}$ and $\bm{\pi}_{v}=\sum_{i=1}^{r}W_{v,i}\bar{\bm{\pi}}_{i}$. Then, the gradients are given by:

• •

Feature loss (all nodes):

[TABLE]

where

[TABLE]

• •

Label loss (training nodes only):

[TABLE]

where

[TABLE]

Closed-form solutions.

Both optimization problems admit closed-form solutions:

• •

Update for $\bar{\mathbf{R}}_{i}$: The optimal solution to the linear minimization problem

[TABLE]

is given by

[TABLE]

where $\nabla_{\bar{\mathbf{R}}_{i}}\phi=\mathbf{U}\,\mathrm{diag}(\lambda_{1},\dots,\lambda_{m})\,\mathbf{U}^{\top}$ is the eigen-decomposition of the gradient, and

[TABLE]

• •

Update for $\bar{\bm{\pi}}_{i}$: Since the objective is linear over the simplex, the solution is the vertex corresponding to the smallest coordinate:

[TABLE]

These updates define the model step in each iteration of the CADO algorithm.

D.3 The Specialized CADO Algorithm

We now summarize the complete CADO algorithm specialized for the node classification setting studied in this paper. This algorithm solves the constrained atomic decomposition problem in (4) using a conditional gradient approach, alternating between updating the embedding vectors $\{\bar{\mathbf{e}}_{i}\}$ and the class-wise models $\{\bar{\bm{\theta}}_{i}=(\bar{\mathbf{R}}_{i},\bar{\bm{\pi}}_{i})\}$.

Embedding step.

In each iteration, the embedding update seeks a direction that reduces the global objective $\phi$ while maintaining feasibility. To make this step efficient, we reparameterize the squared embedding entries as $W_{v,i}=\bar{\epsilon}_{i,v}^{2}$, which allows us to enforce both the diagonal constraint on $\mathbf{L}$ and the convexity condition on $\bm{\theta}_{v}$. The resulting LMO admits a closed-form solution: each row of $\mathbf{W}$ is set to a one-hot vector in the direction of steepest descent.

Model step.

Given the updated embeddings, the model parameters $\bar{\bm{\theta}}_{i}=(\bar{\mathbf{R}}_{i},\bar{\bm{\pi}}_{i})$ are updated by solving LMOs over the model space $\Theta=\mathbb{S}_{\rho_{-},\rho_{+}}\times\Delta$. The gradients of the global objective $\phi$ with respect to both components are derived in closed form based on the structure of the loss functions. These LMOs also admit simple solutions: the covariance matrix $\bar{\mathbf{R}}_{i}$ is updated by projecting the negative gradient onto the spectral box, while the label distribution $\bar{\bm{\pi}}_{i}$ is updated by selecting the coordinate with the smallest gradient value.

Alternating optimization.

The algorithm alternates between these two steps, using a step size $\gamma_{t}=\frac{2}{t+2}$ at iteration $t$ to compute convex combinations of the previous and newly computed atoms. This ensures feasibility at all iterations and convergence under standard assumptions. The resulting procedure is efficient, scalable, and compatible with a wide range of feature and label models.

Final algorithm.

The complete specialized version of the CADO algorithm is presented in Algorithm 2 below.

Appendix E Extended experiments

E.1 Model Parameter Sensitivity

We study the sensitivity of the model’s performance with respect to the weighting parameters $\beta_{g}$, $\beta_{f}$, and $\beta_{l}$, which control the influence of the graph term, feature term, and label term in our unified objective. The results are presented in Figure. 2.

In the left panel, we vary $\beta_{g}$, the weight of the graph term. When $\beta_{g}=0$, the model essentially ignores the graph structure, which causes a significant drop in performance in the Graph + Label configuration. However, the full model (Graph + Feature + Label) maintains high accuracy even at $\beta_{g}=0$, indicating that feature and label information alone can provide a strong signal. As $\beta_{g}$ increases, accuracy improves across all settings, but plateaus after $\beta_{g}\approx 5$, suggesting diminishing returns from overly amplifying the graph signal.

In the center panel, we vary $\beta_{f}$, the weight of the feature term. When $\beta_{f}=0$, the performance of the Graph + Feature configuration drops sharply, as expected. Interestingly, the full model remains relatively robust, highlighting the complementary strength of the graph and label terms. Very large values of $\beta_{f}$ lead to performance degradation in some settings, possibly due to overfitting to noisy feature dimensions.

In the right panel, we sweep $\beta_{l}$, the weight of the label term. At $\beta_{l}=0$, the Graph + Label configuration performs poorly due to the absence of label supervision. Again, the full model is quite robust and achieves high accuracy even with limited label contribution. Increasing $\beta_{l}$ improves performance, but too high values lead to minor declines, likely because the model overemphasizes noisy or limited training labels.

Overall, these results confirm that our model is robust across a wide range of parameter choices and demonstrates strong synergy when all sources of information—graph, features, and labels—are integrated. The default values used in the main experiments strike a good balance across components.

E.2 Data Parameter Sensitivity

Figure˜3 investigates how the performance of our method changes as we vary key data generation parameters, specifically the number of nodes $n$, feature dimension $m$, and the number of clusters $c$.

In the left panel, increasing the number of nodes significantly improves test accuracy across all configurations, as larger graphs provide more structure and statistical signal. The full model (Graph + Feature + Label) consistently achieves the highest accuracy and converges quickly, even with a moderate number of nodes. This highlights the data efficiency of our joint framework, particularly when leveraging all available modalities.

In the center panel, we vary the total feature dimension $m$, which includes both signal and noise components. As $m$ increases, performance generally decreases—especially for methods that rely on features (e.g., Graph + Feature)—due to the increasing influence of noisy or uninformative dimensions. However, the full model (cyan) remains relatively stable and outperforms all other combinations, suggesting that combining features with graph structure and labels helps mitigate the curse of dimensionality.

In the right panel, we increase the number of clusters $c$, making the classification problem more challenging due to finer partitioning and weaker homophily. The performance of all configurations drops, but the full model retains significantly higher accuracy. Notably, using graph-only or label-only configurations fails beyond $c\approx 10$, while feature-based and joint models scale more gracefully. This result supports our theoretical findings: the integration of feature and label information compensates for reduced graph separability as the number of communities increases.

Overall, this experiment confirms the importance of combining multiple modalities. It also demonstrates the robustness of our approach across different data regimes, especially when individual signals become weak or insufficient.