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.

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.

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.

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.

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.

No abstract is available for this article.

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

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.

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.

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.

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.