Keywords: matrix condition number, graph neural networks, AI for numerical methods

1 Introduction↩︎

The paper is concerned with solving a general condition number of a sparse matrix \(A\) of arbitrary size. The condition number \(\kappa(\mathbf{A})\) of a matrix \(\mathbf{A} \in \mathbb{R}^{n \times n}\) quantifies the sensitivity of the solution of a linear system to perturbations in the input data. For a nonsingular matrix \(\mathbf{A}\), the condition number with respect to the \(p\)-norm is defined by \[\kappa_p(\mathbf{A}) = \|\mathbf{A}\|_p \, \|\mathbf{A}^{-1}\|_p. \label{eq:condition95number}\tag{1}\] where the applications of interest lie in \(p=1, 2, +\infty\). Calculating such exact condition number is generally prohibitively expensive for large and sparse matrices. In practice, we rarely do this because it often involves costly steps, such as finding the largest and smallest singular values or computing the matrix inverse, which do not scale well for large, sparse problems. For instance, the spectral condition number \(\kappa_2(\mathbf{A})\) depends on the extreme singular values of \(\mathbf{A}\). In dense linear algebra, these may be obtained by a full singular value decomposition (SVD), but this takes \(\mathcal{O}(n^3)\) floating-point operations. Similarly, forming \(\mathbf{A}^{-1}\) directly also takes \(\mathcal{O}(n^3)\) work in the dense case and is usually avoided. With sparse matrices, things get even trickier because sparse factorizations can lead to significant fill-in.

To address the high computational cost, one often resorts to estimation methods. As referenced in [1], condition number estimation is the process of calculating an inexpensive estimate of a condition number, where the inexpensive usually means that the cost is an order of magnitude less than the exact method. For spectral condition number estimation, iterative methods such as power iteration, Lanczos, Arnoldi, or Golub–Kahan bidiagonalization reduce the cost substantially, typically requiring \(\mathcal{O}(k \cdot \mathrm{nnz}(\mathbf{A}))\) work for \(k\) iterations, up to orthogonalization costs. For decades, efficient condition-number estimation has therefore been recognized as essential in assessing the stability of linear systems and least-squares problems, particularly in conjunction with matrix factorizations such as LU with partial pivoting; early practical strategies focused on inexpensive post-factorization estimates that avoid explicit inversion [2], [3].

For large sparse matrices, the exact computation of the condition number is often impractical because evaluating \(\|\mathbf{A}^{-1}\|_1\) exactly would, in the most direct approach, require either forming \(\mathbf{A}^{-1}\) explicitly or solving \(n\) linear systems to recover its columns. By contrast, the matrix norm \(\|\mathbf{A}\|_1\) is inexpensive to compute: \[\|\mathbf{A}\|_1 = \max_{1 \le j \le n} \sum_{i=1}^n |a_{ij}|,\] which can be evaluated in \(\mathcal{O}(\mathrm{nnz}(\mathbf{A}))\) time for a sparse matrix.

Accordingly, the main difficulty of estimating \(\kappa_1(\mathbf{A}) = \|\mathbf{A}\|_1 \, \|\mathbf{A}^{-1}\|_1\) lies in estimating \(\|\mathbf{A}^{-1}\|_1\). For matrices small enough to be treated by dense linear algebra routines, this quantity may be computed directly from a dense factorization. For larger sparse problems, a standard alternative is to estimate \(\|\mathbf{A}^{-1}\|_1\) by the Hager–Higham iterative algorithm [2], [3]. This procedure avoids explicit inversion by solving a sequence of linear systems \[\mathbf{A} x^{(k)} = \mathbf{b}^{(k)},\] where the right-hand sides \(\mathbf{b}^{(k)}\) are chosen adaptively from the sign pattern of intermediate iterates so as to approximate a vector attaining the maximum in the induced \(1\)-norm. In practice, a sparse LU factorization of \(\mathbf{A}\) is computed once and then reused for all subsequent triangular solves, making the estimator substantially cheaper than explicit inversion.

However, the aforementioned approximate method for condition number estimation mentioned above is limited to the 1-norm case. For the 2-norm (spectral) condition number, which requires estimating the ratio of the largest to smallest singular values, following \[\kappa_2(\mathbf{A}) := \|\mathbf{A}\|_2 \, \|\mathbf{A}^{-1}\|_2 = \frac{\sigma_{\max}(\mathbf{A})}{\sigma_{\min}(\mathbf{A})},\] where one often uses iterative methods such as the power method or Lanczos algorithm (or its bidiagonalization variant) to approximate the extreme singular values [4] and the main challenge is to calculate the smalles singular value \(\sigma_{\min}\). However, direct computation through a full singular value decomposition is generally infeasible for truly large sparse matrices. Both \(\|A\|_2\) and \(\|A^{-1}\|_2\) can be estimated via the power method [5]. An alternative practical way is to exploit Krylov subspace projections, which reduce the high-dimensional problem to a much smaller one while still capturing the behavior of the extreme singular values with good accuracy after a modest number of steps.

In the general nonsymmetric setting, a natural tool for estimating the extreme singular values is Golub–Kahan bidiagonalization [4]. Starting from an initial unit vector, the process constructs orthonormal bases \[\mathbf{U}_{k+1} = [\mathbf{u}_1,\dots,\mathbf{u}_{k+1}] \in \mathbb{R}^{n \times (k+1)}, \qquad \mathbf{V}_{k} = [\mathbf{v}_1,\dots,\mathbf{v}_{k}] \in \mathbb{R}^{n \times k},\] together with a lower bidiagonal matrix \(\mathbf{B}_k \in \mathbb{R}^{(k+1)\times k}\) such that\[\mathbf{A}\mathbf{V}_k = \mathbf{U}_{k+1}\mathbf{B}_k, \qquad \mathbf{A}^\top \mathbf{U}_{k+1} = \mathbf{V}_k \mathbf{B}_k^\top + \alpha_{k+1}\mathbf{v}_{k+1}\mathbf{e}_{k+1}^\top, \label{eq:gk95relation}\tag{2}\] and the projected matrix \(\mathbf{B}_k\) has bidiagonal structure, \[\mathbf{B}_k = \begin{pmatrix} \alpha_1 & & & \\ \beta_2 & \alpha_2 & & \\ & \beta_3 & \ddots & \\ & & \ddots & \alpha_k \\ & & & \beta_{k+1} \end{pmatrix} \in \mathbb{R}^{(k+1)\times k},\] where the coefficients \(\alpha_j > 0\) and \(\beta_j \ge 0\) arise from the recurrence. Each iteration requires one multiplication by \(\mathbf{A}\) and one multiplication by \(\mathbf{A}^\top\), so the dominant cost is typically \(\mathcal{O}(\mathrm{nnz}(\mathbf{A}))\) per step, in addition to orthogonalization. Accordingly, the singular values of the small bidiagonal matrix \(\mathbf{B}_k\) provide Ritz-type approximations to the singular values of \(\mathbf{A}\). In particular, the largest singular value usually converges first, while the smallest singular value may require a somewhat larger subspace dimension, especially when it is poorly separated. Once \(k\) is sufficiently large, one obtains the estimate \[\widehat{\kappa}_2(\mathbf{A}) \approx \frac{\sigma_{\max}(\mathbf{B}_k)}{\sigma_{\min}(\mathbf{B}_k)}.\] In practice, convergence may be monitored by the stabilization of the extreme singular-value approximations or by residual-based stopping criteria. In the special case where \(\mathbf{A}\) is symmetric positive definite, the situation simplifies considerably by replacing \(\sigma_{\max}(\mathbf{B}_k)\) and \(\sigma_{\min}(\mathbf{B}_k)\) with \(\lambda_{\max}(\mathbf{A})\) and \(\lambda_{\min}(\mathbf{A})\), respectively. As such, the Lanczos algorithm may be applied directly to \(\mathbf{A}\), producing a symmetric tridiagonal matrix \(\mathbf{T}_k\) whose extreme eigenvalues approximate those of \(\mathbf{A}\). The spectral condition number is then estimated by the ratio of the largest to the smallest Ritz eigenvalue of \(\mathbf{T}_k\).

