krylov subspace approximation

This study establishes tight theoretical lower bounds on the number of matrix-vector products required to solve linear systems $Ax = b$. By employing minimax analysis, adversarial constructions, and spectral theory, the authors derive the first complexity lower bounds with explicit constants. Their results show that any algorithm permitted to access both $A$ and its transpose must perform at least $\Omega(\kappa \log(1/\varepsilon))$ operations, where $\kappa$ is the condition number and $\varepsilon$ the desired accuracy. In contrast, algorithms restricted to using only $A$ (one-sided methods) require $n$ iterations even for well-conditioned problems, matching known upper bounds. This work reveals a fundamental distinction between one-sided and two-sided methods in their dependence on accuracy and conditioning, and demonstrates that Krylov subspace methods such as conjugate gradient are optimal under specific settings.

To address slow convergence of Krylov solvers for large, sparse, ill-conditioned linear systems, this paper proposes an AI-driven framework for automatic parameter optimization of MCMC-based matrix inverse preconditioners. Our method innovatively integrates graph neural networks (GNNs) and Bayesian optimization: GNNs model preconditioner performance as a function of the coefficient matrix’s graph structure to construct an efficient surrogate model; Bayesian acquisition functions then guide parameter search, eliminating costly manual or grid-based tuning. Evaluated on unseen ill-conditioned systems, our approach achieves superior preconditioning with only 50% of the conventional computational budget, reducing Krylov iterations by approximately 10% on average—significantly improving generalizability and computational efficiency. To the best of our knowledge, this is the first work to jointly leverage GNNs and Bayesian optimization for learning MCMC preconditioner parameters, establishing a scalable, data-efficient paradigm for intelligent preconditioning of sparse linear systems.

This paper investigates theoretical lower bounds for trace estimation of matrix functions, focusing on query complexity and computational cost of the Hutchinson method and Block Krylov subspace techniques. By establishing a tight theoretical connection between the number of Block Krylov iterations and the degree of polynomial approximation to scalar functions, the work derives the first query-complexity lower bound for estimating (mathrm{tr}(W^{-p})) under the Wishart matrix model. Concurrently, it obtains matching upper bounds on the required Krylov iteration count for key matrix functions including (A^{-1/2}) and (A^{-1}). The results uncover a fundamental trade-off among polynomial approximation accuracy, number of random queries, and Krylov iteration cost. This provides the first rigorous complexity characterization of trace estimation algorithms grounded in random matrix theory, thereby filling a long-standing gap in the theoretical understanding of lower bounds for stochastic trace estimators.

Ill-conditioned linear systems frequently cause divergence of Krylov subspace iterative methods (e.g., CG, GMRES), severely limiting their practicality in large-scale scientific computing. To address this, we propose a general, scalable numerical stabilization framework that systematically enhances robustness against ill-conditioning by embedding condition-number-aware preconditioning and residual regularization directly into the Krylov iteration process. The framework preserves the original Krylov structure, is method-agnostic, and has been integrated into the SciPy solver library. Extensive experiments on synthetic benchmarks and real-world high-dimensional ill-conditioned systems—including discretized partial differential equations and inverse problems—demonstrate that our approach achieves stable convergence across all test cases, consistently outperforms the baseline algorithms in convergence rate, and incurs only controlled increases in memory footprint and computational cost. This work provides a broadly applicable, reliable solution for solving large-scale ill-conditioned linear systems.

Solving ill-conditioned linear systems—particularly overdetermined and underdetermined cases—remains challenging due to slow convergence of traditional Krylov methods (e.g., CG, GMRES) and their inability to exploit outlier singular value structures. Method: This paper proposes two novel iterative algorithms—Kaczmarz++ and CD++—featuring adaptive momentum acceleration, Tikhonov-regularized projections, singular-value-aware sampling, and an equation-block information reuse memoization mechanism. Contribution/Results: We theoretically establish that both methods efficiently capture and leverage outlier singular values, achieving Krylov-level convergence for both overdetermined and underdetermined systems within a unified framework—the first such result. Experiments demonstrate that Kaczmarz++ significantly outperforms GMRES and CG across diverse ill-conditioned systems; CD++, while matching CG/GMRES in arithmetic complexity for positive semidefinite systems, exhibits strong competitive performance in benchmark tests.

