Two classes of methods for approximate matrix inversion with convergence orders *p*=3∗2^{k}+1 (Class 1) and *p*=5∗2^{k}−1 (Class 2), *k*≥1 an integer, are given based on matrix multiplication and matrix addition. These methods perform less number of matrix multiplications compared to the known hyperpower method or *p*th-order method for the same orders and can be used to construct approximate inverse preconditioners for solving linear systems. Convergence, error, and stability analyses of the proposed classes of methods are provided. Theoretical results are justified with numerical results obtained by using the proposed methods of orders *p*=7,13 from Class 1 and the methods with orders *p*=9,19 from Class 2 to obtain polynomial preconditioners for preconditioning the biconjugate gradient (BICG) method for solving well- and ill-posed problems. From the literature, methods with orders *p*=8,16 belonging to a family developed by the effective representation of the *p*th-order method for orders *p*=2^{k}, *k* is integer *k*≥1, and other recently given high-order convergent methods of orders *p*=6,7,8,12 for approximate matrix inversion are also used to construct polynomial preconditioners for preconditioning the BICG method to solve the considered problems. Numerical comparisons are given to show the applicability, stability, and computational complexity of the proposed methods by paying attention to the asymptotic convergence rates. It is shown that the BICG method converges very quickly when applied to solve the preconditioned system. Therefore, the cost of constructing these preconditioners is amortized if the preconditioner is to be reused over several systems of same coefficient matrix with different right sides.

Algebraic multigrid (AMG) preconditioners are considered for discretized systems of partial differential equations (PDEs) where unknowns associated with different physical quantities are not necessarily colocated at mesh points. Specifically, we investigate a **Q**_{2}−**Q**_{1} mixed finite element discretization of the incompressible Navier–Stokes equations where the number of velocity nodes is much greater than the number of pressure nodes. Consequently, some velocity degrees of freedom (DOFs) are defined at spatial locations where there are no corresponding pressure DOFs. Thus, AMG approaches leveraging this colocated structure are not applicable. This paper instead proposes an automatic AMG coarsening that mimics certain pressure/velocity DOF relationships of the **Q**_{2}−**Q**_{1} discretization. The main idea is to first automatically define coarse pressures in a somewhat standard AMG fashion and then to carefully (but automatically) choose coarse velocity unknowns so that the spatial location relationship between pressure and velocity DOFs resembles that on the finest grid. To define coefficients within the intergrid transfers, an energy minimization AMG (EMIN-AMG) is utilized. EMIN-AMG is not tied to specific coarsening schemes and grid transfer sparsity patterns, and so it is applicable to the proposed coarsening. Numerical results highlighting solver performance are given on Stokes and incompressible Navier–Stokes problems.

We discuss the convergence of a two-level version of the multilevel Krylov method for solving linear systems of equations with symmetric positive semidefinite matrix of coefficients. The analysis is based on the convergence result of Brown and Walker for the Generalized Minimal Residual method (GMRES), with the left- and right-preconditioning implementation of the method. Numerical results based on diffusion problems are presented to show the convergence.

The uniqueness of multilinear PageRank vectors is discussed, and the new uniqueness condition is given. The new results are better than the one given in the work of Gleich et al. published in SIAM J Matrix Anal Appl. 2015;36;1409-1465. Numerical examples are given to demonstrate the new theoretical results.

For the numerical solution of time-dependent partial differential equations, time-parallel methods have recently been shown to provide a promising way to extend prevailing strong-scaling limits of numerical codes. One of the most complex methods in this field is the “Parallel Full Approximation Scheme in Space and Time” (PFASST). PFASST already shows promising results for many use cases and benchmarks. However, a solid and reliable mathematical foundation is still missing. We show that, under certain assumptions, the PFASST algorithm can be conveniently and rigorously described as a multigrid-in-time method. Following this equivalence, first steps towards a comprehensive analysis of PFASST using blockwise local Fourier analysis are taken. The theoretical results are applied to examples of diffusive and advective type.

Null-space methods for solving saddle point systems of equations have long been used to transform an indefinite system into a symmetric positive definite one of smaller dimension. A number of independent works in the literature have identified that we can interpret a null-space method as a matrix factorization. We review these findings, highlight links between them, and bring them into a unified framework. We also investigate the suitability of using null-space factorizations to derive sparse direct methods and present numerical results for both practical and academic problems.

We consider the numerical solution of a c-stable linear equation in the tensor product space
, arising from a discretized elliptic partial differential equation in
. Utilizing the stability, we produce an equivalent d-stable generalized Stein-like equation, which can be solved iteratively. For large-scale problems defined by sparse and structured matrices, the methods can be modified for further efficiency, producing algorithms of
computational complexity, under appropriate assumptions (with *n*_{s} being the flop count for solving a linear system associated with
). Illustrative numerical examples will be presented.

By employing modulus-based matrix splitting iteration methods as smoothers, we establish modulus-based multigrid methods for solving large sparse linear complementarity problems. The local Fourier analysis is used to quantitatively predict the asymptotic convergence factor of this class of multigrid methods. Numerical results indicate that the modulus-based multigrid methods of the W-cycle can achieve optimality in terms of both convergence factor and computing time, and their asymptotic convergence factors can be predicted perfectly by the local Fourier analysis of the corresponding modulus-based two-grid methods.

