In this paper, we consider the efficient solving of the resulting algebraic system for elliptic optimal control problems with mixed finite element discretization. We propose a block-diagonal preconditioner for the symmetric and indefinite algebraic system solved with minimum residual method, which is proved to be robust and optimal with respect to both the mesh size and the regularization parameter. The block-diagonal preconditioner is constructed based on an isomorphism between appropriately chosen solution space and its dual for a general control problem with both state and gradient state observations in the objective functional. Numerical experiments confirm the efficiency of our proposed preconditioner.

Matrix equations of the kind *A*_{1}*X*^{2}+*A*_{0}*X*+*A*_{−1}=*X*, where both the matrix coefficients and the unknown are semi-infinite matrices belonging to a Banach algebra, are considered. These equations, where coefficients are quasi-Toeplitz matrices, are encountered in certain quasi-birth–death processes as the tandem Jackson queue or in any other processes that can be modeled as a reflecting random walk in the quarter plane. We provide a numerical framework for approximating the minimal nonnegative solution of these equations that relies on semi-infinite quasi-Toeplitz matrix arithmetic. In particular, we show that the algorithm of cyclic reduction can be effectively applied and can approximate the infinite-dimensional solutions with quadratic convergence at a cost that is comparable to that of the finite case. This way, we may compute a finite approximation of the sought solution and of the invariant probability measure of the associated quasi-birth–death process, within a given accuracy. Numerical experiments, performed on a collection of benchmarks, confirm the theoretical analysis.

Finding the maximum eigenvalue of a symmetric tensor is an important topic in tensor computation and numerical multilinear algebra. In this paper, we introduce a new class of structured tensors called *W*-tensors, which not only extends the well-studied nonnegative tensors by allowing negative entries but also covers several important tensors arising naturally from spectral hypergraph theory. We then show that finding the maximum *H*-eigenvalue of an even-order symmetric *W*-tensor is equivalent to solving a structured semidefinite program and hence can be validated in polynomial time. This yields a highly efficient semidefinite program algorithm for computing the maximum *H*-eigenvalue of *W*-tensors and is based on a new structured sums-of-squares decomposition result for a nonnegative polynomial induced by *W*-tensors. Numerical experiments illustrate that the proposed algorithm can successfully find the maximum *H*-eigenvalue of W-tensors with dimension up to 10,000, subject to machine precision. As applications, we provide a polynomial time algorithm for computing the maximum *H*-eigenvalues of large-size Laplacian tensors of hyperstars and hypertrees, where the algorithm can be up to 13 times faster than the state-of-the-art numerical method introduced by Ng, Qi, and Zhou in 2009. Finally, we also show that the proposed algorithm can be used to test the copositivity of a multivariate form associated with symmetric extended *Z*-tensors, whose order may be even or odd.

This paper is devoted to the study of two high-order families of frozen Newton-type methods. The methods are free of bilinear operators, which constitute the main limitation of the classical high-order iterative schemes. Both families are natural generalizations of an efficient third-order method. Although the methods are more demanding, a semilocal convergence analysis is presented using weaker conditions.

We consider subspace iteration (or projection-based) algorithms for computing those eigenvalues (and associated eigenvectors) of a Hermitian matrix that lie in a prescribed interval. For the case that the projector is approximated with polynomials, we present an adaptive strategy for selecting the degree of these polynomials such that convergence is achieved with near-to-optimum overall work without detailed a priori knowledge about the eigenvalue distribution. The idea is then transferred to the approximation of the projector by numerical integration, which corresponds to FEAST algorithm proposed by E. Polizzi in 2009. [E. Polizzi: Density-matrix-based algorithm for solving eigenvalue problems. *Phys. Rev. B* 2009; **79**:115112]. Here, our adaptation controls the number of integration nodes. We also discuss the interaction of the method with search space reduction methods.

We compare some alternative methods for computing solutions of underdetermined linear systems, *A**x*=*b*. Each method involves solving an associated system with a different nonsingular coefficient matrix,
. We obtain bounds on the condition numbers of these nonsingular matrices and test the methods on numerical examples. We discuss implications for computing eigenvector derivatives and make some recommendations.

