Graph-SNMPBB: Projected Gradient for Graph Clustering

• The paper introduces Graph-SNMPBB, which combines a two-factor penalty formulation with graph Laplacian regularization to enhance cluster separability and achieve fast convergence.
• It employs a nonmonotone Barzilai–Borwein projected gradient scheme with two-stage updates, leveraging low-cost curvature estimates and sparse matrix operations.
• Experimental results show that Graph-SNMPBB attains competitive or superior clustering accuracy on benchmark datasets with significantly lower runtime compared to baselines.

Graph-SNMPBB is a projected-gradient algorithm for symmetric nonnegative matrix factorization tailored to graph clustering. It extends SNMPBB by combining a penalty-based two-factor symmetric factorization with graph Laplacian regularization, and it uses a nonmonotone Barzilai–Borwein projected gradient scheme with a two-stage projected update. In the reported formulation, the method is designed for an undirected, nonnegative, symmetric similarity or adjacency matrix and is intended to improve cluster separability while preserving nonnegativity and achieving fast practical convergence. The associated analysis proves global convergence to first-order stationary points under mild assumptions, and the reported experiments show competitive or superior clustering accuracy to strong baselines, including SymANLS, at substantially lower runtime on several benchmarks (Swart et al., 1 Jun 2026).

1. Problem setting and mathematical formulation

Graph-SNMPBB is posed for an input graph represented by an undirected, nonnegative, symmetric matrix $A \in \mathbb{R}^{n \times n}$. In the paper’s notation, $A$ is denoted $V$; the presentation here follows the equivalent notation $A \equiv V$. The degree matrix is $D = \operatorname{diag}(d_1,\ldots,d_n)$, where $d_i = \sum_j A_{ij}$, and the graph regularizer is the normalized symmetric Laplacian

$L_{\mathrm{sym}} = I - D^{-1/2} A D^{-1/2},$

denoted $L$. The target rank $r$ corresponds to the number of clusters. The factor variables are $W \in \mathbb{R}^{n \times r}_+$ and $A$0 (Swart et al., 1 Jun 2026).

Rather than enforcing a single-factor symmetric model directly, Graph-SNMPBB uses a two-factor penalty formulation: $A$1 Here $A$2 penalizes deviation from symmetry by coupling $A$3 and $A$4, while $A$5 imposes graph smoothness through the Laplacian term. The paper contrasts this with the more common single-variable regularized SymNMF objective

$A$6

but its algorithmic development and convergence analysis are based on the two-factor penalty model.

This formulation is central to the identity of Graph-SNMPBB. A common misconception is to view it as merely SymNMF with an added Laplacian penalty. The reported method is more specific: it relies on the penalty coupling $A$7 rather than imposing exact symmetry during optimization. The paper argues that this avoids the slow behavior often observed with vanilla symmetric one-factor projected gradient and gives SNMPBB and Graph-SNMPBB better curvature scaling for updates.

2. Differential structure, feasibility, and stationarity

For fixed $A$8 and $A$9, the gradients of the Graph-SNMPBB objective are

$V$0

and

$V$1

Because $V$2 is symmetric, the derivative of $V$3 with respect to $V$4 is $V$5 (Swart et al., 1 Jun 2026).

The nonnegativity constraints are characterized by KKT conditions in componentwise form. Primal feasibility requires $V$6 and $V$7. Dual feasibility is expressed through the gradients on active constraints in the sense of complementary slackness, and the complementary slackness relations are

$V$8

where $V$9 denotes the Hadamard product.

A practically useful stationarity measure is the projected-gradient mapping

$A \equiv V$0

with elementwise projection $A \equiv V$1. Stationarity is equivalent to $A \equiv V$2 and $A \equiv V$3 for any $A \equiv V$4. This projected-gradient viewpoint is what the algorithm uses operationally for both updates and stopping criteria.

The Laplacian term has a specific clustering interpretation. The trace penalty $A \equiv V$5 encourages smoothness over the graph, so nodes connected by edges with large $A \equiv V$6 have similar encodings in $A \equiv V$7, and through the symmetry penalty, in $A \equiv V$8 as well. This is why the method is presented as a graph clustering algorithm rather than only a matrix factorization routine.

3. Nonmonotone Barzilai–Borwein projected-gradient scheme

Graph-SNMPBB applies the SNMPBB update separately to the $A \equiv V$9-block and the $D = \operatorname{diag}(d_1,\ldots,d_n)$0-block, alternating between them. Each block update uses a two-stage projected step together with a nonmonotone Armijo line search (Swart et al., 1 Jun 2026).