For a general nonsymmetric matrix, however, Arnoldi iteration applied directly to \(\mathbf{A}\) approximates eigenvalues rather than singular values, and these are not sufficient in general to estimate \(\kappa_2(\mathbf{A})\), since \[\|\mathbf{A}\|_2 \neq \max_i |\lambda_i(\mathbf{A})| \quad\text{and}\quad \|\mathbf{A}^{-1}\|_2 \neq \frac{1}{\min_i |\lambda_i(\mathbf{A})|}\] for a nonnormal matrix. Therefore, for the estimation of the spectral condition number in the nonsymmetric case, one typically works with singular-value-oriented procedures such as Golub–Kahan bidiagonalization, or alternatively with Krylov methods applied to the normal equations \(\mathbf{A}^\top \mathbf{A}\), bearing in mind that the latter squares the condition number and may worsen numerical sensitivity.

We propose a data-driven approach based on GNNs for condition number estimation for sparse systems that: (i) Extracts structural and numerical features of matrix data in \(\mathcal{O}(\text{nnz})\) time; (ii) Trains a deep neural network to predict \(\log_{10} \kappa(A)\); (iii) Achieves inference time scales insensitive to, or independent of matrix size. Additionally, a novel graph machine learning approach is used for rapidly estimating the arbitrary-norm condition number of large sparse matrices arising in scientific computing applications. Unlike traditional iterative methods such as the Hager-Higham method [2], [6], our method uses a graph neural network trained on a comprehensive set of matrix-theoretic features to achieve sub-millisecond inference with controllable accuracy. Numerical experiments demonstrate that our approach achieves a significant speedup over classical methods and sparse solvers while maintaining an acceptable relative error. To the best of our knowledge, this is the first work on leveraging graph learning techniques to address condition number estimation. We believe this work can offer unique insights into the application of AI for numerical methods. We demonstrate the effectiveness of our approach with a simple neural network architecture and challenging condition number estimation tasks. The rest of the paper is structured as follows; Section 2 discusses related work on condition number estimation; Section 3 formulates our methodologies, with a focus on the matrix features as part of graph construction, and model architecture and also provide theoretical support for the condition number estimation; Section 4 presents our data generation routines for training and analyzes the simulations on performance; Section 5 illustrates the limitations of our approach, which we leave as future work; Section 6 concludes the paper.

2 Related Work↩︎

Recent advances in artificial intelligence have reshaped the landscape of scientific computing, motivating the integration of learning-based models with simulation, numerical linear algebra, and high-performance computing workflows [7]. Deep learning methods have been increasingly explored for scientific computing tasks such as surrogate modeling, operator learning, and solver acceleration [8], [9]. Motivated by this broader trend, we study the use of graph neural networks to estimate matrix condition numbers for sparse linear systems.

Estimating the condition number is critical for a linear matrix solve. A well-conditioned linear system generally does not require preconditioning, whereas a poorly conditioned one often does. In addition, a well-conditioned linear system can often be addressed with reduced precision in certain steps to improve performance.

The foundational contribution to estimating the condition number in this line of work is due to Cline et al. [10], who introduced a simple yet effective estimator based on the solution of two carefully chosen linear systems. Their approach, which was incorporated into the LINPACK library, produces a lower bound on the one-norm condition number at the cost of only a pair of triangular solves. Although computationally attractive, the method may underestimate the true condition number by one or more orders of magnitude in certain pathological cases.

To estimate matrix condition number, one often employs the power method to compute the spectral norm of a matrix, as introduced by Boyd in 1974 [11]. Following this, a significant improvement appeared shortly thereafter when Hager [2] reformulated the one-norm estimation problem as a convex optimization task. By interpreting \(\|\mathbf{A}^{-1}\|_1\) as the maximum of \(\|A^{-1}x\|_1\) over the unit ball in the one-norm and deriving an efficient iterative procedure (a specialized variant of the simplex method), Hager’s algorithm typically produces substantially tighter and more reliable estimates than the earlier heuristic, while still requiring only a modest number of matrix-vector products after factorization.

Higham [6] later provided a comprehensive survey of condition-number estimation specifically for triangular matrices—the form that naturally arises after LU factorization. The survey systematically compares growth-factor bounds, residual-based techniques, and inverse iteration approaches, offering both theoretical insights and practical recommendations regarding the trade-off between accuracy and computational cost. This work serves as a unifying reference and highlights subtle numerical issues that must be addressed when implementing estimators for upper- or lower-triangular factors.

Building directly on Hager’s optimization framework, Higham [3] subsequently proposed a refined one-norm estimator that addresses several practical limitations of the original algorithm. In particular, Higham provided a rigorous convergence analysis, extended the method to complex matrices, introduced safeguards against premature termination, and developed improved stopping criteria that enhance robustness without sacrificing efficiency. These theoretical and algorithmic refinements established a more reliable foundation for condition-number estimation in floating-point arithmetic.

For practical use, Higham [12] later released portable Fortran 77 implementations (subsequently published as Algorithm 674 in [13]). These routines implement the improved one-norm estimator together with its application to condition-number computation, include careful handling of both real and complex arithmetic, and incorporate extensive testing and error-handling mechanisms. The routines were later adopted—with minor modifications—into LAPACK and remain the de facto standard in many modern numerical libraries, effectively bridging the gap between theoretical development and robust, widely used software.

For the 2-norm condition number estimation, Krylov subspace methods [14]–[16] provide a coherent and computationally efficient framework for condition-number estimation in large sparse problems. Their cost is dominated by sparse matrix–vector multiplications, orthogonalization, and the solution of small projected eigenvalue or singular-value problems, making them far more practical than explicit inversion or full dense decompositions. Based on a randomized method, [17] solves the spectral condition number problem of a real matrix via estimating the smallest singular value by solving the least -squares problem using LSQR [18].

There is existing work [19], [20] that focuses on using GNNs to generate and select preconditioners for sparse linear systems, achieving excellent performance. The improved performance by GNNs often requires well-designed handcrafted features. However, existing research does not provide theoretical insights into feature modeling and selection; rather, it focuses on interpretation based on experience. In Section 3, we offer theoretical insight why we choose these matrix features based on how it effects the matrix condition number range.

Research on predicting condition numbers using machine learning methods remains relatively limited. The only one that published which we are aware is [21], which proposes an approach for predicting the condition number of sparse matrices using support vector regression. Although this method achieves lower accuracy than direct computation or classical estimation techniques, it can still be sufficient for coarse classification tasks, such as determining whether a matrix is well-conditioned or ill-conditioned.

3 Methodology↩︎

Consider the supervised learning problem: given a training set \(\mathcal{D} = \{(\mathbf{A}_i, \kappa^i)\}_{i=1}^N\) where \(\mathbf{A}_i \in \mathbb{R}^{n_i \times n_i}\) are sparse matrices and \(\|\cdot\|\) denotes arbitrary norm by removing subsript and same for \(\kappa\), we seek to construct a rapid estimator by exploiting the factorization \[\kappa(\mathbf{A}) = \|\mathbf{A}\| \cdot \|\mathbf{A}^{-1}\|. \label{eq:cond95factorization}\tag{3}\]

