Comparative Study of Numerical Methods for Sparse Linear Systems

Thesis Abstract

Incomputational science routinely solves large sparse linear systems arising from discretized partial differential equations and data-intensive applications; however, the relative performance of contemporary sparse solvers across diverse problem classes remains unresolved, leading to suboptimal method selection in practice. This study addresses this gap by conducting a systematic, theory-informed comparison of numerical methods for sparse linear systems, focusing on iterative solvers (GMRES, CG, BiCGSTAB) and direct solvers (LU with fill-reducing orderings, multifrontal methods) across structured and unstructured meshes, as well as real-world application matrices from engineering, physics, and data science. The aim is to identify robust, scalable strategies for solver selection and preconditioning that generalize across problem domains. Specific objectives are (i) to evaluate convergence behavior, time-to-solution, and memory footprints for each solver-preconditioner pairing on a benchmark collection of 1,200 sparse matrices with varying sparsity patterns and conditioning; (ii) to quantify the impact of matrix properties (symmetry, positive definiteness, spectral radius, condition number) on solver performance using regression analysis and ANOVA; (iii) to assess scalability on parallel architectures up to 256 cores, including weak and strong scaling analyses; (iv) to develop a decision framework informed by theoretical guarantees from Krylov subspace theory and preconditioning theory that prescribes solver choices for given matrix characteristics; and (v) to validate the framework on three representative application domains finite element discretizations of elasticity problems, reservoir simulations, and graph-based semi-supervised learning tasks. The methodology adopts a mixed-methods, confirmatory comparative design. The population consists of numerical linear systems extracted from open-source benchmarks and industrial collaborations, yielding 1,200 matrices with known metadata. A stratified sampling approach ensures representation across symmetry, definiteness, sparsity patterns, and conditioning ranges. Data collection employs standardized benchmarks on a high-performance computing cluster equipped with Intel Xeon processors and Infiniband interconnect, recording metrics including iteration counts, time-to-solution, memory usage, flop counts, communication overhead, and preconditioner setup times. Instruments include the SuiteSparse Matrix Collection, the PETSc and Trilinos solver libraries, and HPCToolkit for performance profiling. Validity and reliability are enhanced through cross-validation of solver configurations, reproducible random seeds for iterative processes, and sensitivity analyses on preconditioner tolerances. Analytical methods comprise multivariate regression to model time-to-solution as a function of matrix properties and solver settings, ANOVA to detect significant differences among solver families, and nonparametric tests where normality assumptions fail. A meta-analysis synthesizes results across problem classes. A theoretical analytic framework grounded in Krylov subspace theory, spectral properties, and preconditioning theory guides interpretation, with the development of a practical decision framework that maps matrix descriptors to recommended solver-preconditioner combinations. Ethical considerations address transparent reporting of computational resource usage and reproducibility. Expected findings indicate that (a) preconditioned GMRES with algebraic multigrid preconditioners yields robust performance for non-symmetric, poorly conditioned matrices, while conjugate gradient with incomplete Cholesky excels for symmetric positive definite systems; (b) direct solvers remain advantageous for small-to-moderate sized matrices and for matrices with clustered spectra, but their memory demands escalate for larger problems, where multifrontal methods exhibit superior scalability under appropriate fill-reducing orderings; (c) matrix characteristics such as non-normality and spectral conditioning significantly influence convergence rates and parallel efficiency; (d) the proposed decision framework achieves improved average time-to-solution by 15–28% across benchmarks compared with baseline heuristic choices. The study contributes to knowledge by delivering a unified, evidence-based solver selection taxonomy, quantifying the role of matrix properties in solver performance, and providing a scalable methodology for benchmarking sparse linear solvers on heterogeneous architectures. It also offers practical recommendations for software practitioners, including guideline matrices linking problem descriptors to solver configurations and preconditioner choices. The main conclusion is that no single solver universally dominates; instead, a principled, matrix-aware strategy markedly enhances performance and resource efficiency. Recommendations include adopting a matrix-property-driven solver portfolio, implementing adaptive preconditioner tuning, and integrating the decision framework into solver libraries to support automated, problem-aware selection in large-scale simulations and data-driven applications. Suggestions for further research include extending the framework to heterogeneous architectures (GPU-accelerated solvers), exploring machine-learning-based predictor models for solver selection, and broadening benchmarks to incorporate time-dependent and nonlinear systems.