The iSIRA method of Huang, Lin, Lu and Yau (J. Comput. Appl. Math. 346 (2019) 518-531) computes a few smallest positive eigenpairs of large sparse undirected graph Laplacians by integrating a trimming reduction that removes the singularity without factorization, low-rank deflation of converged eigenpairs, and inexact (iteratively solved) shift-invert residual Arnoldi with combinatorial-multigrid preconditioning. The paper closes by posing a fast eigensolver for Laplacians of directed graphs, where the symmetry, positive semi-definiteness, and simple nullspace structure underlying iSIRA all fail, together with a GPU-accelerated implementation; its body additionally leaves the shifted, deflated inner solves - including when the rank-one correction σn1^1^⊤\tfrac{\sigma}{n}\widehat{\mathbf{1}}\widehat{\mathbf{1}}^\top may be dropped - without rigorous criteria. The objective is to develop, analyze, and demonstrate at scale a fast, factorization-free eigensolver for the smallest nonzero eigenpairs of large directed graph Laplacians within the iSIRA framework, with robust, analyzable preconditioned inner solves and a scalable (GPU/parallel) implementation.

Background

The graph Laplacian of a simple, connected, undirected graph with nn vertices is L=D−AL=D-A, where DD is the diagonal matrix of vertex degrees and AA the adjacency matrix; LL is symmetric positive semi-definite, singular, and has a one-dimensional nullspace spanned by the all-ones vector 1n\mathbf{1}_n. Its smallest positive eigenvalues and eigenvectors encode community structure and drive spectral clustering, community detection, and network analysis. Computing a few of these eigenpairs for a large sparse LL (networks with millions of vertices from social media, transaction databases, or sensor systems) is the graph Laplacian eigenvalue problem (GLEP).

Classical eigensolvers for smallest eigenvalues use the shift-invert transformation: Arnoldi's method - a Krylov-subspace projection iteration that builds an orthonormal basis from repeated applications of (L−σI)−1(L-\sigma I)^{-1} - makes eigenvalues near the shift σ\sigma dominant, but each step requires solving a large sparse linear system, classically via an LU or Cholesky factorization. For big-data graph Laplacians this fails twice: LL is singular, so the zero-shift system is ill-posed; and real network matrices have wide-band, nonlocal connectivity that makes factorization memory- and time-prohibitive. Workarounds - a small diagonal (Tikhonov-type) perturbation, or nullspace deflation with exact inner solves - either perturb the target spectrum or still need the factorization; shift-invert iterations with singular operators also require purification, i.e. explicit removal of nullspace components from the iterates (Meerbergen and Spence, 1997).

iSIRA (Huang, Lin, Lu and Yau, 2019) integrates three ingredients into a factorization-free GLEP eigensolver: (i) trimming - the constrained singular system Lz=rL\mathbf{z}=\mathbf{r} with 1n⊤z=0\mathbf{1}_n^\top\mathbf{z}=0 is reduced exactly to an (n−1)×(n−1)(n-1)\times(n-1) nonsingular principal-submatrix system (delete one row and column) with an explicit reconstruction, giving nullspace-orthogonal search directions without factorizing LL; (ii) deflation - a low-rank term δVV⊤\delta VV^\top, with VV holding already-converged eigenvectors, shifts the converged eigenvalues by δ\delta while preserving the rest of the spectrum and the nullspace, so the iteration can move to the next-smallest eigenvalue; and (iii) the shift-invert residual Arnoldi method (SIRA) run inexactly - SIRA expands the search subspace by the residual of the current Ritz pair (an approximate eigenpair obtained from the projected problem) passed through the shifted operator, and the inner linear systems are solved only approximately by iteration, a framework analyzed by Jia and Li (2014). The trimmed inner systems have the form (L^−σI)+σn1^1^⊤+δV^V^⊤(\widehat{L}-\sigma I)+\tfrac{\sigma}{n}\widehat{\mathbf{1}}\widehat{\mathbf{1}}^\top+\delta\widehat{V}\widehat{V}^\top and are preconditioned by combinatorial multigrid (CMG) - near-linear-time preconditioners for symmetric diagonally dominant (SDD) systems built from low-stretch spanning trees and graph sparsification - with the deflation term inverted through the Sherman-Morrison-Woodbury low-rank-update formula; the rank-one term σn1^1^⊤\tfrac{\sigma}{n}\widehat{\mathbf{1}}\widehat{\mathbf{1}}^\top is dropped when nn is 'very large' purely as a heuristic, with no stated criteria. On networks from the Koblenz Network Collection (KONECT) with up to about 2.5 million vertices, iSIRA outperforms MATLAB eigs with diagonal perturbation and a null-space-deflation variant whenever an LU factorization is unavailable.