Observing that \(\|\mathbf{A}\|\) can be computed exactly in \(\mathcal{O}(\mathrm{nnz}(\mathbf{A}))\) time, we reformulate the estimation task: rather than directly predicting \(\kappa(\mathbf{A})\), we learn a surrogate function \(g: \mathcal{M} \to \mathbb{R}_+\) such that \[\label{eq:learning95objective95hybrid} g(\mathbf{A}) \approx \|\mathbf{A}^{-1}\| \quad \text{or} \quad g(\mathbf{A}) \approx \kappa(\mathbf{A}) \quad \forall \mathbf{A} \in \mathcal{M},\tag{4}\] where \(\mathcal{M}\) denotes the space of sparse matrices arising in target applications. The condition number estimate is then recovered via the hybrid formula \[\hat{\kappa}(\mathbf{A}) = \|\mathbf{A}\| \cdot g(\mathbf{A}) \quad \text{or} \quad \hat{\kappa}(\mathbf{A}) =g(\mathbf{A}). \label{eq:hybrid95estimate}\tag{5}\]

Following [13], this scheme is also useful for solving Sylvester equations \(AX + XB = C\), where \(A\) is an \(m\) by \(m\) matrix and \(B\) is an \(n\) by \(n\) matrix. The linear equations may also be presented in the form of \(Px = c\), where \(P =I_n \otimes A + B^T \otimes I_m\) as a Kronecker sum. Following [22] and [23], the estimation of the norm of \(P^{-1}\) is beneficial for assessing the accuracy of a computed solution.

To stabilize numerical conditioning and improve learning dynamics, we adopt the logarithmic transformation \[\tilde{g}(\mathbf{A}) = \log_{10} g(\mathbf{A}), \label{eq:log95transform95hybrid}\tag{6}\] So the training target becomes \[\tilde{\nu}_i = \begin{cases} \log_{10} \|\mathbf{A}_i^{-1}\| = \log_{10} \kappa_i - \log_{10} \|\mathbf{A}_i\|, & \text{if } \hat{\kappa}(\mathbf{A}) = \|\mathbf{A}\| \cdot g(\mathbf{A}) \\[2mm] \log_{10} \|\mathbf{A}_i\| \|\mathbf{A}_i^{-1}\|, & \text{if } \hat{\kappa}(\mathbf{A}) = g(\mathbf{A}) \end{cases}\]

and minimize the empirical risk \[\label{eq:loss95hybrid} \mathcal{L}(\theta) = \frac{1}{N} \sum_{i=1}^N \left( \tilde{g}(\mathbf{A}_i; \theta) - \tilde{\nu}_i \right)^2,\tag{7}\] where \(\theta\) denotes the model parameters. At inference time, the condition number estimate is obtained by exponentiating and combining with the exact norm: \[\label{eq:inference95hybrid} \hat{\kappa}(\mathbf{A}) = \begin{cases} \|\mathbf{A}\| \cdot 10^{\tilde{g}(\mathbf{A}; \theta^*)}, & \text{if norm scaling is used}, \\[4pt] 10^{\tilde{g}(\mathbf{A}; \theta^*)}, & \text{otherwise}. \end{cases}\tag{8}\] where \(\theta^*\) are the optimized parameters. In the following, we refer to the two formulas in 8 as the prediction schemes, namely, Scheme 1 and Scheme 2. The Scheme 1 is used when the norm of matrix is provided, and can be used to estimate the norm of inverse as well. Both schemes can be applied to Demmel condition number estimation.

3.1 Matrix Feature Modeling↩︎

In this section, we discuss the details of obtaining our feature modeling for GNNs. The direct application of deep learning to matrix inputs faces two fundamental challenges: (i) variable matrix dimensions exclude fixed-size network architectures, and (ii) the computational cost of processing \(\mathcal{O}(n^2)\) or \(\mathcal{O}(\text{nnz}(\mathbf{A}))\) entries exceeds the cost of traditional methods. We circumvent these obstacles through a carefully designed feature extraction operator \(\Phi: \mathbb{R}^{n \times n} \to \mathbb{R}^d\) that maps matrices to fixed-dimensional feature vectors while preserving essential conditioning information.

We associate each sparse matrix \(\mathbf{A} \in \mathbb{R}^{n \times n}\) with a compact \(d\)-dimensional feature vector \(\boldsymbol{\phi}(\mathbf{A})\) that encodes a rich collection of structural, numerical, and spectral properties known to be predictive of matrix conditioning and the performance of iterative solvers. We write the mapping \(\boldsymbol{\phi} : \mathbb{R}^{n \times n} \to \mathbb{R}^{d}, \quad \mathbf{A} \mapsto \boldsymbol{\phi}(\mathbf{A})\).

Following this, we propose a linear complexity algorithm with respect to \(\mathrm{nnz}\) for feature extraction. Rather than using all available descriptors, we construct a reduced yet informative representation by selecting a subset of components that capture the dominant scale, extremal behavior, and sparsity characteristics of \(\mathbf{A}\). This design is motivated by the observation that extremal statistics often provide strong proxies for conditioning and worst-case solver behavior, while avoiding redundancy introduced by moment-based descriptors.

The resulting feature vector is formed by concatenating selected components from complementary groups of descriptors, all computable in \(\mathcal{O}(\mathrm{nnz}(\mathbf{A}) + n)\) time: \[\label{eq:phi95concat} \boldsymbol{\phi}(\mathbf{A}) = \boldsymbol{\phi}_{\mathrm{norm}}(\mathbf{A}) \Vert \boldsymbol{\phi}_{\mathrm{struc}}(\mathbf{A}) \Vert \boldsymbol{\phi}_{\mathrm{diag}}(\mathbf{A}) \Vert \boldsymbol{\phi}_{\mathrm{rsp}}(\mathbf{A}) \Vert \boldsymbol{\phi}_{\mathrm{nzval}}(\mathbf{A}) \Vert \boldsymbol{\phi}_{\mathrm{gers}}(\mathbf{A}).\tag{9}\]

The detail of the matrix features used in this paper follows Table ¿tbl:tab:matrix95features?. Here, \(\|\mathbf{A}\|_1 = \max_j \sum_i |a_{ij}|\), \(\|\mathbf{A}\|_\infty = \max_i \sum_j |a_{ij}|\), and \(\|\mathbf{A}\|_F = \sqrt{\sum_{i,j} a_{ij}^2}\). Moreover, \(\mathbf{d} = \mathrm{diag}(\mathbf{A})\), \(\rho_i = \mathrm{nnz}(\mathbf{A}_{i,:})\), \(\mathcal{S} = \{a_{ij} : a_{ij} \neq 0\}\), and \(g_i = \sum_{j \neq i} |a_{ij}|\). The constant \(\epsilon\) is a small number (e.g., \(10^{-10}\)) that protects against underflow.

The assembly of \(\boldsymbol{\phi}(\mathbf{A})\) as in 9 requires only arithmetic operations \(\mathcal{O}(\mathrm{nnz}(\mathbf{A}) + n)\), making feature extraction negligible compared to any subsequent condition-number estimation procedure. These statistics provide a fixed-dimensional input representation for the GNN and enable an end-to-end learning framework for condition-number prediction. Their relevance can also be justified analytically: several components of the feature vector yield rigorous lower or upper bounds for matrix norms, spectral quantities, and condition numbers under appropriate structural assumptions.

In the theoretical statements below, we restrict our attention to 1-norm and 2-norm condition number, explaining how the selected features contribute to the prediction of condition number. We use the unlogged quantities underlying the logarithmic features. The small stabilizing constant \(\epsilon\) is introduced only for numerical robustness in feature computation and is ignored in the analytical bounds. Equivalently, the results should be interpreted in terms of the exact quantities before applying the logarithmic transformation.

