Usefulness of spectral localizations in analysis of various matrix properties, such as stability of dynamical systems, has led us to derive a pseudospectra localization technique using the ideas that come from diagonally dominant matrices. In such way, many theoretical and practical applications of pseudospectra (robust stability, transient behavior, nonnormal dynamics, etc.) can be linked with specific relationships between matrix entries. This allows one to understand certain phenomena that occur in practice better, as we show for the realistic model of soil energetic food web. The novelty of the presented results, therefore, lies not only in new mathematical formulations but also in the conceptual sense because it links stability with empirical data and their uncertainty limitations. Copyright © 2015 John Wiley & Sons, Ltd.

We describe randomized algorithms for computing the dominant eigenmodes of the generalized Hermitian eigenvalue problem *A**x* = *λ**B**x*, with *A* Hermitian and *B* Hermitian and positive definite. The algorithms we describe only require forming operations *A**x*,*B**x* and *B*^{−1}*x* and avoid forming square roots of *B* (or operations of the form, *B*^{1/2}*x* or *B*^{−1/2}*x*). We provide a convergence analysis and a posteriori error bounds and derive some new results that provide insight into the accuracy of the eigenvalue calculations. The error analysis shows that the randomized algorithm is most accurate when the generalized singular values of *B*^{−1}*A* decay rapidly. A randomized algorithm for the generalized singular value decomposition is also provided. Finally, we demonstrate the performance of our algorithm on computing an approximation to the Karhunen–Loève expansion, which involves a computationally intensive generalized Hermitian eigenvalue problem with rapidly decaying eigenvalues. Copyright © 2015 John Wiley & Sons, Ltd.

The low rank solution of the Q-weighted nearest correlation matrix problem is studied in this paper. Based on the property of Q-weighted norm and the Gramian representation, we first reformulate the Q-weighted nearest correlation matrix problem as a minimization problem of the trace function with quadratic constraints and then prove that the solution of the minimization problem the trace function is the stationary point of the original problem if it is rank-deficient. Finally, the nonmonotone spectral projected gradient method is constructed to solve them. Numerical examples illustrate that the new method is feasible and effective. Copyright © 2015 John Wiley & Sons, Ltd.

A new fast algebraic method for obtaining an -approximation of a matrix from its entries is presented. The main idea behind the method is based on the nested representation and the maximum volume principle to select submatrices in low-rank matrices. A special iterative approach for the computation of so-called representing sets is established. The main advantage of the method is that it uses only the hierarchical partitioning of the matrix and does not require special ‘proxy surfaces’ to be selected in advance.

The numerical experiments for the electrostatic problem and for the boundary integral operator confirm the effectiveness and robustness of the approach. The complexity is linear in the matrix size and polynomial in the ranks. The algorithm is implemented as an open-source Python package that is available online. Copyright © 2015 John Wiley & Sons, Ltd.

Among numerous iterative methods for solving the minimal nonnegative solution of an *M*-matrix algebraic Riccati equation, the structure-preserving doubling algorithm (SDA) stands out owing to its overall efficiency as well as accuracy. SDA is globally convergent and its convergence is quadratic, except for the critical case for which it converges linearly with the linear rate 1/2. In this paper, we first undertake a delineatory convergence analysis that reveals that the approximations by SDA can be decomposed into two components: the stable component that converges quadratically and the rank-one component that converges linearly with the linear rate 1/2. Our analysis also shows that as soon as the stable component is fully converged, the rank-one component can be accurately recovered. We then propose an efficient hybrid method, called the *two-phase SDA*, for which the SDA iteration is stopped as soon as it is determined that the stable component is fully converged. Therefore, this two-phase SDA saves those SDA iterative steps that previously have to have for the rank-one component to be computed accurately, and thus essentially, it can be regarded as a quadratically convergent method. Numerical results confirm our analysis and demonstrate the efficiency of the new two-phase SDA. Copyright © 2015 John Wiley & Sons, Ltd.

The existing algorithms for computing the minimal Geršgorin set are designed for small and medium size (irreducible) matrices and based on Perron root computations coupled with bisection method and sampling techniques. Here, we first discuss the drawbacks of the existing methods and present a new approach based on the modified Newton's method to find zeros of the parameter dependent left-most eigenvalue of a *Z*-matrix and a special curve tracing procedure. The advantages of the new approach are presented on several test examples that arise in practical applications. Copyright © 2015 John Wiley & Sons, Ltd.