In this paper, we consider an efficient approach to solve a linear ill-posed inverse problem with total variation regularization, where has a Kronecker product structure, . A numerical scheme is developed by using a variable splitting technique and directly exploiting the Kronecker structure of the problem. Experimental results on image restoration applications illustrate the effectiveness of our proposed method.

For large sparse non-Hermitian positive definite linear systems, we establish exact and inexact quasi-HSS iteration methods and discuss their convergence properties. Numerical experiments show that both iteration methods are effective and robust when they are used either as linear solvers or as matrix splitting preconditioners for the Krylov subspace iteration methods. In addition, these two iteration methods are, respectively, much more powerful than the exact and inexact HSS iteration methods, especially when the linear systems have nearly singular Hermitian parts or strongly dominant skew-Hermitian parts, and they can be employed to solve non-Hermitian indefinite linear systems with only mild indefiniteness.

We suggest and compare different methods for the numerical solution of Lyapunov-like equations with application to control of Markovian jump linear systems. First, we consider fixed-point iterations and associated Krylov subspace formulations. Second, we reformulate the equation as an optimization problem and consider steepest descent, conjugate gradient, and trust-region methods. Numerical experiments illustrate that, for large-scale problems, the trust-region method is more efficient than a direct solution or a standard solver for linear matrix inequalities. The fixed-point approach, however, is superior to the optimization methods. As an application, we consider a networked control system, where the Markov jumps are induced by the wireless communication protocol.

In this paper, we introduce and analyze a new singular value decomposition (SVD) called weighted SVD (WSVD) using a new inner product instead of the Euclidean one. We use the WSVD to approximate the singular values and the singular functions of the Fredholm integral operators. In this case, the new inner product arises from the numerical integration used to discretize the operator. Then, the truncated WSVD (TWSVD) is used to regularize the Nyström discretization of the first-kind Fredholm integral equations. Also, we consider the weighted LSQR (WLSQR) to approximate the solution obtained by the TWSVD method for large problems. Numerical experiments on a few problems are used to illustrate that the TWSVD can perform better than the TSVD.

Recently, there has been much interest in the solution of differential equations on surfaces and manifolds, driven by many applications whose dynamics take place on such domains. Although increasingly powerful algorithms have been developed in this field, many straightforward questions remain, particularly in the area of coupling advanced discretizations with efficient linear solvers. In this paper, we develop a structured refinement algorithm for octahedral triangulations of the surface of the sphere. We explain the composite-grid finite-element discretization of the Laplace–Beltrami operator on such triangulations and extend the fast adaptive composite-grid scheme to provide an efficient solution of the resulting linear system. Supporting numerical examples are presented, including the recovery of second-order accuracy in the case of a nonsmooth solution.

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.

No abstract is available for this article.

]]>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.

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.

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.

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.

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.

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.

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.

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.

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.

Standard numerical algorithms, such as the fast multipole method or
-matrix schemes, rely on low-rank approximations of the underlying kernel function. For high-frequency problems, the ranks grow rapidly as the mesh is refined, and standard techniques are no longer attractive. *Directional* compression techniques solve this problem by using decompositions based on plane waves. Taking advantage of hierarchical relations between these waves' directions, an efficient approximation is obtained. This paper is dedicated to *directional*
-*matrices* that employ local low-rank approximations to handle directional representations efficiently. The key result is an algorithm that takes an arbitrary matrix and finds a quasi-optimal approximation of this matrix as a directional
-matrix using a prescribed block tree. The algorithm can reach any given accuracy, and the approximation requires only
units of storage, where *n* is the matrix dimension, *κ* is the wave number, and *k* is the local rank. In particular, we have a complexity of
if *κ* is constant and
for high-frequency problems characterized by *κ*^{2}∼*n*. Because the algorithm can be applied to arbitrary matrices, it can serve as the foundation of fast techniques for constructing preconditioners.