The theoretical analysis below focuses on the following 11-dimensional subvector, denoted by \(\widetilde{\boldsymbol{\phi}}(\mathbf{A})\in\mathbb{R}^{11}\): \[\label{eq:feature1} \begin{align} \widetilde{\boldsymbol{\phi}}(\mathbf{A})=\Bigl( &\log_{10}(\|\mathbf{A}\|_1+\epsilon),\; \log_{10}(\|\mathbf{A}\|_\infty+\epsilon),\; \log_{10}(\|\mathbf{A}\|_F+\epsilon), \\ &\log_{10}\Bigl(\frac{\|\mathbf{A}\|_1+\epsilon}{\|\mathbf{A}\|_\infty+\epsilon}\Bigr),\; \log_{10}(\min |\mathbf{d}|+\epsilon), \\ &\log_{10}(\max \mathbf{g}+\epsilon),\; \log_{10}(n+1),\; \log_{10}(\max |\mathbf{d}|+\epsilon), \\ &\log_{10}\Bigl(\frac{\max |\mathbf{d}|+\epsilon}{\min |\mathbf{d}|+\epsilon}\Bigr),\; \log_{10}(\max \boldsymbol{\rho}+1), \\ &\log_{10}(\max |\mathcal{S}|+\epsilon) \Bigr)^\top, \end{align}\tag{10}\] where \(\mathbf{d}=\mathrm{diag}(\mathbf{A})\), \(\rho_i=\mathrm{nnz}(\mathbf{A}_{i,:})\), \(\boldsymbol{\rho}=(\rho_1,\ldots,\rho_n)^\top\), \(\mathcal{S}=\{a_{ij}:a_{ij}\neq 0\}\), and \(g_i=\sum_{j\neq i}|a_{ij}|\).

This subvector contains all entries of \(\boldsymbol{\phi}_{\mathrm{norm}}\) and selected extremal statistics from the structural, diagonal, row-sparsity, nonzero-value, and Gershgorin feature groups. These quantities admit direct analytical interpretation: the norm features bound \(\|\mathbf{A}\|_2\) and yield a lower bound for \(\kappa_1(\mathbf{A})\); the diagonal and Gershgorin features provide spectral bounds under symmetry; and the diagonal range, row sparsity, and maximum nonzero magnitude yield explicit condition-number intervals under a computable dominance margin.

The dominance margin is critical for condition-number control. The following remark gives a simple example illustrating why such a margin is necessary.

The results above show that the selected 11-dimensional subvector \(\widetilde{\boldsymbol{\phi}}(\mathbf{A})\) has explicit analytical relevance to matrix conditioning: it provides norm bounds, spectral bounds under symmetry, and condition-number intervals under diagonal-dominance-type assumptions.

3.2 Graph Neural Network Architecture↩︎

To leverage the structural information encoded in the sparsity pattern of \(\mathbf{A}\), we represent each matrix as an attributed graph \[\mathcal{G} = (\mathcal{V}, \mathcal{E}, \mathbf{X}, \boldsymbol{\phi}),\] where:

• \(\mathcal{V} = \{1, \ldots, n\}\) is the set of nodes corresponding to matrix rows/columns;

• \(\mathcal{E} = \{(i,j) : a_{ij} \neq 0\}\) is the directed edge set induced by the nonzero pattern of \(\mathbf{A}\);

• \(\mathbf{X} \in \mathbb{R}^{n \times 2}\) is the node feature matrix with rows \[\label{eq:node95features} \mathbf{x}_i = \bigl( \log_{10}(|a_{ii}| + \epsilon),\; \log_{10}(\rho_i + 1) \bigr)^\top,\tag{22}\] where \(\rho_i = \mathrm{nnz}(\mathbf{A}_{i,:})\) denotes the row sparsity;

• \(\boldsymbol{\phi}(\mathbf{A}) \in \mathbb{R}^{d}\) is the global feature vector defined in 9 .

The graph convolutional layers operate on the unweighted connectivity pattern specified by \(\mathcal{E}\); that is, the message-passing graph uses the sparsity structure of \(\mathbf{A}\) but does not assign additional learned or matrix-value-dependent edge attributes. One may seek improved performance via including edge attributes \(e_{ij}=\log_{10}(|a_{ij}|+\epsilon)\) but in the our GCN implementation these attributes are not directly used by the message-passing layers. The GCN branch uses the sparsity pattern through edge index while numerical magnitude information is incorporated through node features and graph-level statistical features. Numerical magnitude information is incorporated through the node features in 22 and the global matrix descriptor \(\boldsymbol{\phi}(\mathbf{A})\).

Additionally, when required by Scheme 1 in the inference formula 8 , we store the exact matrix norm \(\|\mathbf{A}\|_p\) as graph-level metadata. In particular, this corresponds to \(\|\mathbf{A}\|_1\) for the 1-norm condition number and \(\|\mathbf{A}\|_2\) for the 2-norm condition number.

3.2.1 Message Passing Layers↩︎

We employ a \(K\)-layer graph convolutional network (GCN) [24] to learn node embeddings from the unweighted sparsity connectivity of \(\mathbf{A}\). The GCN branch uses the edge set \(\mathcal{E}\) induced by the nonzero pattern, while matrix-value information is incorporated through the node features in 22 and the global descriptor \(\boldsymbol{\phi}(\mathbf{A})\).

The initial node embedding is obtained by a linear encoder followed by a nonlinear activation: \[\mathbf{h}_i^{(0)} = \sigma\left( \mathbf{W}_{\mathrm{node}}\mathbf{x}_i + \mathbf{b}_{\mathrm{node}} \right), \label{eq:node95encoder}\tag{23}\] where \(\mathbf{W}_{\mathrm{node}}\in\mathbb{R}^{\ell\times 2}\), \(\mathbf{b}_{\mathrm{node}}\in\mathbb{R}^{\ell}\), \(\ell\) is the hidden dimension, and \(\sigma(\cdot)\) denotes the ReLU activation.

Following the standard GCN normalization, we augment the graph with self-loops: \[\tilde{\mathcal{E}} = \mathcal{E} \cup \{(i,i):i\in\mathcal{V}\}.\] Let \(\tilde{\mathcal{N}}^{-}(i) = \{j:(j,i)\in\tilde{\mathcal{E}}\}\) be the set of source nodes that send messages to node \(i\), and let \(\tilde{d}_i = |\tilde{\mathcal{N}}^{-}(i)|\) be the corresponding degree. The \(k\)-th GCN layer is then written as \[\mathbf{h}_i^{(k+1)} = \sigma\left( \sum_{j\in \tilde{\mathcal{N}}^{-}(i)} \frac{1}{\sqrt{\tilde{d}_i\tilde{d}_j}}\, \mathbf{W}^{(k)}\mathbf{h}_j^{(k)} + \mathbf{b}^{(k)} \right), \qquad k=0,\ldots,K-1, \label{eq:gcn95propagation}\tag{24}\] where \(\mathbf{W}^{(k)}\in\mathbb{R}^{\ell\times\ell}\) and \(\mathbf{b}^{(k)}\in\mathbb{R}^{\ell}\) are learnable parameters. This propagation rule corresponds to message passing over the unweighted graph induced by the sparsity pattern of \(\mathbf{A}\).

3.2.2 Graph-Level Aggregation and Prediction↩︎