To address the high-throughput, numerical stability, and computational efficiency requirements for matrix exponential evaluation in generative AI, this paper proposes a novel adaptive Taylor series evaluation algorithm. Unlike classical approaches—such as Paterson–Stockmeyer or Padé approximants combined with scaling-and-squaring—the method jointly optimizes the Taylor expansion order and scaling factor to minimize floating-point operations under a prescribed relative error tolerance (<1e−12), while integrating dynamic error control within a scaling-evaluation-squaring framework. Theoretical analysis guarantees numerical stability and achieves a Pareto improvement in both accuracy and asymptotic complexity. Empirically, on large-scale generative tasks—including diffusion models and continuous-time flow matching—the algorithm achieves an average 2.3× speedup over state-of-the-art methods. The implementation is open-sourced and validated across mainstream GPU platforms.

This work addresses the limited adaptability of traditional kernel-based extended dynamic mode decomposition (kEDMD), which relies on manually prescribed kernel functions and hyperparameters. The authors propose a learnable, weighted multi-kernel kEDMD framework that, for the first time, enables end-to-end gradient-based optimization of both kernel functions and their parameters to automatically approximate the Koopman operator of nonlinear dynamical systems. Inspired by dictionary learning, the method constructs a flexible combination of kernels and employs weight analysis to prune redundant components, thereby balancing model expressiveness with parsimony. Experiments on benchmark systems—including the Duffing oscillator and the Kuramoto–Sivashinsky equation—demonstrate that the proposed approach effectively learns high-quality kernel representations, significantly improving the accuracy and robustness of Koopman operator approximation.

This work addresses the challenge of constructing globally consistent Koopman eigenfunction representations for continuous-time dynamical systems exhibiting singularities—such as multistability, limit cycles, or separatrices—when only sparse, local observations are available. Conventional approaches struggle to achieve this efficiently. Leveraging the algebraic structure that non-zero Koopman eigenfunctions form a multiplicative group, the authors propose generating an expanded feature space via polynomial combinations of a small set of principal eigenfunctions. They further introduce a cross-singularity matching and continuation strategy that substantially enriches the repertoire of usable eigenfunctions. This framework enables high-fidelity, globally coherent modeling of dynamics from sparse data and significantly enhances the representation of key observables in complex systems.

This work addresses the high computational cost and limited flexibility of traditional deformable body subspace simulation, which relies on sequential time integration and struggles to support optimization tasks under geometric variations. The authors propose a novel approach based on Dynamic Mode Decomposition (DMD) to construct a low-rank Koopman operator that enables direct future-state prediction via matrix operations, yielding an efficient model for deformable dynamics. The key innovation lies in achieving, for the first time, shared Koopman dynamics across varying shapes and mesh resolutions—overcoming the limitation of existing DMD methods that are confined to fixed geometry and discretization. The resulting framework supports trajectory prediction with logarithmic-linear time complexity, allowing large numbers of time steps to be skipped while preserving accuracy, thereby significantly accelerating simulation and enabling efficient control, initial condition estimation, and shape optimization.

This work proposes the first fast deterministic algorithm for linear systems with arbitrary rectangular displacement structures—including Toeplitz-like, Vandermonde-like, and Cauchy-like matrices—addressing the limitation of existing fastest methods that rely on randomization and lack deterministic guarantees. The approach reformulates the structured linear system as a univariate polynomial modular equation and solves it via three deterministic steps: computing a vector M-Padé approximant basis, eliminating extraneous variables by solving a linear system over the polynomial ring, and recovering the solution through a synchronized M-Padé approximation. The algorithm deterministically computes both the solution to the linear system and a basis for its nullspace in $\tilde{O}(\alpha^{\omega-1}(m+n))$ operations over the base field, offering broader applicability and stronger reliability than prior randomized methods.