In this paper, we present a method for fast summation of long-range potentials on 3D lattices with multiple defects and having non-rectangular geometries, based on rank-structured tensor representations. This is a significant generalization of our recent technique for the grid-based summation of electrostatic potentials on the rectangular *L* × *L* × *L* lattices by using the canonical tensor decompositions and yielding the *O*(*L*) computational complexity instead of *O*(*L*^{3}) by traditional approaches. The resulting lattice sum is calculated as a Tucker or canonical representation whose directional vectors are assembled by the 1D summation of the generating vectors for the shifted reference tensor, once precomputed on large *N* × *N* × *N* representation grid in a 3D bounding box. The tensor numerical treatment of defects is performed in an algebraic way by simple summation of tensors in the canonical or Tucker formats. To diminish the considerable increase in the tensor rank of the resulting potential sum, the *ϵ*-rank reduction procedure is applied based on the generalized reduced higher-order SVD scheme. For the reduced higher-order SVD approximation to a sum of canonical/Tucker tensors, we prove the stable error bounds in the relative norm in terms of discarded singular values of the side matrices. The required storage scales linearly in the 1D grid-size, *O*(*N*), while the numerical cost is estimated by *O*(*N**L*). The approach applies to a general class of kernel functions including those for the Newton, Slater, Yukawa, Lennard-Jones, and dipole-dipole interactions. Numerical tests confirm the efficiency of the presented tensor summation method; we demonstrate that a sum of millions of Newton kernels on a 3D lattice with defects/impurities can be computed in seconds in Matlab implementation. The tensor approach is advantageous in further functional calculus with the lattice potential sums represented on a 3D grid, like integration or differentiation, using tensor arithmetics of 1D complexity. Copyright © 2015 John Wiley & Sons, Ltd.

The geometric mean of two matrices is considered from a computational viewpoint. Several numerical algorithms based on different properties and representations of the geometric mean are discussed and analyzed. It is shown that most of the algorithms can be classified in terms of the rational approximations of the inverse square root function. A review of relevant applications is given. Copyright © 2015 John Wiley & Sons, Ltd.

The numerical computation of Lagrangian invariant subspaces of large-scale Hamiltonian matrices is discussed in the context of the solution of Lyapunov equations. A new version of the low-rank alternating direction implicit method is introduced, which, in order to avoid numerical difficulties with solutions that are of very large norm, uses an inverse-free representation of the subspace and avoids inverses of ill-conditioned matrices. It is shown that this prevents large growth of the elements of the solution that may destroy a low-rank approximation of the solution. A partial error analysis is presented, and the behavior of the method is demonstrated via several numerical examples. Copyright © 2015 John Wiley & Sons, Ltd.

We introduce a numerical method for the numerical solution of the *Lur'e equations*, a system of matrix equations that arises, for instance, in linear-quadratic infinite time horizon optimal control. We focus on small-scale, dense problems. Via a Cayley transformation, the problem is transformed to the discrete-time case, and the structural infinite eigenvalues of the associated matrix pencil are deflated. The deflated problem is associated with a symplectic pencil with several Jordan blocks of eigenvalue 1 and even size, which arise from the nontrivial Kronecker chains at infinity of the original problem. For the solution of this modified problem, we use the *structure-preserving doubling algorithm*. Implementation issues such as the choice of the parameter *γ* in the Cayley transform are discussed. The most interesting feature of this method, with respect to the competing approaches, is the absence of arbitrary rank decisions, which may be ill-posed and numerically troublesome. The numerical examples presented confirm the effectiveness of this method. Copyright © 2015 John Wiley & Sons, Ltd.

The symmetric Lanczos method is commonly applied to reduce large-scale symmetric linear discrete ill-posed problems to small ones with a symmetric tridiagonal matrix. We investigate how quickly the nonnegative subdiagonal entries of this matrix decay to zero. Their fast decay to zero suggests that there is little benefit in expressing the solution of the discrete ill-posed problems in terms of the eigenvectors of the matrix compared with using a basis of Lanczos vectors, which are cheaper to compute. Similarly, we show that the solution subspace determined by the LSQR method when applied to the solution of linear discrete ill-posed problems with a nonsymmetric matrix often can be used instead of the solution subspace determined by the singular value decomposition without significant, if any, reduction of the quality of the computed solution. Copyright © 2015 John Wiley & Sons, Ltd.

