ABSTRACT We aim to accelerate the restarted generalized minimal residual (GMRES) method for the solutions of linear systems by combining two types of techniques. On the one hand, mixed precision GMRES algorithms, which use lower precision in certain steps of the inner cycles, offer significant reductions in computational and memory costs. On the other hand, augmented GMRES algorithms, which recycle information on the eigenvalues of the matrix between restarts by incorporating additional vectors into the Krylov basis, can significantly speed up the convergence. In this work, we investigate how to combine mixed precision and augmentation, in order to cumulate the reduced per‐iteration cost of the former with the reduced number of iterations of the latter. We first explore the GMRES with deflated restarting (GMRES‐DR) variant, which we demonstrate possesses limited mixed precision opportunities. Indeed, GMRES‐DR can exploit a preconditioner constructed in low precision, but requires a flexible paradigm to also apply it in low precision; moreover, the matrix–vector product and orthonormalization steps must both be kept in high precision, as otherwise the method stagnates at low accuracy. We explain that this is because GMRES‐DR is based on some algebraic simplifications that are only valid in exact arithmetic, but fail to hold in finite precision. This observation leads us to investigate a more general augmented GMRES framework (AugGMRES) that avoids making these simplifications. AugGMRES is much more resilient to the use of low precision, does not require a flexible paradigm, and successfully converges to high accuracy even when low precision is used for all inner operations. Our numerical experiments illustrate the robustness and fast convergence of AugGMRES on a range of sparse matrices, including ill‐conditioned real‐life ones.
ABSTRACT This paper focuses on approximating a quaternion matrix and its singular values via the quaternion randomized block Krylov iteration (qRBKI) algorithm. Specifically, for qRBKI with block size , we establish novel gap‐based error bounds for the rank‐ matrix approximation in both Frobenius and spectral norms, explicitly demonstrating that qRBKI outperforms the quaternion randomized subspace iteration (qRSI) in approximation quality. For , leveraging simulated starting block vectors and quaternion polynomials, we show that the approximation error depends on the minimal relative singular value gap under the ‐order neighbouring structure. The near‐zero‐gap‐induced dependence can be eliminated by perturbing the input matrix with a low‐intensity random quaternion Gaussian or uniform matrix. Finally, numerical experiments on synthetic data and color image applications validate our theoretical findings and illustrate the effectiveness of qRBKI.
ABSTRACT Reservoir computing is analyzed through the algebraic structure generated by finite‐width reservoirs on fixed input blocks. For linear and quadratic logistic activations, the induced input‐to‐state map admits an exact finite Volterra representation, yielding an explicit characterization of the associated hypothesis class. It is shown that reservoir width imposes rank constraints on the admissible Volterra kernels, which in the quadratic case reduce to classical matrix‐rank limitations and lead to approximation estimates in terms of spectral decay and best low‐rank approximation error. A Volterra‐Barron complexity measure is introduced for finite‐block polynomial functionals, and corresponding approximation rates are established. Numerical experiments successfully validate the theoretical predictions.
ABSTRACT In this paper, we propose an ‐splitting iteration method for solving the all‐at‐once system of weakly nonlinear equations. The system arises from the discrete semi‐linear diffusion equation with a Caputo fractional derivative. We employ L1 formula and the standard central difference scheme to discretize time and space, respectively. We then introduce the ‐splitting iteration method to tackle the all‐at‐once system of weakly nonlinear equations. Furthermore, we establish local and global convergence theories for the proposed iteration method under certain assumptions. Numerical results demonstrate that this method is effective for solving the all‐at‐once nonlinear system.
Krylov subspace methods solve large sparse linear systems by building a sequence of polynomial approximations to from successive matrix-vector products. In finite precision, the number of numerically independent directions that can be extracted from this sequence is bounded by the intrinsic information dimension , defined as the index at which the Krylov basis matrix becomes numerically rank-deficient: , where is the accumulated finite-precision noise floor. For sequences of elliptic PDE discretizations with spectral density, the normalized Gram matrix of any Chebyshev basis block decays off-diagonal at a rate , established via two integrations by parts applied to the Chebyshev moment integral; this gives (monomial basis) or (Chebyshev basis), where is the spectral condition number. Two consequences follow directly from this bound. First, any restarted Krylov method requires restart length - to reproduce unrestarted convergence; the common choices or produce 8- overhead or practical stagnation for . Second, -step GCR with Chebyshev basis and Forward Gauss-Seidel Gram solve requires only blocks of inter-block orthogonalization history, achieving storage and global synchronizations per iteration independent of iteration count; the 20- gap between this truncation depth and the restart minimum is explained by the same decay rate.
ABSTRACT Multilinear PageRank is an extension of the famous PageRank to higher order Markov chain. The minimal polynomial extrapolation method is an efficient way to speed up the fixed‐point iteration for multilinear PageRank computation. However, the storage requirements of the minimal polynomial extrapolation method bring us an enormous burden when the size of the problem is very large. Thus, it is interesting to seek new technologies that can reduce memory consumption while achieving acceleration simultaneously. On the other hand, it is only necessary to rank some of the most important elements of multilinear PageRank, in practice, rather than all of them. In this work, we first propose a partially updating extrapolation method for speeding up multilinear PageRank. Second, to further enhance the accuracy of some less important elements in the computed solution, we present a partially updating fixed‐point iteration for multilinear PageRank, and propose an improved partially updating extrapolation method. The idea is to use some important elements (i.e., we are interested in) to refine the accuracy of the less important ones. Third, theoretical analysis is given to justify the validity of the partially updating fixed‐point iteration, and the convergence behavior of the improved partially updating extrapolation method is investigated. Comprehensive numerical experiments on some real‐world data sets show that the proposed methods can reduce memory usage and accelerate computations. In particular, they are competitive candidates for large‐scale and difficult multilinear PageRank problems.
ABSTRACT The need to minimize a linear combination of an expression that involves an ‐norm of a linear transformation of the computed solution and the ‐norm of the residual error arises in image restoration as well as in statistics. A solution method that is suited for large‐scale problems is the iterative maximization‐minimization method, which determines an approximate solution in a generalized Krylov subspace. This paper describes a preconditioner for accelerating the convergence of the iterative method. The design of the preconditioner is inspired by the iterated Tikhonov regularization method. Numerical examples in image deblurring and computerized tomography show that the preconditioner significantly reduces the required CPU time and gives computed solutions of about the same quality as the unpreconditioned iterative method.
Tensor optimization problems with rank regularization have attracted significant attentions in recent years due to their extensive applications in various fields. In this paper, we focus on the CANDECOMP/PARAFAC (CP) rank regularized tensor optimization problem and equivalently transform it into a column sparse regularized problem, that is, the problem with a regularization to describe the column sparsity of CP factor matrices. Moreover, we study the relationships between the original problem and the transformed problem in the sense of global and local minimizers. To solve the nonconvex and nonsmooth transformed problem, we design an inertial block coordinate descent (iBCD) algorithm and establish its global convergence. Finally, we apply the proposed iBCD algorithm to solve tensor CP decomposition and low-CP-rank tensor completion problems. Numerical experiments on both synthetic data and real-world images validate the promising performance of our proposed method compared with several excellent methods for solving the CP rank regularized tensor optimization problem.
Sparsity models with -norm constraint are widely used in image processing, statistical regression, wireless communications, and other fields, where numerous thresholding-based methods are employed. However, when the acquisition of linear measurement is time-consuming or the measurement matrix is too large to store in memory, many thresholding-based methods that require the whole measurement matrix become impractical. To address this issue, we propose a randomized preconditioned hard thresholding pursuit algorithm for finding a sparse solution to a linear system of equations with -norm constraint. The proposed method is also applicable to general settings. It is well-known that Gaussian random matrices satisfy the restricted isometry property (RIP) with high probability, where the RIP of the measurement matrix is crucial for a thresholding-based algorithm to recover a true sparse solution effectively. Taking Gaussian measurement matrices as an illustrative example, we quantitatively analyze how the preconditioner affects the RIP of the measurement submatrix chosen uniformly at random. The linear convergence of the proposed method is also established under some assumptions. Finally, numerical examples based on synthetic and real-world data demonstrate the efficiency of the proposed method.
The classical implicit Runge-Kuta (IRK) method is a powerful method with some desirable accuracy and stability properties but is rarely used in practical applications because of its high computational cost and difficulties in preconditioning. In this paper, we introduce a space-time coupled IRK scheme that offers high-order accuracy and stability, and also has a high degree of parallelism in both space and time when used with two-level tensor-structure-preserving overlapping Schwarz preconditioners. The convergence of the proposed method is studied numerically and we show that the convergence rate depends only mildly on the mesh size, the time step size, the number of processors, and the window size. We compare the parallel performance of the proposed method with the classical method in terms of the strong scalability, and the window-size-scaled weak scalability. The numerical results indicate that the proposed method outperforms the classical method when the number of processors is large and the space-only parallelization of the classical method is a limiting factor. Moreover, in terms of the total compute time, we show numerically that higher order space-time IRK outperforms the lower order space-time IRK when a suitable window size is chosen for solving the problem in the entire space-time domain with similar accuracy.
For a class of optimal control problems constrained with certain space-fractional diffusion equations, by making use of the right rectangular rule for the cost function and the implicit finite difference scheme with the shifted Gr & uuml;nwald formula for the constraint equation along with Lagrange multiplier approach, we obtain specially structured block two-by-two linear systems. We construct the circulant-based and -matrix-based approximate block preconditioning matrices for the coefficient matrices of the discrete linear systems and analyze spectral properties of the corresponding preconditioned matrices. Theoretical results indicate that except for a small number of outliers the eigenvalues of the preconditioned matrices are clustered around 1. Numerical experiments show that these structured preconditioners can significantly improve the convergence behavior of the Krylov subspace methods.
ABSTRACT In this article, we investigate the multilevel block ‐Toeplitz matrix with , where each level possesses a block ‐Toeplitz structure and the innermost level remains unstructured. The necessary and sufficient conditions for a matrix to be a multilevel block ‐Toeplitz matrix are derived. Next, matrix expansion techniques are utilized to facilitate the application of Vandermonde factorization in solving the structured mapping problems associated with this kind of matrices. We consider the above problems whose coefficient matrices possess the multilevel block ‐Toeplitz structure, where the solutions belong to the Jordan algebra, the Lie algebra and the automorphism group. In particular, the analytical expressions for the solutions to these problems are derived. Furthermore, a complete analytical solution is provided for the open problem proposed by Mackey et al. [SIAM J. Matrix Anal. Appl., 2008, 29(4), 1389‐1410], which is a kind of structured mapping problem with unstructured coefficient matrices. Finally, several algorithms are proposed, and corresponding numerical examples are presented to verify the rationality and feasibility of the results obtained in this study.
A quaternion conjugate gradient normal equation error (QCGNE) algorithm is presented to compute the numerical solution of quaternion equations, which are widely used in color face recognition, video compression, and color image encryption. QCGNE can obtain the quaternion solution, the least Frobenius norm solution, and solve the optimal approximation problem. Compared with the real CGNE, which requires matrix dimension expansion, the proposed structure-preserving QCGNE can directly compute in the quaternion space. Therefore, it can reduce storage overhead and computational complexity, and maintain the intrinsic structural properties. To accelerate the convergence of QCGNE, two preconditioned QCGNE (PQCGNE) are proposed. Numerical examples illustrate that QCGNE and PQCGNE have the fastest convergence rate and the shortest computational time compared to the existing three algorithms. Experimental results indicate that QCGNE is feasible and much more effective for dealing with the color image encryption problem.
ABSTRACT In this paper, we propose and analyze two block Gauss‐Seidel methods, one randomized and one deterministic, to solve large linear least squares problems. We prove the convergence of the proposed methods for full‐rank overdetermined systems and provide numerical experiments to show their effectiveness compared to existing methods. In contrast to other block Kaczmarz and block coordinate descent algorithms, the proposed methods are especially advantageous when it comes to their implementation.
We present a framework for the simultaneous solution of forward () and adjoint () linear systems by reformulating the coupled problem as a single augmented nonsymmetric saddle point system. For this formulation, we derive conditions on a weight matrix under which the augmented matrix possesses a real positive spectrum, thereby permitting the use of a conjugate gradient-like iteration in a transformed inner-product setting. We also develop a Schur complement-based preconditioning strategy implemented through incomplete QR factorization of , avoiding explicit formation of the normal equations. The numerical study comprises fourteen sparse test problems from the SuiteSparse Matrix Collection, drawn from a broad range of application domains, structural classes, and conditioning regimes. The principal experiment is a matched preconditioner comparison in which NspCG, GLSQR, and MINRES share the same incomplete QR factor, so that iteration-phase runtime differences can be assessed without confounding by different factorization costs. In that setting, NspCG attains the smallest mean iteration count and the shortest mean iteration-phase time, with a clear advantage over MINRES and an essentially tied per-iteration cost with GLSQR. A brief secondary comparison with QMR and GLSQR under incomplete LU is reported for completeness, but is interpreted as a comparison of complete solver pipelines rather than of Krylov recurrences in isolation.
In this paper, we assess the performance of adaptive and nested factorized sparse approximate inverses as smoothers in multilevel V-cycles, when smoothing is performed following the Chebyshev iteration of the fourth kind, for the efficient solution of linear systems arising from a conforming discretization of higher-order partial differential equations via the partition of unity method (PUM). To this end, we consider the (anisotropic) biharmonic and triharmonic equations in two and three dimensions and discretize these problems with a - or -regular PUM and polynomial exactness up to degree 5. We adapt existing adaptive algorithms for the construction of sparse approximate inverses to the native block structure of matrices arising in the PUM. Additionally, we also present a simplified formulation of the Chebyshev iteration of the fourth kind.
Recently, a greedy quasi-block coordinate descent (GQBCD) method was proposed for solving overdetermined linear least-squares problems as published in [Applied Mathematics Letters, 171(2025), 109675]. The GQBCD method is based on a random partition and the working block is determined by a greedy strategy at each step. In this paper, we introduce a modified version of the greedy quasi block coordinate descent method, referred to as MGQBCD, and provide a theoretical analysis of its convergence properties. The proposed approach enhances the original GQBCD framework by generating updated approximations through error minimization over an expanded subspace formed by two working blocks. Numerical experiments confirm both the feasibility and the improved performance of the MGQBCD method compared to the original GQBCD, particularly in terms of accuracy and convergence behavior.
The locally optimal block preconditioned conjugate gradient (LOBPCG) method is a popular solver for large and sparse Hermitian eigenvalue problems. However, recently proposed alternatives for its single-vector version LOPCG indicate certain problematic cases with less accurate preconditioners and clustered target eigenvalues. Therein LOPCG suffers from convergence delays with strongly oscillating reduction of residual norms. By interpreting LOPCG as a truncated generalized Davidson iteration, we propose utilizing some specific information from the convergence history to recover the major loss caused by truncation. Two new schemes are presented: LOPCGa, where the trial subspace is augmented and managed by an angle-based criterion; and TPCGa, which employs augmented two-term recurrences with a residual-based update strategy to reduce dependence on the threshold choice. Numerical experiments demonstrate that both schemes outperform LOPCG and its recent competitors regarding the number of required steps and the total computational time.
When utilizing the finite element method for simulating SAW devices, a large-scale complex linear system of equations is generated. The efficiency of device simulation primarily depends on the storage requirement and computational speed for solving the linear system. However, the convergence rate of an iterative solver is confined by the ill-conditioned nature of the coefficient matrix, which is non-Hermitian and indefinite. In this paper, we propose several approximations for the Schur complement matrix. These approximations can be generated by the finite element method instead of matrix-matrix multiplication and matrix inversion operations. Based on them, we develop the lower triangular preconditioners for the coefficient matrix and analyze the eigenvalue distribution of the preconditioned matrices as well. Theoretical analysis reveals that the preconditioned matrix is diagonalizable, with its eigenvalues confined to a small rectangular region, thus ensuring rapid convergence of iterative solvers. The theoretical considerations are supported by numerical results.
In 2021, B & ouml;ttcher, Gasca, Grudsky, and Kozak showed that the limit set of the spectra of tetradiagonal Toeplitz matrices consists of one, two, or three analytic arcs in the complex plane. Building on this result, in two recent papers, the authors constructed and rigorously justified a uniform asymptotic expansion for all eigenvalues in the case where the limit set is a single arc. In the present work, we take the next natural step in addressing the asymptotic approximation of the corresponding eigenvectors in this single-arc setting. The resulting formulas explicitly reveal the structure of the eigenvectors. We also provide numerical examples that illustrate the high accuracy and linear-time computability of the proposed approximations.