The final node embeddings \(\{\mathbf{h}_i^{(K)}\}_{i=1}^n\) are mapped to a graph-level representation using permutation-invariant mean and max readouts: \[\mathbf{z}_{\mathrm{mean}} = \frac{1}{n} \sum_{i=1}^n \mathbf{h}_i^{(K)}, \qquad \mathbf{z}_{\max} = \max_{1\leq i\leq n} \mathbf{h}_i^{(K)}, \label{eq:graph95pooling}\tag{25}\] where the maximum is taken element-wise. In parallel, the global matrix descriptor \(\boldsymbol{\phi}(\mathbf{A})\) is encoded by a separate fully connected branch: \[\mathbf{u} = \sigma\left( \mathbf{W}_{\mathrm{global}}\boldsymbol{\phi}(\mathbf{A}) + \mathbf{b}_{\mathrm{global}} \right), \label{eq:global95encoder}\tag{26}\] with \(\mathbf{W}_{\mathrm{global}}\in\mathbb{R}^{\ell\times d}\) and \(\mathbf{b}_{\mathrm{global}}\in\mathbb{R}^{\ell}\).

The graph-level representation is formed by concatenating the local structural readouts and the encoded global descriptor: \[\mathbf{q}(\mathbf{A}) = \left[ \mathbf{z}_{\mathrm{mean}}; \mathbf{z}_{\max}; \mathbf{u} \right] \in\mathbb{R}^{3\ell}. \label{eq:combined95representation}\tag{27}\] This representation is processed by a multilayer prediction head: \[\begin{align} \mathbf{s}^{(0)} &= \mathbf{q}(\mathbf{A}),\\ \mathbf{s}^{(m)} &= \sigma\left( \mathrm{Drop}_q \left( \mathbf{W}_{\mathrm{out}}^{(m)} \mathbf{s}^{(m-1)} + \mathbf{b}_{\mathrm{out}}^{(m)} \right) \right), \qquad m=1,\ldots,L-1,\\ \tilde{g}(\mathbf{A};\theta) &= \mathbf{W}_{\mathrm{out}}^{(L)} \mathbf{s}^{(L-1)} + \mathbf{b}_{\mathrm{out}}^{(L)}. \end{align} \label{eq:gnn95predictor95general}\tag{28}\] Here \(\mathrm{Drop}_q\) denotes dropout with rate \(q\) [25].

For the prediction scheme 1, the scalar output \(\tilde{g}(\mathbf{A};\theta)\) is then used in the corresponding inference formula: \[\hat{\kappa}_p(\mathbf{A}) = \|\mathbf{A}\|_p \cdot 10^{\tilde{g}(\mathbf{A};\theta^*)}, \label{eq:final95estimate}\tag{29}\] where \(\|\mathbf{A}\|_p\) is supplied explicitly in 29 , thus the neural network is trained to approximate the remaining inverse-norm factor in logarithmic scale.

While for prediction scheme 2, the scalar output follows: \[\hat{\kappa}_p(\mathbf{A}) = 10^{\tilde{g}(\mathbf{A};\theta^*)}.\]

The resulting two-stream architecture combines a local GCN branch, which captures the sparsity connectivity of \(\mathbf{A}\), with a global MLP branch, which incorporates matrix-level norm, diagonal, magnitude, and sparsity statistics.

4 Numerical Experiments↩︎

We evaluate the performance of our method against standard approaches for condition number estimation. Unless otherwise noted, all methods are implemented in PyTorch. Experiments and timing measurements were performed on a server equipped with 1 TB of RAM, dual Intel Xeon Gold 6330 processors (56 cores and 112 threads per processor, 2.00 GHz base frequency), and four NVIDIA A100 GPUs (80 GB PCIe each). All computations were executed on the a single GPU.

While condition-number estimation for sparse matrices has well-established CPU-based approaches, to the best of our knowledge, no widely adopted general-purpose GPU routine is available for this task in mainstream sparse linear algebra libraries. In particular, the optimized solution in Python is lacking. The mainstream Python GPU libraries do not currently provide a dedicated sparse-matrix condition-number API. PyTorch exposes torch.linalg.cond for dense tensors as a direct call for calculating condition number, while CuPy provides sparse SVD and Krylov subroutines that can be used to construct condition-number estimators, but not a direct sparse condition-number routine. In PyTorch, as of early 2026, torch.linalg.cond is not documented as a sparse-specialized condition-number estimator. According to its documented behavior, it computes condition numbers using standard dense linear algebra operations: singular-value computations for the \(2\)-norm-related cases, and matrix norm plus matrix inverse computations for norms such as the Frobenius, nuclear, \(1\)-, and \(\infty\)-norms. Although CuPy ² provides cuSPARSE-backed sparse matrix operations, our preliminary experiments showed that it was consistently slower than PyTorch for the algorithms considered in this study, across all tested data scales. We therefore restrict the subsequent evaluation to PyTorch in 2.12.1 with CUDA12.6.

To ensure a comprehensive and fair evaluation, we consider several implementations of the competing approaches. For the computation of the 1-norm condition number, we consider the following methods: \(\tikz[baseline=(char.base)]{ \node[shape=circle,draw,inner sep=1pt] (char) {1}; }\) computing the exact value using torch.linalg.cond on GPU; \(\tikz[baseline=(char.base)]{ \node[shape=circle,draw,inner sep=1pt] (char) {2}; }\) using scipy.sparse.linalg.norm to compute the matrix 1-norm, scipy.sparse.linalg.splu to perform linear solves required for matrix inversion, and scipy.sparse.linalg.onenormest³ to estimate the 1-norm of the matrix inverse; \(\tikz[baseline=(char.base)]{ \node[shape=circle,draw,inner sep=1pt] (char) {3}; }\) applying the Hager–Higham algorithm with torch.linalg.lu_factor for LU solves, while keeping the remaining operations implemented in PyTorch on GPU. The PyTorch-based implementations for 1-norm condition number evaluation are performed on both the CPU and the GPU. Methods \(\tikz[baseline=(char.base)]{ \node[shape=circle,draw,inner sep=1pt] (char) {2}; }\) and \(\tikz[baseline=(char.base)]{ \node[shape=circle,draw,inner sep=1pt] (char) {3}; }\) are referred to as Higham-SciPy and Hager–Higham, respectively, in the following.

For the computation of the 2-norm condition number, we consider two approaches: \(\tikz[baseline=(char.base)]{ \node[shape=circle,draw,inner sep=1pt] (char) {4}; }\) computing the exact value using torch.linalg.cond on GPU; \(\tikz[baseline=(char.base)]{ \node[shape=circle,draw,inner sep=1pt] (char) {5}; }\) a GPU-based Golub–Kahan method, i.e., a condition number estimate based on Ritz singular values computed via Golub–Kahan bidiagonalization; \(\tikz[baseline=(char.base)]{ \node[shape=circle,draw,inner sep=1pt] (char) {6}; }\) SciPy’s sparse singular value decomposition (scipy.sparse.linalg.svds), which extracts the largest and smallest singular values of the sparse matrix. Methods \(\tikz[baseline=(char.base)]{ \node[shape=circle,draw,inner sep=1pt] (char) {1}; }\) and \(\tikz[baseline=(char.base)]{ \node[shape=circle,draw,inner sep=1pt] (char) {4}; }\) are referred to as the exact approaches. Our proposed approach is referred to as GNN in the following. Our Golub–Kahan method takes the Golub–Kahan Lanczos bidiagonalization for estimating the 2-norm condition number of a matrix. To enhance numerical stability, full reorthogonalization is applied during the iterations using the Modified Gram–Schmidt procedure. The sparse matrices are stored as PyTorch sparse tensors to enable GPU-accelerated computations.

4.1 Training Protocols↩︎