Large-scale reservoir simulations are extremely time-consuming because of the solution of large-scale linear systems arising from the Newton or Newton–Raphson iterations. The problem becomes even worse when highly heterogeneous geological models are employed. This paper introduces a family of multi-stage preconditioners for parallel black oil simulations, which are based on the famous constrained pressure residual preconditioner. Numerical experiments demonstrate that our preconditioners are robust, efficient, and scalable. Copyright © 2015 John Wiley & Sons, Ltd.

We perform a spectral analysis of the preconditioned Hermitian/skew-Hermitian splitting (PHSS) method applied to multilevel block Toeplitz linear systems in which the coefficient matrix *T*_{n}(*f*) is associated with a Lebesgue integrable matrix-valued function *f*. When the preconditioner is chosen as a Hermitian positive definite multilevel block Toeplitz matrix *T*_{n}(*g*), the resulting sequence of PHSS iteration matrices *M*_{n} belongs to the generalized locally Toeplitz class. In this case, we are able to compute the symbol *ϕ*(*f*,*g*) describing the asymptotic eigenvalue distribution of *M*_{n}when *n**∞* and the matrix size diverges. By minimizing the infinity norm of the spectral radius of the symbol *ϕ*(*f*,*g*), we are also able to identify effective PHSS preconditioners *T*_{n}(*g*) for the matrix *T*_{n}(*f*). A number of numerical experiments are presented and commented, showing that the theoretical results are confirmed and that the spectral analysis leads to efficient PHSS methods. Copyright © 2015 John Wiley & Sons, Ltd.

We present an algebraic structured preconditioner for the iterative solution of large sparse linear systems. The preconditioner is based on a multifrontal variant of sparse LU factorization used with nested dissection ordering. Multifrontal factorization amounts to a partial factorization of a sequence of logically dense frontal matrices, and the preconditioner is obtained if structured factorization is used instead. This latter exploits the presence of low numerical rank in some off-diagonal blocks of the frontal matrices. An algebraic procedure is presented that allows to identify the hierarchy of the off-diagonal blocks with low numerical rank based on the sparsity of the system matrix. This procedure is motivated by a model problem analysis, yet numerical experiments show that it is successful beyond the model problem scope. Further aspects relevant for the algebraic structured preconditioner are discussed and illustrated with numerical experiments. The preconditioner is also compared with other solvers, including the corresponding direct solver. Copyright © 2015 John Wiley & Sons, Ltd.

We construct, analyze, and implement SSOR-like preconditioners for non-Hermitian positive definite system of linear equations when its coefficient matrix possesses either a dominant Hermitian part or a dominant skew-Hermitian part. We derive tight bounds for eigenvalues of the preconditioned matrices and obtain convergence rates of the corresponding SSOR-like iteration methods as well as the corresponding preconditioned GMRES iteration methods. Numerical implementations show that Krylov subspace iteration methods such as GMRES, when accelerated by the SSOR-like preconditioners, are efficient solvers for these classes of non-Hermitian positive definite linear systems. Copyright © 2015 John Wiley & Sons, Ltd.

This paper presents parallel preconditioners and multigrid solvers for solving linear systems of equations arising from stochastic polynomial chaos formulations of the diffusion equation with random coefficients. These preconditioners and solvers are extensions of the preconditioner developed in an earlier paper for strongly coupled systems of elliptic partial differential equations that are norm equivalent to systems that can be factored into an algebraic coupling component and a diagonal differential component. The first preconditioner, which is applied to the norm equivalent system, is obtained by sparsifying the inverse of the algebraic coupling component. This sparsification leads to an efficient method for solving these systems at the large scale, even for problems with large statistical variations in the random coefficients. An extension of this preconditioner leads to stand-alone multigrid methods that can be applied directly to the actual system rather than to the norm equivalent system. These multigrid methods exploit the algebraic/differential factorization of the norm equivalent systems to produce variable-decoupled systems on the coarse levels. Moreover, the structure of these methods allows easy software implementation through re-use of robust high-performance software such as the Hypre library package. Two-grid matrix bounds will be established, and numerical results will be given. Copyright © 2015 John Wiley & Sons, Ltd.