The problem of determining matrix inertia is very important in many applications, for example, in stability analysis of dynamical systems. In the (point-wise) H-matrix case, it was proven that the diagonal entries solely reveal this information. This paper generalises these results to the block H-matrix cases for 1, 2, and matrix norms. The usefulness of the block approach is illustrated on 3 relevant numerical examples, arising in engineering.

The Jacobi, Gauss-Seidel and successive over-relaxation methods are well-known basic iterative methods for solving system of linear equations. In this paper, we extend those basic methods to solve the tensor equation
, where
is an *m*th-order *n*−dimensional symmetric tensor and *b* is an *n*-dimensional vector. Under appropriate conditions, we show that the proposed methods are globally convergent and locally r-linearly convergent. Taking into account the special structure of the Newton method for the problem, we propose a Newton-Gauss-Seidel method, which is expected to converge faster than the above methods. The proposed methods can be extended to solve a general symmetric tensor equations. Our preliminary numerical results show the effectiveness of the proposed methods.

It is shown that the matrix sequence generated by Euler's method starting from the identity matrix converges to the principal *p*th root of a square matrix, if all the eigenvalues of the matrix are in a region including the one for Newton's method given by Guo in 2010. The convergence is cubic if the matrix is invertible. A modification version of Euler's method using the Schur decomposition is developed. Numerical experiments show that the modified algorithm has the overall good numerical behavior.

It is well known that the standard algorithm for the mixed least squares–total least squares (MTLS) problem uses the QR factorization to reduce the original problem into a standard total least squares problem with smaller size, which can be solved based on the singular value decomposition (SVD). In this paper, the MTLS problem is proven to be closely related to a weighted total least squares problem with its error-free columns multiplied by a large weighting factor. A criterion for choosing the weighting factor is given; and for the sake of stability in solving the MTLS problem, the Cholesky factorization-based inverse (Cho-INV) iteration and Rayleigh quotient iteration are also considered. For large-scale MTLS problems, numerical tests show that Cho-INV is superior to the standard QR-SVD method, especially for the case with big gap between the desired and undesired singular values and the case when the coefficient matrix has much more error-contaminated columns. Rayleigh quotient iteration behaves more efficient than QR-SVD for most cases and fails occasionally, and in some cases, it converges much faster than Cho-INV but still less efficient due to its higher computation cost.

Image registration is a central problem in a variety of areas involving imaging techniques and is known to be challenging and ill-posed. Regularization functionals based on hyperelasticity provide a powerful mechanism for limiting the ill-posedness. A key feature of hyperelastic image registration approaches is their ability to model large deformations while guaranteeing their invertibility, which is crucial in many applications. To ensure that numerical solutions satisfy this requirement, we discretize the variational problem using piecewise linear finite elements, and then solve the discrete optimization problem using the Gauss–Newton method. In this work, we focus on computational challenges arising in approximately solving the Hessian system. We show that the Hessian is a discretization of a strongly coupled system of partial differential equations whose coefficients can be severely inhomogeneous. Motivated by a local Fourier analysis, we stabilize the system by thresholding the coefficients. We propose a Galerkin-multigrid scheme with a collective pointwise smoother. We demonstrate the accuracy and effectiveness of the proposed scheme, first on a two-dimensional problem of a moderate size and then on a large-scale real-world application with almost 9 million degrees of freedom.

A fast algorithm for enclosing the solution of the nonsymmetric algebraic Riccati equation arising in transport theory is proposed. The equation has a special structure, which is taken into account to reduce the complexity. By exploiting the structure, the enclosing process involves only quadratic complexity under a reasonable assumption. The algorithm moreover verifies the uniqueness and minimal positiveness of the enclosed solution. Numerical results show the efficiency of the algorithm.

Sparse symmetric indefinite linear systems of equations arise in numerous practical applications. In many situations, an iterative method is the method of choice but a preconditioner is normally required for it to be effective. In this paper, the focus is on a class of incomplete factorization algorithms that can be used to compute preconditioners for symmetric indefinite systems. A limited memory approach is employed that incorporates a number of new ideas with the goal of improving the stability, robustness, and efficiency of the preconditioner. These include the monitoring of stability as the factorization proceeds and the incorporation of pivot modifications when potential instability is observed. Numerical experiments involving test problems arising from a range of real-world applications demonstrate the effectiveness of our approach.

In this paper, we study a deblurring algorithm for distorted images by random impulse response. We propose and develop a convex optimization model to recover the underlying image and the blurring function simultaneously. The objective function is composed of 3 terms: the data-fitting term between the observed image and the product of the estimated blurring function and the estimated image, the squared difference between the estimated blurring function and its mean, and the total variation regularization term for the estimated image. We theoretically show that under some mild conditions, the resulting objective function can be convex in which the global minimum value is unique. The numerical results confirm that the peak-to-signal-noise-ratio and structural similarity of the restored images by the proposed algorithm are the best when the proposed objective function is convex. We also present a proximal alternating minimization scheme to solve the resulting minimization problem. Numerical examples are presented to demonstrate the effectiveness of the proposed model and the efficiency of the numerical scheme.

