SNMPBB: Projected-Gradient for Symmetric NMF

• SNMPBB is a projected-gradient algorithm for symmetric NMF that reformulates the problem with a penalty-coupled two-factor model to improve optimization efficiency.
• It employs Barzilai–Borwein spectral steps and a two-stage projected update, achieving significant runtime improvements while maintaining convergence guarantees.
• Extensions like Graph-SNMPBB and LAI-SNMPBB adapt the method for graph clustering and large-scale inputs, offering practical benefits in clustering accuracy and scalability.

SNMPBB is a projected-gradient algorithm for symmetric nonnegative matrix factorization (Symmetric NMF) that combines a penalty-based two-factor formulation, Barzilai–Borwein spectral step sizes, and a nonmonotone projected line search. It was introduced as the first adaptation of nonmonotone projected Barzilai–Borwein methods to Symmetric NMF, with extensions to graph clustering and large-scale low-rank approximate-input settings through Graph-SNMPBB and LAI-SNMPBB, respectively (Swart et al., 1 Jun 2026). In this formulation, the longstanding view that projected-gradient methods are intrinsically too slow for Symmetric NMF is replaced by a more specific claim: the poor behavior of earlier methods was largely a consequence of the formulation and globalization strategy rather than of first-order methods per se.

1. Problem setting and reformulation

Symmetric NMF seeks a low-rank nonnegative factorization of a symmetric nonnegative matrix $V$ by approximating it as $WW^T$, with $W \in \mathbb{R}_+^{n\times r}$ and $r \ll n$. The baseline symmetric objective is

$f(W,W)=\frac{1}{2}\|V-WW^T\|^2.$

This objective is quartic in $W$, hence nonconvex and harder to optimize than the standard two-factor NMF objective. The setting is especially important when $V$ is a graph similarity or adjacency matrix, because the symmetric factorization aligns naturally with clustering structure (Swart et al., 1 Jun 2026).

SNMPBB does not optimize the quartic one-variable objective directly. Instead, it solves a penalty-coupled two-factor problem,

$$f(W,H;\lambda)=\frac{1}{2}\|V-WH\|^2+\frac{\lambda}{2}\|W-H^T\|^2,$$

subject to

$$W\in \mathbb{R}_+^{m\times r}, \qquad H\in \mathbb{R}_+^{r\times n},$$

with $m=n$ in the symmetric case. The parameter $WW^T$0 penalizes deviation from symmetry. This reformulation preserves the computational advantages of alternating optimization in two variables while driving $WW^T$1 and $WW^T$2 toward agreement.

A central design claim is that the naive projected-gradient extension that enforces $WW^T$3 after every step performs poorly in practice, possibly because the symmetric assignments reduce the degrees of freedom. SNMPBB replaces that rigid coupling with a soft quadratic penalty. This suggests that the main algorithmic innovation lies less in inventing a new descent primitive than in selecting a formulation in which projected spectral steps remain informative.

2. Core algorithmic structure

SNMPBB adapts the nonmonotone projected Barzilai–Borwein framework to the penalty objective by alternating between convex quadratic subproblems in $WW^T$4 and $WW^T$5. With $WW^T$6 fixed, the gradient with respect to $WW^T$7 is

$WW^T$8

and with $WW^T$9 fixed, the gradient with respect to $W \in \mathbb{R}_+^{n\times r}$0 is

$W \in \mathbb{R}_+^{n\times r}$1

Projection onto the nonnegative orthant is elementwise:

$W \in \mathbb{R}_+^{n\times r}$2

For each alternating subproblem, SNMPBB uses a two-stage projected update. First, from the current iterate $W \in \mathbb{R}_+^{n\times r}$3, it computes a projected gradient trial point

$W \in \mathbb{R}_+^{n\times r}$4

The local Lipschitz constant of $W \in \mathbb{R}_+^{n\times r}$5 is

$W \in \mathbb{R}_+^{n\times r}$6

and $W \in \mathbb{R}_+^{n\times r}$7 is omitted in practice because it is usually small relative to $W \in \mathbb{R}_+^{n\times r}$8.

Second, SNMPBB constructs a projected Barzilai–Borwein search direction

$W \in \mathbb{R}_+^{n\times r}$9

The Barzilai–Borwein step size uses Frobenius-inner-product secant information. With

$r \ll n$0

the spectral step is

$r \ll n$1

safeguarded by

$r \ll n$2

The paper characterizes this two-stage mechanism as yielding a “quasi-Newton-like update using only first-order information” (Swart et al., 1 Jun 2026).

3. Nonmonotone line search and stationarity

A defining feature of SNMPBB is that it is nonmonotone: it does not require $r \ll n$3 at every iteration. Instead, it uses a reference envelope

$r \ll n$4