The following paper discusses the application of a multigrid-in-time scheme to Least Squares Shadowing (LSS), a novel sensitivity analysis method for chaotic dynamical systems. While traditional sensitivity analysis methods break down for chaotic dynamical systems, LSS is able to compute accurate gradients. Multigrid is used because LSS requires solving a very large Karush–Kuhn–Tucker system constructed from the solution of the dynamical system over the entire time interval of interest. Several different multigrid-in-time schemes are examined, and a number of factors were found to heavily influence the convergence rate of multigrid-in-time for LSS. These include the iterative method used for the smoother, how the coarse grid system is formed and how the least squares objective function at the center of LSS is weighted. Copyright © 2014 John Wiley & Sons, Ltd.

No abstract is available for this article.

]]>We extend the balancing domain decomposition by constraints (BDDC) method to flows in porous media discretised by mixed-hybrid finite elements with combined mesh dimensions. Such discretisations appear when major geological fractures are modelled by one-dimensional or two-dimensional elements inside three-dimensional domains. In this set-up, the global problem and the substructure problems have a symmetric saddle-point structure, containing a ‘penalty’ block due to the combination of meshes. We show that the problem can be reduced by means of iterative substructuring to an interface problem, which is symmetric and positive definite. The interface problem can thus be solved by conjugate gradients with the BDDC method as a preconditioner. A parallel implementation of this algorithm is incorporated into an existing software package for subsurface flow simulations. We study the performance of the iterative solver on several academic and real-world problems. Numerical experiments illustrate its efficiency and scalability. Copyright © 2015 John Wiley & Sons, Ltd.

High-order strongly stable time integration method enables use of large time steps and is applicable also for differential–algebraic problems, without any order reduction. At each time step, frequently a large-scale linear algebraic system must be solved. To solve the arising block matrix systems, an efficient preconditioning method is presented and analysed. The method is applied for the solution of a consolidation problem arising in poroelasticity. Copyright © 2015 John Wiley & Sons, Ltd.

Our method is a kind of exact interpolation scheme by Achi Brandt *et al*. In the exact interpolation scheme, for **x** being the fine-level approximation of the solution, the coarse-space *V* = *V*(**x**) is constructed so that **x** satisfies **x**∈*V*. We achieve it simply by adding the vector **x** as a first column of the prolongator. (The columns of the prolongator *P* form a computationally relevant basis of the coarse-space *V* = Range(*P*).) The advantages of this construction become obvious when solving non-linear problems. The cost of enriching the coarse-space *V* by the current approximation **x** is a single dense column of the prolongator that has to be updated each iteration. Our method can be used for multilevel acceleration of virtually any iterative method used for solving both linear and non-linear systems. Copyright © 2015 John Wiley & Sons, Ltd.

In this paper, the idea of auxiliary space multigrid methods is introduced. The construction is based on a two-level block factorization of local (finite element stiffness) matrices associated with a partitioning of the domain into overlapping or non-overlapping subdomains. The two-level method utilizes a coarse-grid operator obtained from additive Schur complement approximation. Its analysis is carried out in the framework of auxiliary space preconditioning and condition number estimates for both the two-level preconditioner and the additive Schur complement approximation are derived. The two-level method is recursively extended to define the auxiliary space multigrid algorithm. In particular, so-called Krylov cycles are considered. The theoretical results are supported by a representative collection of numerical tests that further demonstrate the efficiency of the new algorithm for multiscale problems. Copyright © 2014 John Wiley & Sons, Ltd.

A cheap symmetric stiffness-based preconditioning of the Hessian of the dual problem arising from the application of the finite element tearing and interconnecting domain decomposition to the solution of variational inequalities with varying coefficients is proposed. The preconditioning preserves the structure of the inequality constraints and affects both the linear and nonlinear steps, so that it can improve the rate of convergence of the algorithms that exploit the conjugate gradient steps or the gradient projection steps. The bounds on the regular condition number of the Hessian of the preconditioned problem, which are independent of the coefficients, are given. The related stiffness scaling is also considered and analysed. The improvement is demonstrated by numerical experiments including the solution of a contact problem with variationally consistent discretization of the non-penetration conditions. The results are relevant also for linear problems. Copyright © 2015 John Wiley & Sons, Ltd.

