Nonlinear spectral clustering with C++ GraphBLAS
arXiv:2605.26975v1 [cs.DC] 26 May 2026
Dimosthenis Pasadakis & Olaf Schenk & Verner Vlacic & Albert-Jan Yzelman
Abstract—Nonlinear reformulations of the spectral clustering method have gained a lot of recent attention due to their increased numerical benefits and their solid mathematical background. However, the estimation of the multiple nonlinear eigenvectors is associated with an increased computational cost. We present an implementation of a direct multiway spectral clustering algorithm in the p-norm, for p ∈ (1, 2], using a novel C++ GraphBLAS API. The key operations are expressed in linear algebraic terms and are executed over the resulting sparse matrices and dense vectors, parameterized in the algebra pertinent to the computation. We demonstrate the effectiveness and accuracy of our shared-memory algorithm on several artificial test cases. Our numerical examples and comparative results against competitive methods indicate that the proposed implementation attains high quality clusters in terms of the balanced graph cut metric. The strong scaling capabilities of our algorithm are showcased on a range of datasets with up to 8 million nodes and 48 million edges. Index Terms—Algebraic programming, C++ GraphBLAS, graph p-Laplacian, spectral clustering
an example of which is a generalized semiring under which an SpMV multiplication takes place. The recently introduced C++11 implementation of GraphBLAS [7] has showcased impressive results on the speed-up of algorithms based on SpMVs [8]. We express the key operations of the method introduced in [4] as well as the k-means discretization of the resulting eigenvectors in this C++ GraphBLAS API. This allows us to leverage its auto-parallelisation capabilities, and furnish the first, according to our best knowledge, p-norm spectral clustering algorithm applicable to large-scale data for shared-memory machines. II. A C++ G RAPH BLAS p- SPECTRAL CLUSTERING ALGORITHM
For an undirected weighted graph G(V, E, W) where V is the set of n nodes, E the set of edges, and W the weighted adjacency matrix, estimating a set of k p-eigenvectors on the Grassmann manifold Gr can be expressed as
I. I NTRODUCTION Spectral clustering is a popular community detection method k X n X wij |uℓi − uℓj |p that can be applied to any kind of data with a suitable similarity , p ∈ (1, 2]. (1) minimize Fp (U) = 2∥uℓ ∥pp U∈Gr(k,n) metric between them forming a graphical structure. At its core ℓ=1 i,j=1 lies the computation of the mutually orthogonal eigenvectors Let ℓ = 1, 2, . . . , k denote the eigenvector indices, and, at of the graph Laplacian, a symmetric and positive semi-definite the minimizer, the columns of U = (u , . . . , u ) approximate 1 k matrix, which are treated as the spectral coordinates of the the eigenvectors associated with the smallest k eigenvalues graph, and are subsequently discretized using distance based of the p-Laplacian operator ∆ . For i ∈ V the p-Laplacian p algorithms [1]. This eigenspectrum computation offers ample operator is defined as (∆ u) = P p j∈V wij ϕp (ui − uj ) , with i room for parallelization, with both shared and distributed ϕp : R →pRPbeing ϕp (x) = |x|p−1 sign(x), and the p-norm is memory implementations widely used [2]. Nonlinear variants n p ∥u∥p = p i=1 |ui | , for u ∈ R. of the method in the p-norm, for p ∈ (1, 2], that have We use the Riemannian optimization software package been proposed lead to a minimization of balanced graph cut ROPTLIB [9] to minimize (1) for progressively smaller values metrics, and an increase in the accuracy of the final clustering of p using Newton’s method on the Grassmann manifold, assignment [3]. Recently in [4], p-spectral clustering was cast where the solution of the linearized Newton subproblems as a nonlinear unconstrained optimization problem on the is handled by a truncated conjugate gradient scheme. The Grassmann manifold [5], by approximating the constraint for description of the optimization problem is accomplished p-orthogonality with an analogous one for 2-orthogonality. This in ROPTLIB by specifying the function EucGrad which approach is not applicable to large-scale data, due to the large computes the gradient of Fp (U) as well as the function number of multiplications of the graph adjacency matrix with EucHessianEta which computes η 7→ Hη = (Hℓ η ℓ )kℓ=1 the computed eigenvectors that are required for convergence, for arbitrary η ∈ Rk×n , where the collection of matrices especially as the value of p tends to 1. H1 , . . . , Hk ∈ Rn×n corresponds to the Hessian of Fp (U). For GraphBLAS is a standard [6] for expressing graph compuillustration, we include our C++ GraphBLAS implementation tations in the language of linear algebra. Its core concepts are of EucHessianEta in Algorithm 1. Here the subroutines (i) algebraic containers, which correspond to sparse matrices ROPTLIBtoGRB and GRBtoROPTLIB serve for I/O from and and vectors, (ii) algebraic operators, describing sparse matrixto the ROPTLIB data structures. The C++ GraphBLAS API vector (SpMV) multiplications, and (iii) algebraic relations, leverages the algebraic structure of the ring of real numbers to Dimosthenis Pasadakis, and Olaf Schenk are with the Advanced parallelize the SpMV operation (the grb::vxm primitive). Computing Laboratory at the Institute of Computing, Università della Svizzera italiana (USI), Lugano, Switzerland. email: {dimosthenis.pasadakis, olaf.schenk}@usi.ch. Verner Vlacic, and Albert-Jan Yzelman are with the Computing Systems Lab, Huawei Zurich Research Center, Switzerland. email: {verner.vlacic, albertjan.yzelman}@huawei.com.
III. N UMERICAL R ESULTS In order to demonstrate the effectiveness of the C++ GraphBLAS API for implementing the p-spectral clustering
Algorithm 1 The function EucHessianEta. η, a k × n matrix (D[ℓ])kℓ=1 , where each D[ℓ] = diag(Hℓ ) Input: (H[ℓ])kℓ=1 , where each H[ℓ] = diag(Hℓ ) − Hℓ Output: r, the result of η 7→ Hη 1: grb::Semiring<grb::operators::add<double>, grb::operators::mul<double>, grb::identities::zero, grb::identities::one> reals ring; 2: std::vector<grb::Vector<double>> grb eta, grb res; 3: grb::Vector<double> v, w; 4: ROPTLIBtoGRB(η, grb eta); 5: for ℓ = 1 to k do 6: grb::set(v, 0); 7: grb::vxm(v, grb eta[ℓ], H[ℓ], reals ring); 8: grb::eWiseApply(w, grb eta[ℓ]),D[ℓ], grb::operators::mul<double>()); 9: grb::eWiseApply(grb res[ℓ], w, v, grb::operators::subtract<double>()); 10: GRBtoROPTLIB(grb res, r); 11: return r
(a)
(b)
Fig. 1: Strong scaling of the C++ GraphBLAS components of the algorithm for the Delaunay graphs. a) Results for the mid-scale cases of node size n ∈ [216 , 219 ], b) Results for the large-scale cases of node size n ∈ [220 , 223 ]. The runtime is normalized versus single-thread execution.
smallest case (r = 16) was ∼ 300 sec, and that of the largest one (r = 23) ∼ 20 hours. A breakdown of the runtime shows that only the GraphBLAS components of the algorithm exhibit excellent weak scalability for the large-scale graphs. IV. C ONCLUSION & O UTLOOK
Method
Del. 16
Del. 17
Del. 18
Del. 19
Spec pMulti GrB-pGrass
0.129 −6.21% −8.12%
0.089 −3.11% −6.43%
0.062 −2.21% −4.56%
0.045 −2.45% −4.19%
TABLE I. Results in terms of the balanced graph cut metric RCut. We report the baseline (Spec) RCut, and in percentage the reduction of the cut that the methods pMulti and GrB-pGrass (ours) achieved.
method, in Section III-A we report on the quality of the graph cuts obtained, and in Section III-B we present its parallel performance. For our experiments we select 8 matrices from the SuiteSparse matrix collection [10] with an increasing number of nodes n = 2r and edges m ≈ 6 ∗ 2r , for r = 16, . . . , 23, corresponding to Delaunay triangulations in a unit square. A. Quality of graph cuts We identify four clusters Ci , i = 1, 2, 3, 4, and compute the value of the balanced graph cut metric P4 i ,Ci ) RCut(C1 , C2 , C3 , C4 ) = i=1 cut(C . The results for the |Ci | mid-scale cases with r = {16, 17, 18, 19} are summarized in Table I. We compare our method (GrB-pGrass) against traditional spectral clustering (Spec) [1] and against the first full eigenvector analysis of p-Laplacian leading to direct multiway clustering (pMulti) [3]. B. Parallel performance The strong scaling results of the developed algorithm are illustrated in Figure 1. In both plots, the dashed red line indicates the ideal scalability. We utilize up to 32 threads for the mid-scale cases (Figure 1a), and up to 88 threads for the large-scale Delaunay graphs with r = {20, 21, 22, 23} (Figure 1b). On average, the parallel execution of the algorithm is 5.5× faster than its sequential variant for the mid-scale tests, and 6.4× faster for the large-scale cases. The run-time of the
In this work, we have expressed they key operations of a multiway p-spectral clustering algorithm in the C++ GraphBLAS API. This enabled accurate parallel clustering of large-scale graphs on a shared-memory machine. We intend to further explore the potential gains of expressing graph partitioning and clustering algorithms in linear algebraic terms. ACKNOWLEDGMENT D.P. and O.S. acknowledge the support of the joint DFG 470857344 and SNSF - 204817 project. R EFERENCES [1] U. Luxburg, “A tutorial on spectral clustering,” Statistics and Computing, vol. 17, no. 4, p. 395–416, Dec. 2007. [2] S. T. Wierzchoń and M. A. Kłopotek, Spectral Clustering. Cham: Springer International Publishing, 2018, pp. 181–259. [3] D. Luo, H. Huang, C. Ding, and F. Nie, “On the eigenvectors of pLaplacian,” Machine Learning, vol. 81, no. 1, pp. 37–51, 2010. [4] D. Pasadakis, C. L. Alappat, O. Schenk, and G. Wellein, “Multiway p-spectral graph cuts on Grassmann manifolds,” Machine Learning, vol. 111, no. 2, pp. 791–829, Feb 2022. [5] A. Edelman, T. A. Arias, and S. T. Smith, “The geometry of algorithms with orthogonality constraints,” SIAM J. Matrix Anal. Appl., vol. 20, no. 2, p. 303–353, Apr. 1999. [6] J. Kepner, D. Bade, A. Buluc, J. Gilbert, T. Mattson, and H. Meyerhenke, “Graphs, matrices, and the GraphBLAS: Seven good reasons,” Procedia Computer Science, vol. 51, pp. 2453–2462, 2015. [7] A. N. Yzelman, D. Di Nardo, J. M. Nash, and W. J. Suijlen, “A C++ GraphBLAS: specification, implementation, parallelisation, and evaluation,” 2020, preprint. [Online]. Available: http://albert-jan.yzelman. net/PDFs/yzelman20.pdf [8] A. Scolari and A. Yzelman, “Effective implementation of the high performance conjugate gradient benchmark on GraphBLAS,” in 2023 IEEE International Parallel and Distributed Processing Symposium Workshops (IPDPSW). Los Alamitos, CA, USA: IEEE Computer Society, May 2023, pp. 216–225. [9] W. Huang, P.-A. Absil, K. A. Gallivan, and P. Hand, “Roptlib: An objectoriented C++ library for optimization on Riemannian manifolds,” ACM Trans. Math. Softw., vol. 44, no. 4, Jul. 2018. [10] T. A. Davis and Y. Hu, “The university of Florida sparse matrix collection,” ACM Trans. Math. Softw., vol. 38, no. 1, dec 2011.