For a $D = \operatorname{diag}(d_1,\ldots,d_n)$1-update with $D = \operatorname{diag}(d_1,\ldots,d_n)$2 fixed, the first projected gradient trial is

$D = \operatorname{diag}(d_1,\ldots,d_n)$3

The scale $D = \operatorname{diag}(d_1,\ldots,d_n)$4 is a Lipschitz step for the $D = \operatorname{diag}(d_1,\ldots,d_n)$5-subproblem, and $D = \operatorname{diag}(d_1,\ldots,d_n)$6 is computed from the spectral norm of the $D = \operatorname{diag}(d_1,\ldots,d_n)$7 matrix $D = \operatorname{diag}(d_1,\ldots,d_n)$8.

The Barzilai–Borwein curvature information is then formed from

$D = \operatorname{diag}(d_1,\ldots,d_n)$9

The BB1 and BB2 steps are

$d_i = \sum_j A_{ij}$0

using the Frobenius inner product for matrices. The paper uses the bounded BB1 variant,

$d_i = \sum_j A_{ij}$1

The second-stage projected direction is

$d_i = \sum_j A_{ij}$2

Acceptance is controlled by a Grippo–Lucidi style nonmonotone Armijo condition. The reference sequence is

$d_i = \sum_j A_{ij}$3

Starting from $d_i = \sum_j A_{ij}$4, backtracking is performed until

$d_i = \sum_j A_{ij}$5

The update is then

$d_i = \sum_j A_{ij}$6

with $d_i = \sum_j A_{ij}$7 used in practice.

The $d_i = \sum_j A_{ij}$8-update is completely analogous, with Lipschitz scale

$d_i = \sum_j A_{ij}$9

This reflects the left and right linear operators appearing in $L_{\mathrm{sym}} = I - D^{-1/2} A D^{-1/2},$0. The resulting alternation yields a block projected-gradient algorithm whose practicality comes from low-cost $L_{\mathrm{sym}} = I - D^{-1/2} A D^{-1/2},$1 curvature estimates and sparse matrix multiplications.

4. Initialization, stopping criteria, and implementation details

Preprocessing begins by constructing $L_{\mathrm{sym}} = I - D^{-1/2} A D^{-1/2},$2 from data, for example by a Gaussian kernel, retaining $L_{\mathrm{sym}} = I - D^{-1/2} A D^{-1/2},$3 nearest neighbors per node, and symmetrizing. The normalized Laplacian $L_{\mathrm{sym}} = I - D^{-1/2} A D^{-1/2},$4 is then computed (Swart et al., 1 Jun 2026).

The rank $L_{\mathrm{sym}} = I - D^{-1/2} A D^{-1/2},$5 is chosen as the number of clusters. For initialization, the paper uses random entries in $L_{\mathrm{sym}} = I - D^{-1/2} A D^{-1/2},$6 scaled by $L_{\mathrm{sym}} = I - D^{-1/2} A D^{-1/2},$7, where $L_{\mathrm{sym}} = I - D^{-1/2} A D^{-1/2},$8, and sets $L_{\mathrm{sym}} = I - D^{-1/2} A D^{-1/2},$9. It also notes that SVD-based warm starts using $L$0 and $L$1 can help on graph tasks, although they are more costly.

The reported hyperparameter guidance is explicit. The symmetry penalty is set as $L$2, and the graph regularization weight as $L$3. For nonmonotone line search, the paper uses $L$4, $L$5 with examples such as a constant $L$6, and a relaxation factor $L$7. The BB safeguards use bounds such as $L$8 and $L$9, with a reset when $r$0.

Stopping can be based on any of three criteria: $r$1 or

$r$2

or a KKT residual

$r$3

For application to a new graph dataset, the paper gives a concrete implementation pattern: build $r$4 via a KNN Gaussian kernel

$r$5

symmetrize via $r$6, choose $r$7 by a median distance heuristic, alternate the $r$8 and $r$9 updates while exploiting sparsity, and output labels by row-wise argmax of $W \in \mathbb{R}^{n \times r}_+$0 or by $W \in \mathbb{R}^{n \times r}_+$1-means on rows of $W \in \mathbb{R}^{n \times r}_+$2 after optional row normalization.

5. Complexity, convergence theory, and computational profile

For sparse $W \in \mathbb{R}^{n \times r}_+$3 and $W \in \mathbb{R}^{n \times r}_+$4, the computational structure is dominated by sparse matrix products and low-rank Gram matrices. The paper lists the following principal costs (Swart et al., 1 Jun 2026):