In this paper, we study tensors with time dimension which arises in many scientific and engineering applications such as time series gene expression analysis and video analysis. In these applications, we are interested in determining a set of components interacting closely and consistently together over a period of time. The main aim of this paper is to develop a numerical method to compute such constrained CANDECOMP/PARAFAC (CP) decompositions. We make use of the total variation regularization to constrain the time dimension factor in the decomposition in order to obtain a piecewise constant function with respect to time points. The components of the other dimensions corresponding to these time points are closely related. For example, in time series gene expression analysis, a set of genes may regulate a biological process together during a specific time period; in video analysis, a set of image pixels may refer to a foreground object in the video frames. We employ ADMM to solve the resulting optimization problem, and study its convergence. Numerical examples on synthetic and real data sets are used to demonstrate that the proposed total variation based CP decomposition model can provide more accurate and interesting results. Copyright © 2015 John Wiley & Sons, Ltd.

In this paper, we extend the inexact Uzawa algorithm of Hu and Zou to the nonsymmetric generalized saddle point problem (SPP). The techniques used here are similar with those in the paper of Bramble *et al.*, where the convergence of Uzawa type algorithm for solving nonsymmetric case depends on the spectrum of the preconditioners involved. The main contributions of this paper focus on two new linear Uzawa type algorithms for nonsymmetric generalized saddle point problems and their convergence properties. Our exact Uzawa algorithm can always converge without any prior estimate on the spectrum of two preconditioned subsystems involved. The convergence of our inexact Uzawa algorithm does not need the prior spectrum estimate of Schur complement, which may not be easy to achieve in practical applications. Numerical results of the algorithm on the Navier–Stokes problem are also presented. Copyright © 2015 John Wiley & Sons, Ltd.

In this work, we consider the numerical solution of large nonlinear eigenvalue problems that arise in thermoacoustic simulations involved in the stability analysis of large combustion devices. We briefly introduce the physical modeling that leads to a nonlinear eigenvalue problem that is solved using a nonlinear fixed point iteration scheme. Each step of this nonlinear method requires the solution of a complex non-Hermitian linear eigenvalue problem. We review a set of state-of-the-art eigensolvers and discuss strategies to recycle spectral information from one nonlinear step to the next. More precisely, we consider the Jacobi–Davidson algorithm, the Implicitly Restarted Arnoldi method, the Krylov–Schur solver and its block-variant, and the subspace iteration method with Chebyshev acceleration. On a small test example, we study the relevance of the different approaches and illustrate on a large industrial test case the performance of the parallel solvers best suited to recycle spectral information for large-scale thermoacoustic stability analysis. Copyright © 2015 John Wiley & Sons, Ltd.

In this paper, we study the convergence property of a block improvement method (BIM) for the bi-quadratic polynomial optimization problem over unit spheres. We establish the global convergence of the method generally and establish its linear convergence rate under the second-order sufficient condition. We also extend the BIM to inhomogeneous polynomial optimization problems over unit spheres. Numerical results reported in this paper show that the method is promising. Copyright © 2015 John Wiley & Sons, Ltd.

Discrete representations of the Helmholtz operator generally give rise to extremely difficult linear systems from an iterative solver perspective. This is due in part to the large oscillatory near null space of the linear system. Typical iterative methods do not effectively reduce error components in the subspace associated with this near null space. Traditional coarse grids used within multilevel solvers also cannot capture these components because of their oscillatory nature. While the shifted Laplacian is a popular method allowing multigrid preconditioners to achieve improved convergence, it is still unsatisfactory for higher-frequency problems. This paper analyzes the shifted Laplacian through polynomial smoothers to show its effect on the spectrum of the discrete Helmholtz operator. This analysis reveals desirable features of the shifted Laplacian preconditioners as well as limitations. Shifted Laplacian techniques are first extended to smoothed aggregation algebraic multigrid. Motivated by analysis, we propose augmenting the algebraic multigrid shifted Laplacian with a two-grid error correction step. This additional step consists of the generalized minimal residual iterations applied to a projected version of the unshifted equations on a single auxiliary coarse grid. Results indicate that augmentation can improve the convergence significantly and can reduce the overall solve time on two-dimensional and three-dimensional problems. Copyright © 2015 John Wiley & Sons, Ltd.