with $r \ll n$5. When $r \ll n$6, this reduces to a monotone Armijo rule; when $r \ll n$7, temporary increases are allowed relative to the current function value as long as the step is acceptable with respect to the running envelope.

The line search seeks $r \ll n$8 such that

$r \ll n$9

where $f(W,W)=\frac{1}{2}\|V-WW^T\|^2.$0. Backtracking begins from $f(W,W)=\frac{1}{2}\|V-WW^T\|^2.$1. After acceptance, the update is

$f(W,W)=\frac{1}{2}\|V-WW^T\|^2.$2

with fixed relaxation factor $f(W,W)=\frac{1}{2}\|V-WW^T\|^2.$3; the paper reports that $f(W,W)=\frac{1}{2}\|V-WW^T\|^2.$4 worked best experimentally (Swart et al., 1 Jun 2026).

The stationarity analysis is expressed through the projected-gradient mapping

$f(W,W)=\frac{1}{2}\|V-WW^T\|^2.$5

Two identities are central:

$f(W,W)=\frac{1}{2}\|V-WW^T\|^2.$6

and

$f(W,W)=\frac{1}{2}\|V-WW^T\|^2.$7

These relations connect projected-step norms to first-order optimality under nonnegativity constraints. They also explain why the method can use projected spectral steps as both a search mechanism and a convergence diagnostic.

4. Convergence theory

The convergence theory proceeds by viewing SNMPBB as alternating minimization over two convex quadratic subproblems. When $f(W,W)=\frac{1}{2}\|V-WW^T\|^2.$8 is fixed, the $f(W,W)=\frac{1}{2}\|V-WW^T\|^2.$9-subproblem remains convex, and its gradient is Lipschitz continuous. For feasible $W$0,

$W$1

with

$W$2

An analogous statement holds for the $W$3-subproblem (Swart et al., 1 Jun 2026).

The nonmonotone reference sequence satisfies

$W$4

so $W$5 is nonincreasing and bounded below. The line search therefore terminates finitely, and accepted step sizes admit a uniform lower bound away from zero whenever the current iterate is not stationary. The main theoretical conclusion is

$W$6

and every accumulation point of the sequence generated by SNMPBB is a first-order stationary point.

The same reasoning extends to Graph-SNMPBB because the added Laplacian term preserves convexity in $W$7 and modifies the Lipschitz constant through the largest eigenvalue of $W$8. For LAI-SNMPBB, the notable theoretical point is different: the low-rank approximate-input model changes the gradient itself, but not the gradient difference used in the Barzilai–Borwein secant formula. Hence the curvature information required by the BB step is preserved exactly.

5. Extensions: graph regularization and low-rank approximate input

Graph-SNMPBB adds graph Laplacian regularization to the penalty objective:

$W$9

where the normalized graph Laplacian is

$V$0

and $V$1 is the diagonal degree matrix with

$V$2

The $V$3-gradient is unchanged, while the $V$4-gradient becomes

$V$5

The rationale is that plain Symmetric NMF often yields soft cluster assignments, whereas the Laplacian term encourages structure aligned with graph geometry. Hyperparameter guidance in the paper places $V$6 around $V$7 and $V$8 around $V$9 (Swart et al., 1 Jun 2026).

LAI-SNMPBB is designed for very large problems. A random sketch $$f(W,H;\lambda)=\frac{1}{2}\|V-WH\|^2+\frac{\lambda}{2}\|W-H^T\|^2,$$0 is used to form

$$f(W,H;\lambda)=\frac{1}{2}\|V-WH\|^2+\frac{\lambda}{2}\|W-H^T\|^2,$$1

followed by the small eigendecomposition

$$f(W,H;\lambda)=\frac{1}{2}\|V-WH\|^2+\frac{\lambda}{2}\|W-H^T\|^2,$$2

Defining

$$f(W,H;\lambda)=\frac{1}{2}\|V-WH\|^2+\frac{\lambda}{2}\|W-H^T\|^2,$$3

the matrix is approximated by

$$f(W,H;\lambda)=\frac{1}{2}\|V-WH\|^2+\frac{\lambda}{2}\|W-H^T\|^2,$$4

This replaces multiplication by $$f(W,H;\lambda)=\frac{1}{2}\|V-WH\|^2+\frac{\lambda}{2}\|W-H^T\|^2,$$5 with low-rank multiplication through $$f(W,H;\lambda)=\frac{1}{2}\|V-WH\|^2+\frac{\lambda}{2}\|W-H^T\|^2,$$6 and $$f(W,H;\lambda)=\frac{1}{2}\|V-WH\|^2+\frac{\lambda}{2}\|W-H^T\|^2,$$7, reducing the listed complexity from $$f(W,H;\lambda)=\frac{1}{2}\|V-WH\|^2+\frac{\lambda}{2}\|W-H^T\|^2,$$8 to $$f(W,H;\lambda)=\frac{1}{2}\|V-WH\|^2+\frac{\lambda}{2}\|W-H^T\|^2,$$9.