• $W \in \mathbb{R}^{n \times r}_+$5: $W \in \mathbb{R}^{n \times r}_+$6
• $W \in \mathbb{R}^{n \times r}_+$7: $W \in \mathbb{R}^{n \times r}_+$8
• $W \in \mathbb{R}^{n \times r}_+$9 with sparse $A$00: $A$01
• $A$02: $A$03
• $A$04: $A$05
• $A$06 with sparse $A$07: $A$08, approximately $A$09 for normalized $A$10

Accordingly, the per-inner-update cost is

$A$11

typically dominated by sparse multiplies on large sparse graphs. Memory usage includes storage of $A$12, $A$13, and sparse $A$14 and $A$15, with implementations often storing $A$16 as $A$17 by transposition for cache efficiency. The sparse products $A$18 and $A$19 parallelize naturally, while the $A$20 matrices $A$21 and $A$22 are comparatively cheap.

The convergence theory is blockwise. Under fixed $A$23, $A$24 is convex in $A$25 and $A$26 is Lipschitz with constant

$A$27

Under fixed $A$28, $A$29 is convex in $A$30 and $A$31 is Lipschitz with constant

$A$32

Using the nonmonotone Armijo sequence $A$33 and the estimate

$A$34

the paper shows that the line search terminates finitely, the sequence $A$35 is nonincreasing and bounded below, and $A$36. Hence any accumulation point $A$37 is first-order stationary. The paper situates this analysis within the nonmonotone spectral projected-gradient framework of Birgin–Martínez–Raydan and Li–Shi, and states that adding the Laplacian term preserves the convexity and Lipschitz properties needed for the $A$38-block and therefore extends the same global convergence conclusions to Graph-SNMPBB.

6. Clustering interpretation, empirical behavior, and limitations

The clustering output is derived from the soft assignment encoded in $A$39 or $A$40. A node $A$41 may be labeled by row-wise argmax,

$A$42

or by applying $A$43-means to the rows of $A$44, often after row normalization. The latter is described as often yielding more robust partitions. This places Graph-SNMPBB in close relation to spectral clustering objectives, but with a nonnegative factorization model and an explicit symmetry-penalty mechanism (Swart et al., 1 Jun 2026).

The empirical results reported in the paper distinguish between synthetic symmetric NMF experiments and graph clustering experiments. On synthetic dense data, SNMPBB attains residuals similar to SymANLS and SymNewton but with up to approximately $A$45 speedups, with advantages growing at higher ranks. For graph clustering, the paper reports that on the geometric rings (“bullseye”) dataset with $A$46 and $A$47, Graph-SNMPBB converges much faster than SNMPBB without Laplacian, SymANLS, PGD, and SymNewton, and more reliably recovers the nonlinear ring structure. On six real-world benchmarks—ORL faces, COIL20, Isolet1, MNIST subset, Reuters-21578 top-20 classes, and TDT2 top-20 classes—Graph-SNMPBB matches or exceeds SymANLS final accuracy on five of six datasets while converging substantially faster. The reported advantages are especially pronounced on high-dimensional sparse text datasets such as Reuters and TDT2, consistent with the $A$48 per-iteration cost relative to the $A$49 cost stated for SymANLS.

The paper also specifies recommended use-cases and failure modes. Graph-SNMPBB relies on a reasonable similarity graph: if the graph is very noisy or the KNN construction is poor, too small a $A$50 yields soft clusters, whereas too large a $A$51 can over-smooth and blur boundaries. Likewise, a very large $A$52 can dominate the objective and force $A$53 prematurely, while too small a $A$54 weakens the intended soft symmetry. On extremely sparse or flat-curvature problems, the backtracking line search may struggle; the recommended mitigations are bounded backtracking and conservative $A$55 bounds. The method is described as best suited to medium-to-large sparse graphs, moderate numbers of clusters, and datasets for which spectral clustering is already competitive.

A related variant, LAI-SNMPBB, replaces $A$56 by a low-rank approximation obtained through randomized sketching with thin QR and a small eigendecomposition. The paper states that the BB curvature information is unchanged by the approximation error because gradient differences cancel the constant bias term, and reports that LAI-SNMPBB outperforms LAI-SymPGNCG on 34 SuiteSparse matrices in both runtime and residual quality. This suggests that Graph-SNMPBB belongs to a broader algorithmic family in which nonmonotone projected BB updates are used as a scalable alternative to alternating least-squares and Newton-type methods, particularly when sparsity and graph structure can be exploited.