This article discusses a more general and numerically stable Rybicki Press algorithm, which enables inverting and computing determinants of covariance matrices, whose elements are sums of exponentials. The algorithm is true in exact arithmetic and relies on introducing new variables and corresponding equations, thereby converting the matrix into a banded matrix of larger size. Linear complexity banded algorithms for solving linear systems and computing determinants on the larger matrix enable linear complexity algorithms for the initial semi-separable matrix as well. Benchmarks provided illustrate the linear scaling of the algorithm. Copyright © 2015 John Wiley & Sons, Ltd.

Combining the modified matrix–vector equation approach with the technique of Lyapunov majorant function and the Banach fixed point theorem, we obtain improved rigorous perturbation bounds for the LU and QR factorizations with normwise perturbation in the given matrix. Each of the improved rigorous perturbation bounds is a rigorous version of the first-order perturbation bound derived by the matrix–vector equation approach in the literature, and we present their explicit expressions. These bounds are always tighter than those given by Chang and Stehlé in the paper entitled “Rigorous perturbation bounds of some matrix factorizations”. This fact is illustrated by numerical examples. Copyright © 2015 John Wiley & Sons, Ltd.

In various applications, for instance, in the detection of a Hopf bifurcation or in solving separable boundary value problems using the two-parameter eigenvalue problem, one has to solve a generalized eigenvalue problem with 2 × 2 operator determinants of the form

We present efficient methods that can be used to compute a small subset of the eigenvalues. For full matrices of moderate size, we propose either the standard implicitly restarted Arnoldi or Krylov–Schur iteration with shift-and-invert transformation, performed efficiently by solving a Sylvester equation. For large problems, it is more efficient to use subspace iteration based on low-rank approximations of the solution of the Sylvester equation combined with a Krylov–Schur method for the projected problems. Copyright © 2015 John Wiley & Sons, Ltd.

An efficient technique based on low-rank separated approximations is proposed for computation of three-dimensional integrals arising in the energy deposition model that describes ion-atomic collisions. Direct tensor-product quadrature requires grids of size 4000^{3} that is unacceptable. Moreover, several of such integrals have to be computed simultaneously for different values of parameters. To reduce the complexity, we use the structure of the integrand and apply numerical linear algebra techniques for the construction of low-rank approximation. The resulting algorithm is 10^{3} faster than spectral quadratures in spherical coordinates used in the original DEPOSIT code. The approach can be generalized to other multidimensional problems in physics. Copyright © 2015 John Wiley & Sons, Ltd.

Using the Schur complement of a matrix, we propose a computational framework for performing constrained maximum likelihood estimation in which the unknown parameters can be partitioned into two sets. Under appropriate regularity conditions, the corresponding estimating equations form a non-linear system of equations with constraints. Solving this system is typically accomplished via methods which require computing or estimating a Hessian matrix. We present an alternative algorithm that solves the constrained non-linear system in block coordinate descent fashion. An explicit form for the solution is given. The overall algorithm is shown in numerical studies to be faster than standard methods that either compute or approximate the Hessian as well as the classical Nelder–Mead algorithm. We apply our approach to a motivating problem of evaluating the effectiveness of Road Safety Policies. This includes several numerical studies on simulated data. Copyright © 2015 John Wiley & Sons, Ltd.

The null linear discriminant analysis method is a competitive approach for dimensionality reduction. The implementation of this method, however, is computationally expensive. Recently, a fast implementation of null linear discriminant analysis method was proposed in the paper of Sharma *et al.* In this method, a random matrix multiplication with scatter matrix is used to produce the orientation matrix. In this paper, we consider whether the random matrix can be replaced by an arbitrary full rank matrix of appropriate size, and focus on the necessary and sufficient condition to guarantee full column rank of the orientation matrix. We investigate how to choose the matrix properly, such that the two criteria of the null linear discriminant analysis method are satisfied. Furthermore, we give a necessary and sufficient condition to guarantee full column rank of the orientation matrix, and describe the geometric characterization of this condition. Numerical experiments justify our theoretical analysis. Copyright © 2015 John Wiley & Sons, Ltd.