We minimize the loss function 7 augmented with \(L^2\) regularization: \[\label{eq:reg95loss} \mathcal{L}_\text{reg}(\theta) = \mathcal{L}(\theta) + \lambda \sum_{L} \|\mathbf{W}^{(L)}\|_F^2,\tag{30}\] with regularization parameter \(\lambda = 10^{-5}\). We use the Adam algorithm optimizer [26] with learning rate \(\alpha = 10^{-3}\) and default momentum parameters \(\beta_1 = 0.9\), \(\beta_2 = 0.999\). Training proceeds for up to \(100\) epochs with mini-batch size \(32\), employing early stopping with patience \(20\) based on validation set performance.

The learning rate is adaptively reduced by a factor of \(0.5\) when validation loss plateaus for \(10\) consecutive epochs, implementing the ReduceLROnPlateau schedule [27]. Gradient clipping with threshold \(\tau = 1.0\) is applied to prevent exploding gradients during training. Dropout regularization [25] with a rate of \(0.1\) is applied after the first two hidden layers of the prediction head to mitigate overfitting prone to model training.

For the GNN-based estimator, we report two timing metrics. GNN (total) denotes the full pipeline runtime, including conversion of the sparse matrix into the graph representation required by the network, exact computation of \(\|A\|_2\), and subsequent model inference. GNN (inference) denotes only the pure forward-pass time of the trained GNN, excluding preprocessing and norm computation.

4.2 Evaluation Metrics↩︎

For the reference result on the exact value of the condition number, we use the function torch.linalg.cond to compute the condition number \(\kappa_p(A) = \|\mathbf{A}\|_p \cdot \|\mathbf{A}^{-1}\|_p\) depending on the chosen norm. When \(p=2\) (i.e., 2-norm), the library employs singular value decomposition to factorize \(A\), and defines the 2-norm condition number as the ratio of the largest to the smallest singular value, i.e., \(\kappa_2(A) = \sigma_{\max} / \sigma_{\min}\). This approach delivers high computational stability and is thus the default and most commonly used method. The underlying computations rely on LAPACK routines for numerical linear algebra. For other norms \(p\), it requires explicitly computing the inverse of \(A\) and thus is generally less numerically stable than the 2-norm computation based on SVD.

We assess both accuracy and computational efficiency using the following metrics:

• For each test matrix \(\mathbf{A}\), let \(\kappa_1(\mathbf{A})\) denote the ground truth condition number and \(\hat{\kappa}_1(\mathbf{A})\) the estimate produced by a given method. We define Logarithmic relative error: \[\mathrm{LRE}(\mathbf{A}) = \frac{\left| \log_{10} \hat{\kappa}(\mathbf{A}) - \log_{10} \kappa(\mathbf{A}) \right|}{\left| \log_{10} \kappa(\mathbf{A}) \right|}. \label{eq:lre}\tag{31}\] This normalized variant accounts for the magnitude of the true condition number, with values reported as percentages. In practice, for zero-division protection, a small number is added, e.g. \(10"^{16}\), to the denominator.

  We aggregate these metrics across the test set via their mean and maximum: \[\overline{\mathrm{LRE}} = \frac{1}{N_{\mathrm{test}}} \sum_{i=1}^{N_{\mathrm{test}}} \mathrm{LRE}(\mathbf{A}_i), \quad \mathrm{LRE}_{\max} = \max_{1 \leq i \leq N_{\mathrm{test}}} \mathrm{LRE}(\mathbf{A}_i). \label{eq:aggregated95metrics}\tag{32}\]

  Unlike standard relative error in linear space, the logarithmic metrics 32 is well-suited to condition numbers spanning multiple orders of magnitude. For instance, a linear relative error of \(100\%\) for \(\kappa_1 = 10^6\) (\(\hat{\kappa}_1 = 2 \times 10^6\)) translates to \(\mathrm{LRE} \approx 30\%\), correctly reflecting that the estimate remains within the same order of magnitude.

• To evaluate computational efficiency, we measure wall-clock time for each method over \(R = 4\) independent runs per test matrix, reporting:

  1. Mean inference time: \[T_{\mu}(\mathbf{A}) = \mathrm{mean}\left\{ t_1(\mathbf{A}), \ldots, t_R(\mathbf{A}) \right\}, \label{eq:median95time}\tag{33}\] where \(t_r(\mathbf{A})\) is the time for the \(r\)-th run. The median is preferred over the mean to mitigate outliers from system-level noise. For task 1), the runtime for GNN also additionally includes the trivial computation of \(\|\mathbf{A}\|_1\).

  2. Speedup factor: \[S_{\mathcal{M}_1 \to \mathcal{M}_2}(\mathbf{A}) = \frac{T_{\mu}^{\mathcal{M}_1}(\mathbf{A})}{T_{\mu}^{\mathcal{M}_2}(\mathbf{A})}, \label{eq:speedup}\tag{34}\] quantifying how many times faster method \(\mathcal{M}_2\) is compared to exact method \(\mathcal{M}_1\) (performed by numpy.linalg.cond(A)). We primarily report \(S_{\mathrm{Exact} \to \mathrm{GNN}}\) and \(S_{\mathrm{Higham} \to \mathrm{GNN}}\).

  For the GNN, we decompose the total time as \[T_{\mathrm{GNN}}(\mathbf{A}) = \begin{cases} T_{\mathrm{feature}}(\mathbf{A}) + T_{\mathrm{inference}}(\mathbf{A}) + T_{\mathrm{norm}}(\mathbf{A}), & \text{if } \hat{\kappa}(\mathbf{A}) = \|\mathbf{A}\| \cdot g(\mathbf{A}), \\ T_{\mathrm{feature}}(\mathbf{A}) + T_{\mathrm{inference}}(\mathbf{A}), & \text{if } \hat{\kappa}(\mathbf{A}) = g(\mathbf{A}). \end{cases}\] where \(T_{\mathrm{feature}}\) is the time to construct the graph and extract features \(\boldsymbol{\phi}(\mathbf{A})\), \(T_{\mathrm{inference}}\) is the forward pass through the GNN, and \(T_{\mathrm{norm}}\) is the overhead of computing \(\|\mathbf{A}\|_1\) and applying 29 .

4.3 Numerical Experiments↩︎

Following the above formulation, we focus on the two schemes on prediction of (1) \(g(\mathbf{A}) \approx \|\mathbf{A}^{-1}\|_1\) and (2) \(g(\mathbf{A}) \approx \kappa(\mathbf{A})\), \(\forall \mathbf{A} \in \mathcal{M}\), respectively. To mitigate the effects of environmental variability, we independently reran all competing methods for each scheme. We run each competing algorithm four times while deprecating the first warm-up time, and get their average performance as recorded score.

4.3.1 Parameterized Matrix Families↩︎

To ensure the trained model generalizes across diverse sparse linear systems, we construct a heterogeneous training corpus consisting of matrices from five representative problem classes as follows.