We generalize the matrix Kronecker product to tensors and propose the tensor Kronecker product singular value decomposition that decomposes a real *k*-way tensor
into a linear combination of tensor Kronecker products with an arbitrary number of *d* factors. We show how to construct
, where each factor
is also a *k*-way tensor, thus including matrices (*k*=2) as a special case. This problem is readily solved by reshaping and permuting
into a *d*-way tensor, followed by a orthogonal polyadic decomposition. Moreover, we introduce the new notion of general symmetric tensors (encompassing symmetric, persymmetric, centrosymmetric, Toeplitz and Hankel tensors, etc.) and prove that when
is structured then its factors
will also inherit this structure.

In this paper, we compare two block triangular preconditioners for different linearizations of the Rayleigh–Bénard convection problem discretized with finite element methods. The two preconditioners differ in the nested or nonnested use of a certain approximation of the Schur complement associated to the Navier–Stokes block. First, bounds on the generalized eigenvalues are obtained for the preconditioned systems linearized with both Picard and Newton methods. Then, the performance of the proposed preconditioners is studied in terms of computational time. This investigation reveals some inconsistencies in the literature that are hereby discussed. We observe that the nonnested preconditioner works best both for the Picard and for the Newton cases. Therefore, we further investigate its performance by extending its application to a mixed Picard–Newton scheme. Numerical results of two- and three-dimensional cases show that the convergence is robust with respect to the mesh size. We also give a characterization of the performance of the various preconditioned linearization schemes in terms of the Rayleigh number.

No abstract is available for this article.

]]>Tikhonov regularization is one of the most popular approaches to solving linear discrete ill-posed problems. The choice of the regularization matrix may significantly affect the quality of the computed solution. When the regularization matrix is the identity, iterated Tikhonov regularization can yield computed approximate solutions of higher quality than (standard) Tikhonov regularization. This paper provides an analysis of iterated Tikhonov regularization with a regularization matrix different from the identity. Computed examples illustrate the performance of this method.

The treatment of the stochastic linear quadratic optimal control problem with finite time horizon requires the solution of stochastic differential Riccati equations. We propose efficient numerical methods, which exploit the particular structure and can be applied for large-scale systems. They are based on numerical methods for ordinary differential equations such as Rosenbrock methods, backward differentiation formulas, and splitting methods. The performance of our approach is tested in numerical experiments.

The contour integral-based eigensolvers are the recent efforts for computing the eigenvalues inside a given region in the complex plane. The best-known members are the Sakurai–Sugiura method, its stable version CIRR, and the FEAST algorithm. An attractive computational advantage of these methods is that they are easily parallelizable. The FEAST algorithm was developed for the generalized Hermitian eigenvalue problems. It is stable and accurate. However, it may fail when applied to non-Hermitian problems. Recently, a dual subspace FEAST algorithm was proposed to extend the FEAST algorithm to non-Hermitian problems. In this paper, we instead use the oblique projection technique to extend FEAST to the non-Hermitian problems. Our approach can be summarized as follows: (a) construct a particular contour integral to form a search subspace containing the desired eigenspace and (b) use the oblique projection technique to extract desired eigenpairs with appropriately chosen test subspace. The related mathematical framework is established. Comparing to the dual subspace FEAST algorithm, we can save the computational cost roughly by a half if only the eigenvalues or the eigenvalues together with their right eigenvectors are needed. We also address some implementation issues such as how to choose a suitable starting matrix and design-efficient stopping criteria. Numerical experiments are provided to illustrate that our method is stable and efficient.

The finite difference discretization of the spatial fractional diffusion equations gives discretized linear systems whose coefficient matrices have a diagonal-plus-Toeplitz structure. For solving these diagonal-plus-Toeplitz linear systems, we construct a class of diagonal and Toeplitz splitting iteration methods and establish its unconditional convergence theory. In particular, we derive a sharp upper bound about its asymptotic convergence rate and deduct the optimal value of its iteration parameter. The diagonal and Toeplitz splitting iteration method naturally leads to a diagonal and circulant splitting preconditioner. Analysis shows that the eigenvalues of the corresponding preconditioned matrix are clustered around 1, especially when the discretization step-size *h* is small. Numerical results exhibit that the diagonal and circulant splitting preconditioner can significantly improve the convergence properties of GMRES and BiCGSTAB, and these preconditioned Krylov subspace iteration methods outperform the conjugate gradient method preconditioned by the approximate inverse circulant-plus-diagonal preconditioner proposed recently by Ng and Pan (M.K. Ng and J.-Y. Pan, SIAM J. Sci. Comput. 2010;32:1442-1464). Moreover, unlike this preconditioned conjugate gradient method, the preconditioned GMRES and BiCGSTAB methods show *h*-independent convergence behavior even for the spatial fractional diffusion equations of discontinuous or big-jump coefficients.