The approximation-quality statement carried into the paper is

$$W\in \mathbb{R}_+^{m\times r}, \qquad H\in \mathbb{R}_+^{r\times n},$$0

Most importantly for the algorithmic design, when $$W\in \mathbb{R}_+^{m\times r}, \qquad H\in \mathbb{R}_+^{r\times n},$$1 is fixed the approximate-input gradient satisfies

$$W\in \mathbb{R}_+^{m\times r}, \qquad H\in \mathbb{R}_+^{r\times n},$$2

so the bias is constant with respect to $$W\in \mathbb{R}_+^{m\times r}, \qquad H\in \mathbb{R}_+^{r\times n},$$3. Therefore the gradient difference used in the BB denominator obeys

$$W\in \mathbb{R}_+^{m\times r}, \qquad H\in \mathbb{R}_+^{r\times n},$$4

which is precisely the sense in which Barzilai–Borwein curvature information is preserved.

6. Empirical behavior, use cases, and limitations

The empirical results reported for SNMPBB overturn the earlier expectation that projected-gradient methods are intrinsically uncompetitive for Symmetric NMF. On synthetic data, the abstract reports that SNMPBB achieves a 6 times speedup over SymANLS for similar residuals, with the advantage growing at higher ranks. Representative dense synthetic results include $$W\in \mathbb{R}_+^{m\times r}, \qquad H\in \mathbb{R}_+^{r\times n},$$5, where SNMPBB attains runtime $$W\in \mathbb{R}_+^{m\times r}, \qquad H\in \mathbb{R}_+^{r\times n},$$6 and residual $$W\in \mathbb{R}_+^{m\times r}, \qquad H\in \mathbb{R}_+^{r\times n},$$7, compared with SymANLS runtime $$W\in \mathbb{R}_+^{m\times r}, \qquad H\in \mathbb{R}_+^{r\times n},$$8 and residual $$W\in \mathbb{R}_+^{m\times r}, \qquad H\in \mathbb{R}_+^{r\times n},$$9, and SymNewton runtime $m=n$0 and residual $m=n$1 (Swart et al., 1 Jun 2026).

For graph clustering, Graph-SNMPBB is reported to match or exceed SymANLS accuracy across six real-world benchmarks—ORL, COIL20, Isolet1, MNIST, Reuters-21578, and TDT2—while requiring substantially less time. On the bullseye geometric dataset, the graph-regularized variant reaches higher clustering accuracy faster than competing methods, whereas plain SNMPBB is materially weaker for clustering. This establishes an important methodological distinction: good symmetric factorization performance does not automatically imply good clustering performance, and the Laplacian term is operationally significant rather than cosmetic.

For large-scale approximate-input settings, LAI-SNMPBB outperforms LAI-SymPGNCG on 34 SuiteSparse matrices in both runtime and residual quality. At performance ratio $m=n$2, it obtains the better residual on about 70% of the problems. On the Web of Science dataset, the reported mean residuals are $m=n$3 for LAI-SNMPBB and $m=n$4 for LAI-SymPGNCG, with mean ARI $m=n$5 and $m=n$6, respectively; the practical advantage shown in the paper is primarily the time-to-solution profile rather than a uniformly dominant clustering metric.

The operational guidance is correspondingly differentiated. SNMPBB is the default choice for symmetric nonnegative low-rank approximation when the objective is factorization quality with low first-order cost. Graph-SNMPBB is the preferred variant for graph clustering. LAI-SNMPBB is intended for very large matrices where multiplication by the full input is the bottleneck. Initialization can use random $m=n$7 scaled by

$m=n$8

where $m=n$9 is the average value of entries of $WW^T$00, with $WW^T$01; an SVD-based initialization is also described for graph problems, though it was not used in the experiments because of cost.

The main limitations are parameter sensitivity and preprocessing dependence. The method depends on $WW^T$02, $WW^T$03, $WW^T$04, $WW^T$05, $WW^T$06, $WW^T$07, and the BB safeguards $WW^T$08. Graph-based performance depends heavily on similarity construction and $WW^T$09-nearest-neighbor sparsification. In LAI-SNMPBB, over-solving the inner projected-gradient subproblems can be harmful; the reported recommendation is to cap inner iterations at 3 to 5, which the paper interprets as a form of implicit regularization. These caveats do not undercut the central result, but they do locate the method within the usual trade space of first-order nonconvex optimization: low iteration cost, strong scalability, and favorable empirical speed, balanced against formulation and hyperparameter sensitivity.