• Partial Differential Equation Discretizations:

  • 2D Poisson Equation. We consider the Poisson problem \[-\Delta u = f\] on the unit square with homogeneous Dirichlet boundary conditions. The equation is discretized on a uniform \(m \times m\) grid using the standard 5-point finite difference stencil on interior grid points (\(n = m^2\) unknowns).

    The resulting matrix \(\mathbf{A} \in \mathbb{R}^{n \times n}\) is sparse, symmetric positive definite (SPD), and corresponds to nearest-neighbor coupling in the grid. Up to a constant scaling factor, its diagonal entries are equal to \(4\) and off-diagonal entries are \(-1\).

    Its condition number satisfies \(\kappa_2(\mathbf{A}) = \Theta(m^2) = \mathcal{O}(h^{-2})\) [28].

  • Anisotropic Diffusion Equation. We consider the strongly orthotropic anisotropic diffusion problem \[-\nabla \cdot (K \nabla u) = f, \qquad K = \operatorname{diag}(\varepsilon, 1),\] where \(\varepsilon \ll 1\) induces strong anisotropy along the coordinate axes (principal directions aligned with the grid). The operator is discretized using a standard 5-point finite difference scheme on a uniform grid.

    The resulting matrix remains sparse and SPD, while its conditioning deteriorates significantly as \(\varepsilon \to 0\). In our dataset, \(\varepsilon\) is sampled from a log-uniform distribution \[\varepsilon \sim 10^{\mathcal{U}(-8,-3)},\] producing highly anisotropic systems.

  • High-Contrast Diffusion Problems. We consider variable-coefficient diffusion operators of the form \[-\nabla \cdot (\kappa(x,y)\nabla u),\] where the coefficient \(\kappa(x,y)\) exhibits strong spatial variation.

    The operator is discretized on a uniform grid using a conservative 5-point scheme with harmonic averaging of \(\kappa\) at cell interfaces, ensuring symmetry and positive definiteness of the resulting matrix.

    At each grid point, coefficients are sampled as \(\kappa_i \in [1, 10^c]\) with \(c \sim \mathcal{U}(6,13)\), producing matrices with large coefficient contrast. A uniform scaling factor (e.g., \(h^{-2}\)) is omitted, as it does not affect the spectral properties or condition number of the system.

    These matrices are representative of challenging problems arising in heterogeneous media and are typically highly ill-conditioned.

  • 2D Convection–Diffusion Equation. We consider the steady-state convection–diffusion problem \[-\varepsilon \nabla^2 u + \boldsymbol{\beta} \cdot \nabla u = f, \qquad \boldsymbol{\beta} = (\beta_x, \beta_y),\] on the unit square with homogeneous Dirichlet boundary conditions.

    The diffusion term is discretized using the standard 5-point stencil, while the convection term is discretized using a first-order upwind scheme. The resulting matrix \(\mathbf{A} \in \mathbb{R}^{n \times n}\) is sparse with at most five nonzero entries per row.

    Due to the upwind discretization, the matrix is generally nonsymmetric whenever \(\boldsymbol{\beta} \neq \mathbf{0}\). The scheme improves stability in convection-dominated regimes while explicitly breaking symmetry.

• Synthetic Random Matrices:

  • Random Sparse SPD Matrices. We generate random sparse symmetric positive definite matrices by first sampling a sparse matrix \(\mathbf{B}\) with independent entries and symmetrizing it: \[\mathbf{A}_0 = \mathbf{B} + \mathbf{B}^\top.\] Positive definiteness is enforced by adding a sufficiently large diagonal shift: \[\mathbf{A} = \mathbf{A}_0 + \alpha \mathbf{I},\] where \(\alpha\) is chosen based on the row-wise absolute sums to ensure strict diagonal dominance.

    An additional small diagonal perturbation is applied to influence the spectrum. This procedure yields sparse SPD matrices spanning a broad range of conditioning, from moderately to highly ill-conditioned regimes.

  • Ill-Conditioned SPD Matrices via Diagonal Scaling. To further generate extremely ill-conditioned matrices, we apply diagonal congruence transformations of the form \[\mathbf{A} = \mathbf{S} \mathbf{D} \mathbf{S},\] where \(\mathbf{D}\) is a sparse SPD matrix constructed as above and \(\mathbf{S} = \operatorname{diag}(s_i)\) is a diagonal scaling matrix with logarithmically distributed entries.

    This transformation preserves symmetry and positive definiteness while significantly amplifying spectral variation, producing matrices with very large condition numbers.

  • Symmetric Tridiagonal Matrices. We include symmetric tridiagonal matrices \[\mathbf{A} = \operatorname{tridiag}(-\alpha, 2, -\alpha) \in \mathbb{R}^{n \times n},\] where \(\alpha \sim \mathcal{U}(0.1, 0.9)\). These matrices are strictly diagonally dominant and hence SPD. Their eigenvalues are given analytically by \[\lambda_k = 2 - 2\alpha \cos\left(\frac{k\pi}{n+1}\right), \quad k=1,\dots,n,\] facilitating direct validation of spectral estimation methods.

We limit the variation in our training to increments of 2,000 (i.e., from 1,000 to 3,000), because increasing the gap would necessarily require more data. However, in this paper, we focus on model performance when trained with a limited amount of data across different machines or GPU scales. To verify that our method is insensitive to discrepancies in matrix sizes between the training and test sets, we restrict the training matrices to dimensions in the range of 1,000 and 3,000, and the test matrices to 500 and 5,000. We shuffled the final collection with a fixed random seed and split it into training, validation, and testing sets. The matrices for training, validation, and testing are sampled with equal probability from each class. Our training set consists only of size \(1,000\) matrices and a validation set of size \(100\) with the same matrix dimension as the training set, and a testing set of size \(100\). The generation routines ensure that the matrices used for training and testing exhibit a wide range of condition numbers. The distribution of condition numbers for the training and testing set are shown in Figure 2. The statistical information of the generated matrices for training, validation, and testing is shown in Table 2.

For Scheme 1, where \(\hat{\kappa}(\mathbf{A}) = |\mathbf{A}| \cdot 10^{\tilde{g}(\mathbf{A}; \theta^)}\), the learning task becomes \(\tilde{g}(\mathbf{A}; \theta^) \approx \log_{10} |\mathbf{A}^{-1}|\). The runtime and accuracy of the competing methods are reported in Table 3. For both the 1-norm and 2-norm condition numbers, the GNN approach consistently outperforms the Hager–Higham methods and the exact approach in terms of the average runtime. The speedup of GNN over Hager-Higham approach is close to \(9\)–\(11\times\). For the 1-norm condition number, GNN achieves lower \(\overline{\mathrm{LRE}}\) and \(\mathrm{LRE}_{\max}\) than the Hager–Higham method by an order of magnitude. For the 2-norm condition number, GNN is approximately \(5\)–\(10\times\) faster than the Golub–Kahan method, while also achieving substantially smaller logarithmic relative errors. Across both norm types, GNN maintains \(\overline{\mathrm{LRE}} < 0.5\) for all test samples. The two implementations of Hager–Higham show a constructive profile, while the version deployed in SciPy, though obtaining a slightly smaller \(\mathrm{LRE}_{\max}\), achieves longer runtime and higher \(\overline{\mathrm{LRE}}\), which is also observed in Table 3-that empirical results of Scheme 2.

For Scheme 2, the task focuses on \(\hat{\kappa}(\mathbf{A}) = 10^{\tilde{g}(\mathbf{A}; \theta^)}\), the task becomes a straightforward precision where \(\tilde{g}(\mathbf{A}; \theta^) \approx \log_{10} \bigl(|\mathbf{A}|,|\mathbf{A}^{-1}|\bigr)\). The corresponding runtime and accuracy are reported in Table 3. Similar to the results for the 1-norm condition number in Scheme 1, GNN again achieves the best runtime performance. For 1-norm and 2-norm condition number prediction, GNN attains \(\overline{\mathrm{LRE}} < 0.5\) for all test samples, which is slightly better than the Hager–Higham method. For 2-norm condition number prediction, GNN achieves \(\overline{\mathrm{LRE}} < 0.5\) for \(100\%\) of the test samples, whereas the Golub–Kahan method achieves this threshold for only \(56\%\) of samples with \(20\) iterations and \(67\%\) with \(60\) iterations. All competing methods maintain \(\overline{\mathrm{LRE}} < 1\) across all test samples.