For directed graphs the natural Laplacian L=Dout−AL=D_{\mathrm{out}}-A (out-degree matrix minus adjacency) is nonsymmetric: eigenvalues may be complex, eigenvectors are not orthogonal, the matrix need not be SDD, and the nullspace is one-dimensional only for strongly connected digraphs - in general its dimension equals the number of sink strongly connected components. Perron-Frobenius theory still guarantees that for a strongly connected digraph all nonzero eigenvalues have positive real part, so 'smallest nonzero eigenvalues' remain meaningful (the directed spectral theory of Chung, 2006), and they govern consensus and synchronization dynamics, Markov-chain mixing, and spectral analysis of web, citation, and neural networks. Numerical work on directed Laplacians has so far targeted linear solves and pseudo-inverses (e.g., Boley, 2021) or symmetrized/Hermitian surrogates, not factorization-free computation of smallest nonzero eigenpairs.

Problem Statement

The iSIRA paper closes by posing, for Laplacians arising from directed graphs, the question of exploring new features and developing a fast eigensolver, noting that all properties and advantages of the undirected Laplacian used in that work no longer apply; it further announces the development of a GPU version of iSIRA, while its body leaves the inner-solve preconditioning analysis - robustness for nonzero shifts and multiple low-rank deflation terms, and when the rank-one correction can be dropped - without rigorous criteria. The problem is:

Develop, analyze, and demonstrate at scale a fast, factorization-free eigensolver for computing a few smallest nonzero eigenvalues and associated eigenvectors of large sparse Laplacian matrices arising from directed graphs, preserving the iSIRA design of inexact inner solves, explicit singularity handling, and deflation. The objective comprises three jointly linked components:

  1. Nonsymmetric theory. Replace the symmetry-based reductions of iSIRA - the trimming reduction of the constrained singular shifted-deflated system to a nonsingular solvable system, the deflation spectral-shift property, and the consistency/orthogonality conditions - by valid analogues for directed-graph Laplacians (e.g., L=Dout−AL=D_{\mathrm{out}}-A), covering kernels whose dimension is governed by the strongly-connected-component structure and possibly complex, clustered spectra near zero; and state, then consistently use, the ordering convention under which 'smallest nonzero eigenvalues' is defined (e.g., by real part or by modulus).

  2. Inner solves and preconditioning. Design and analyze scalable preconditioners for the trimmed, shifted, deflated inner systems (L^−σI)+σn1^1^⊤+δV^V^⊤(\widehat{L}-\sigma I)+\tfrac{\sigma}{n}\widehat{\mathbf{1}}\widehat{\mathbf{1}}^\top+\delta\widehat{V}\widehat{V}^\top and their nonsymmetric counterparts, robust for nonzero shifts σ\sigma and multiple low-rank deflation terms at very large nn, including a rigorous, proved criterion for when the rank-one correction σn1^1^⊤\tfrac{\sigma}{n}\widehat{\mathbf{1}}\widehat{\mathbf{1}}^\top may be neglected.

  3. Scale and implementation. Provide a parallel or GPU-accelerated implementation exploiting the concurrency of the inner iterative solves and the outer Arnoldi framework, evaluated on real large-scale directed networks against reasonable baselines (e.g., implicitly restarted shift-invert Arnoldi with regularization or symmetrization, and dense or factorization-based solvers where feasible), reporting eigenpair residuals, nullspace/purification accuracy, and runtime/memory.

The three components are jointly required for the one objective stated above; a GPU-accelerated implementation of the undirected iSIRA is the paper's explicitly announced immediate goal, but an undirected implementation alone does not resolve the directed-graph objective.

References

  1. Wei-Qiang Huang, Wen-Wei Lin, Henry Horng-Shing Lu, Shing-Tung Yau. iSIRA: Integrated shift-invert residual Arnoldi method for graph Laplacian matrices from big data. Journal of Computational and Applied Mathematics 346 (2019) 518-531. DOI: 10.1016/j.cam.2018.07.031, https://doi.org/10.1016/j.cam.2018.07.031 (preprint: arXiv:1710.07080, https://arxiv.org/abs/1710.07080)
  2. Daniel Boley. On fast computation of directed graph Laplacian pseudo-inverse. Linear Algebra and its Applications 623 (2021) 128-148. DOI: 10.1016/j.laa.2020.10.018, https://doi.org/10.1016/j.laa.2020.10.018
  3. Fan Chung. The Diameter and Laplacian Eigenvalues of Directed Graphs. The Electronic Journal of Combinatorics 13(1) (2006), article N4. DOI: 10.37236/1142, https://doi.org/10.37236/1142
  4. Zhongxiao Jia, Cen Li. Inner iterations in the shift-invert residual Arnoldi method and the Jacobi-Davidson method. Science China Mathematics 57(8) (2014) 1733-1752. DOI: 10.1007/s11425-014-4791-5, https://doi.org/10.1007/s11425-014-4791-5
  5. Karl Meerbergen, Alastair Spence. Implicitly restarted Arnoldi with purification for the shift-invert transformation. Mathematics of Computation 66(218) (1997) 667-689. DOI: 10.1090/s0025-5718-97-00844-2, https://doi.org/10.1090/s0025-5718-97-00844-2
  6. Jerome Kunegis. KONECT: The Koblenz Network Collection. Proceedings of the 22nd International Conference on World Wide Web (WWW Companion) (2013) 1343-1350. DOI: 10.1145/2487788.2488173, https://doi.org/10.1145/2487788.2488173