Comparing the empirical results of Scheme 1 and Scheme 2, GNN exhibits a relatively smaller LRE and less runtime in Scheme 1 for both 1-and 2-norm cases. However, the 2-norm of a matrix is also a tricky task for estimating so Scheme 2 is more practical for 2-norm condition number estimation. Across the tests, the dominant bottleneck in the GNN approach lies in its inference phase.

In both 1-norm prediction schemes on CPU, Higham-SciPy exhibits close runtime compared to Torch Cond and twice the runtime achieved by the Hager-Higham method. While on GPU, Torch Cond and Hager-Higham perform similarly in speed, while being much faster than Higham-Scipy, regardless of the parameter \(t\). Among the competing algorithms on CPU deployment, the Hager-Higham achieves the fastest speed, followed by the GNN approach. With respect to both CPU and GPU deployment, we found that GNN achieves a significant decrease in runtime compared to competing algorithms on GPU.

For the 2-norm prediction, it is observed that the Golub–Kahan method with increasing iterations achieves a significant increase in runtime, though with an observed error prone is decreased.

4.3.2 SuiteSparse Matrix Collection↩︎

We built our dataset from the SuiteSparse Matrix Collection. We focused on real, square, sparse matrices with dimensions between 2,000 and 10,000. After applying structural filters on size and sparsity, we converted each candidate matrix to compressed sparse row format.

To keep only numerically well-behaved matrices, we removed any matrix for which the condition number estimation failed, returned non-finite values, or exceeded \(10^{20}\). This eliminated cases of singular or extremely ill-conditioned data. Same as the last test, we shuffled the final collection with a fixed random seed and split it into a 70% training set and a 30% test set. Eventually, we trained on 210 matrices (using 30 matrices among them for validation) and testing on 70 matrices. The information of training, validation, and testing set are presented in Table 4. The distribution of condition numbers for the training and testing set are shown in Figure 3.

For the ground-truth 2-norm condition number, we computed the reference 1-norm condition number as \(\kappa_1(A)=\|A|_1\|A^{-1}\|_1\), where \(\|A\|_1\) was computed directly and \(\|A^{-1}\|_1)\) was obtained by solving \(AX=I\) and evaluating the 1-norm of the resulting inverse. For the ground-truth 2-norm condition number, we computed the reference 2-norm condition number as \(\kappa_2(A)=\sigma_{\max}(A)/\sigma_{\min}(A)\), where the singular values were obtained using dense singular value decomposition.

For the 1-norm condition number estimation, we compared our GNN estimators against several standard baselines: exact 1-norm computation in PyTorch, SciPy’s 1-norm estimator with different block sizes (\(t=2, t=10\)), and a PyTorch implementation of the Hager–Higham estimator (both CPU and GPU versions when available). The GNN estimators follows the same two schemes–one predicting \(\log_{10}\|A^{-1}\|_1\) and reconstructing \(\kappa_1(A)\) from the known \(\|A\|_1\), and another that directly predicts \(\log_{10}\kappa_1(A)\). For the 2-norm condition number estimation, we benchmarked the proposed GNN estimators against several representative 1-norm and 2-norm baselines as illustrated in the last experiment. In addition, we evaluated two GNN prediction Scheme: one predicting \(\log_{10}\|A^{-1}\|\) and reconstructing \(\kappa(A)\) from the known \(\|A\|\), and another that directly predicts \(\log_{10}\kappa(A)\).

The empirical results on CPU and GPU are illustrated in Table 5. The GNN achieves the highest speedup on both CPU and GPU for the 1-norm and 2-norm cases. In terms of estimation accuracy evaluated by \(\overline{\mathrm{LRE}}\), the GNN estimators achieve \(\overline{\mathrm{LRE}} < 1\) in more than 80% of cases and \(\overline{\mathrm{LRE}} < 0.5\) in over 70% of cases. Although these results fall behind the Hager-Higham and Higham-SciPy approaches on 1-norm case, they are slightly better than the Golub–Kahan approach and worse than the power method on 2-norm case. Nevertheless, the GNN retains clear advantages in computational speed. The degraded accuracy compared to our previous simulation may be attributed to the limited training data and the distribution shift between the training and test sets.

Besides, from Table 5, SciPy svds fails to achieve \(\overline{\mathrm{LRE}} < 1\) over 85%. It fails on the test mainly because it requires to compute the smallest singular value using \(\texttt{svds}(A,\;\texttt{which="SM"})\), which is much harder and less stable than computing the largest singular value. If the matrix is ill-conditioned, then the smallest singular value \(\sigma_{\min}(A)\) is very close to zero, so tiny numerical errors can dominate the result. In addition, SciPy’s sparse SVD method with the default ARPACK solver works through an eigenvalue problem related to \(A^{T}A \text{or} A^{T}\), whose condition number satisfies approximately \(\kappa(A^{T}A) = \kappa(A)^2\). Therefore, a matrix that is already badly conditioned becomes even more numerically difficult internally. As a result, the iterative solver may fail to converge, return an inaccurate estimate of \(\sigma_{\min}(A)\), or trigger an exception. Thus, the failure is strongly related to the condition number, but it may also reflect solver instability rather than the matrix being exactly singular.

5 Limitations↩︎

Our approach relies on the quality of the training data: the closer the distributions of the training and test sets, the better the accuracy of the model. Ideally, the test data is expected to follow a distribution similar to that used to train the model; however, this is not strictly required, as we have evaluated on the test of SuiteSparse Matrix Collection, where training and testing dataset are different, and our GNN model’s generalization ability is demonstrated. Future studies will investigate which factors (e.g., matrix sizes or matrix element magnitude) are most sensitive and contribute to performance loss on testing data. In this paper, we focus on the formulation of our approach and demonstrate the capabilities of GNNs for condition number estimation. We have not optimized the model architecture, hyperparameters, or training routine, so readers may achieve further refined results by tuning the hyperparameters. Besides, as observed in the experiments, the dominant runtime for condition number prediction for GNN lies in the inference. Therefore, potential speedup can be achieved with ONNX [29].

6 Conclusion↩︎

We proposed a rapid condition number estimation approach using graph neural networks for sparse linear systems. We propose two prediction schemes, one by decomposing the condition number into an exactly computable forward norm and a neural-network-predicted inverse norm, and another by directly predicting the matrix condition number. To fully leverage the graph neural network power for matrix-based prediction, the neural network architecture combines graph convolutional layers to capture sparsity patterns with global feature encoders for statistical properties, enabling effective generalization. Future work will explore optimized GNN architectures and end-to-end feature engineering. Extensive experiments for the two schemes show that, with a limited amount of data, our method achieves orders of magnitude speedup over exact computation and compares favorably against the Hager-Higham and Golub–Kahan methods at a certain scale in terms of speed while maintaining competitive accuracy. Additionally, we study the selected matrix features theoretically and empirically for their contributions to model performance.

We believe that our work can be beneficial to existing work on precision tuning (e.g. [30]), since precision usage is highly dependent on condition number and a faster estimation of it can significantly boost the applications.

7 Acknowledegement↩︎

We would like to thank Cunyi Li for carefully proofreading the manuscript.

References↩︎

---

1. Supported by the Charles University Research Centre program No. UNCE/24/SCI/005 and the European Union (ERC, inEXASCALE, 101075632). Views and opinions expressed are those of the authors only and do not necessarily reflect those of the European Union or the European Research Council.↩︎

2. version 14.0.1, https://cupy.dev/↩︎

3. A Block Algorithm for Matrix 1-Norm Estimation, note that it is to compute a lower bound of the 1-norm of a sparse array (matrix).↩︎