Krylov Space Methods Research Articles

Converging Time-Dependent Density Functional Theory Calculations in Five Iterations with Minimal Auxiliary Preconditioning.

Eigenvalue problems and linear systems of equations involving large symmetric matrices are commonly solved in quantum chemistry using Krylov space methods, such as the Davidson algorithm. The preconditioner is a key component of Krylov space methods that accelerates convergence by improving the quality of new guesses at each iteration. We systematically design a new preconditioner for time-dependent density functional theory (TDDFT) calculations based on the recently introduced TDDFT-ris semiempirical model by retuning the empirical scaling factor and the angular momenta of a minimal auxiliary basis. The final preconditioner produced includes up to d-functions in the auxiliary basis and is named "rid". The rid preconditioner converges excitation energies and polarizabilities in 5-6 iterations on average, a factor of 2-3 faster than the conventional diagonal preconditioner, without changing the converged results. Thus, the rid preconditioner is a broadly applicable and efficient preconditioner for TDDFT calculations.

Convection experiments with the exponential time integration scheme

We present the construction and preliminary tests of the exponential time integration schemes for the Euler equations. The equations are written in terms of potential temperature and Exner function. The time stepping is accomplished with the EPI2 and EPI3 schemes which have been used in the past for very accurate, long time step integration of the shallow water equations on the sphere. The stability of exponential integrators of this class is ensured by approximation of the oscillatory term using an absolutely stable exponential formula. Consequently, the schemes are not subject to the CFL condition imposed by the linear oscillatory part of the system and they also eliminate the phase errors associated with all the known implicit time stepping algorithms. All calculations are performed using the phipm algorithm based on the Krylov space methods. It is shown that the spectral radius of the Jacobian and consequently the Krylov space dimension are very sensitive to the presence of spatial noise. Efficiency and accuracy increase after the Shapiro filter is applied every two time steps. The method is investigated using several standardized and well accepted benchmark tests, including convective bubbles and a cold density current proposed by Straka and collaborators. Numerical experiments show that a stable integration is obtained for the wave Courant number of the order of 60 to 250. We present a detailed discussion of the development of a cold density current capturing the spiraling Prandtl vortex and the subsequent appearance of the Kelvin–Helmholtz billows. Our results from the exponential method compare very well with those obtained from the accurate schemes reported in the recent literature. The algorithm accurately handles both advection as well as acoustic and gravity waves. These results warrant the further development of the scheme for modeling atmospheric processes on a wide range of scales.

Krylov Methods for Large Sparse Systems: A Comprehensive Overview

In this paper are analyzed behavior and properties for different Krylov methods applied in different categories of problems. These categories often include PDEs, econometrics and network models, which are represented by large sparse systems. For our empirical analysis are taken into consideration size, the density of non-zero elements, symmetry/un-symmetry, eigenvalue distribution, also well/ill-conditioned and random systems. Convergence, approximation error and residuals are compared for the full version of methods, some restarted methods and preconditioned methods. Two preconditioners are considered respectively, ILU(0) and IC(0) by using at least five preconditioning techniques. In each case, empirical results show which technique is best to use based on properties of the system and are backed up by general theoretical information already found on Krylov space methods.

Block Modified Gram--Schmidt Algorithms and Their Analysis

New block modified Gram--Schmidt (BMGS) methods for the Q-R factorization of a full column rank matrix $X \in {{R}}^{m \times n}$, $m \geq n$, are considered. Such methods factor $X$ into $Q \in {{R}}^{m \times n}$ and an upper triangular $R \in {{R}}^{m \times n}$ such that $X = QR,$ where, in exact arithmetic, $Q$ is left orthogonal (i.e., $Q^T\;{Q}=I_n$). Gram--Schmidt-based algorithms play an important role in the implementation of Krylov space methods, such as GMRES, Arnoldi, and Lanczos. For these applications, the left orthogonal factor $Q$ is needed and the matrix is produced either one column at a time or one block of columns at a time. For block implementations of Krylov methods, a block of columns of $X$ is introduced at each step, and a new block of columns of $Q$ must then be produced. That is a task for which BMGS methods are ideally suited. However, for these Krylov methods to converge properly, the BMGS algorithms need to have numerical behavior that is similar to that of modified Gram--Schmidt. To design such BMGS algorithms, we build upon the block Householder representation of Schreiber and Van Loan [SIAM J. Sci. Stat. Comput., 10 (1989), pp. 53--57] and an observation by Charles Sheffield analyzed by Paige [SIAM J. Matrix Anal. Appl., 31 (2009), pp. 565--583] about the relationship between modified Gram--Schmidt and Householder Q-R factorization. Our new BMGS algorithms exploit the Sheffield framework so that they share a similar relationship to Householder Q-R and thus have error analysis properties similar to modified Gram--Schmidt. The last BMGS algorithm developed is based entirely upon matrix multiplications and the “tall, skinny” Q-R (TSQR) factorization---two operations that have been studied extensively for cache-based architectures and distributed architectures. It is shown that if the TSQR part of the BMGS algorithm satisfies error analysis properties connected to the Sheffield structure, then so does the entire BMGS algorithm. Thus new criteria for when a BMGS algorithm has error analysis properties similar to those of MGS are proposed.

A Charge-Conservative Finite Element Method for Inductionless MHD Equations. Part II: A Robust Solver

In [L. Li, M. Ni, and W. Zheng, SIAM J. Sci. Comput., 41 (2019), pp. B796--B815] a charge-conservative finite element method is proposed for solving inductionless and incompressible magnetohydrodynamic (MHD) equations. The purpose of this paper is to propose a robust solver for the discrete problem. Using the framework of field-of-values-equivalence, we first study the preconditioned Krylov space method for the continuous problem in the setting of Hilbert spaces. The algebraic preconditioner for the discrete problem is then obtained by representing the preconditioner for the continuous problem in finite element spaces. By three numerical examples, the optimality of the solver to the number of unknowns is demonstrated for both stationary and time-dependent MHD problems.

Accelerating the iterative solution of convection–diffusion problems using singular value decomposition

SummaryThe discretization of convection–diffusion equations by implicit or semi‐implicit methods leads to a sequence of linear systems usually solved by iterative linear solvers such as the generalized minimal residual method. Many techniques bearing the name of recycling Krylov space methods have been proposed to speed up the convergence rate after restarting, usually based on the selection and retention of some Arnoldi vectors. After providing a unified framework for the description of a broad class of recycling methods and preconditioners, we propose an alternative recycling strategy based on a singular value decomposition selection of previous solutions and exploit this information in classical and new augmentation and deflation methods. The numerical tests in scalar nonlinear convection–diffusion problems are promising for high‐order methods.

Accelerating molecular property calculations with nonorthonormal Krylov space methods.

We formulate Krylov space methods for large eigenvalue problems and linear equation systems that take advantage of decreasing residual norms to reduce the cost of matrix-vector multiplication. The residuals are used as subspace basis without prior orthonormalization, which leads to generalized eigenvalue problems or linear equation systems on the Krylov space. These nonorthonormal Krylov space (nKs) algorithms are favorable for large matrices with irregular sparsity patterns whose elements are computed on the fly, because fewer operations are necessary as the residual norm decreases as compared to the conventional method, while errors in the desired eigenpairs and solution vectors remain small. We consider real symmetric and symplectic eigenvalue problems as well as linear equation systems and Sylvester equations as they appear in configuration interaction and response theory. The nKs method can be implemented in existing electronic structure codes with minor modifications and yields speed-ups of 1.2-1.8 in typical time-dependent Hartree-Fock and density functional applications without accuracy loss. The algorithm can compute entire linear subspaces simultaneously which benefits electronic spectra and force constant calculations requiring many eigenpairs or solution vectors. The nKs approach is related to difference density methods in electronic ground state calculations and particularly efficient for integral direct computations of exchange-type contractions. By combination with resolution-of-the-identity methods for Coulomb contractions, three- to fivefold speed-ups of hybrid time-dependent density functional excited state and response calculations are achieved.

Auxiliary master equation approach within matrix product states: Spectral properties of the nonequilibrium Anderson impurity model

Within the recently introduced auxiliary master equation approach it is possible to address steady state properties of strongly correlated impurity models, small molecules or clusters efficiently and with high accuracy. It is particularly suited for dynamical mean field theory in the nonequilibrium as well as in the equilibrium case. The method is based on the solution of an auxiliary open quantum system, which can be made quickly equivalent to the original impurity problem. In its first implementation a Krylov space method was employed. Here, we aim at extending the capabilities of the approach by adopting matrix product states for the solution of the corresponding auxiliary quantum master equation. This allows for a drastic increase in accuracy and permits us to access the Kondo regime for large values of the interaction. In particular, we investigate the nonequilibrium steady state of a single impurity Anderson model and focus on the spectral properties for temperatures T below the Kondo temperature T_K and for small bias voltages phi. For the two cases considered, with T=T_K/4 and T=T_K/10 we find a clear splitting of the Kondo resonance into a two-peak structure for phi close above T_K. In the equilibrium case (phi=0) and for T=T_K/4, the obtained spectral function essentially coincides with the one from numerical renormalization group.

Fast computation of orthonormal basis for RBF spaces through Krylov space methods

In recent years, in the setting of radial basis function, the study of approximation algorithms has particularly focused on the construction of (stable) bases for the associated Hilbert spaces. One of the ways of describing such spaces and their properties is the study of a particular integral operator and its spectrum. We proposed in a recent work the so-called WSVD basis, which is strictly connected to the eigen-decomposition of this operator and allows to overcome some problems related to the stability of the computation of the approximant for a wide class of radial kernels. Although effective, this basis is computationally expensive to compute. In this paper we discuss a method to improve and compute in a fast way the basis using methods related to Krylov subspaces. After reviewing the connections between the two bases, we concentrate on the properties of the new one, describing its behavior by numerical tests.

COMPLEXITY OF FIRST-ORDER METHODS FOR DIFFERENTIABLE CONVEX OPTIMIZATION

This is a short tutorial on complexity studies for differentiable convex optimization. A complexity study is made for a class of problems, an that obtains information about the problem at a given point, and a stopping rule for algorithms. These three items compose a scheme, for which we study the performance of algorithms and problem complexity. Our problem classes will be quadratic minimization and convex minimization in ℝn. The oracle will always be first order. We study the performance of steepest descent and Krylov spacemethods for quadratic function minimization and Nesterov’s approach to the minimization of differentiable convex functions.

A robust all-at-once multigrid method for the Stokes control problem

In this paper we present an all-at-once multigrid method for a distributed Stokes control problem (velocity tracking problem). For solving such a problem, we use the fact that the solution is characterized by the optimality system (Karush-Kuhn-Tucker-system). The discretized optimality system is a large-scale linear system whose condition number depends on the grid size and on the choice of the regularization parameter forming a part of the problem. Recently, block-diagonal preconditioners have been proposed, which allow to solve the problem using a Krylov space method with convergence rates that are robust in both, the grid size and the regularization parameter or cost parameter. In the present paper, we develop an all-at-once multigrid method for a Stokes control problem and show robust convergence, more precisely, we show that the method converges with rates which are bounded away from one by a constant which is independent of the grid size and the choice of the regularization or cost parameter.

Experiments on a Parallel Nonlinear Jacobi–Davidson Algorithm

The Jacobi–Davidson (JD) algorithm is very well suited for the computation of a few eigen-pairs of large sparse complex symmetric nonlinear eigenvalue problems. The performance of JD crucially depends on the treatment of the so-called correction equation, in particular the preconditioner, and the initial vector. Depending on the choice of the spectral shift and the accuracy of the solution, the convergence of JD can vary from linear to cubic. We investigate parallel preconditioners for the Krylov space method used to solve the correction equation.We apply our nonlinear Jacobi–Davidson (NLJD) method to quadratic eigenvalue problems that originate from the time-harmonic Maxwell equation for the modeling and simulation of resonating electromagnetic structures.

Large scale micro finite element analysis of 3D bone poroelasticity

In this paper, a solver for poroelasticity problems related to osteoporotic human bones is discussed. Osteoporosis is a major health problem that compromises the integrity of bones. A good understanding of the disease requires an accurate simulation of the physics. For that purpose, a finite element solver based on Biot's consolidation equations has been developed. A mixed formulation is used to discretize the geometries taken from medical imaging. The resulting indefinite linear systems are solved by Krylov space methods supplemented by variants of Schur complement-based block preconditioners.

Eigenvalue computations based on IDR

The induced dimension reduction (IDR) method, which has been introduced as a transpose-free Krylov space method for solving nonsymmetric linear systems, can also be used to determine approximate eigenvalues of a matrix or operator. The IDR residual polynomials are the products of a residual polynomial constructed by successively appending linear smoothing factors and the residual polynomials of a two-sided (block) Lanczos process with one right-hand side and several left-hand sides. The Hessenberg matrix of the OrthoRes version of this Lanczos process is explicitly obtained in terms of the scalars defining IDR by deflating the smoothing factors. The eigenvalues of this Hessenberg matrix are approximations of eigenvalues of the given matrix or operator.

Eigenvalue Computations Based on IDR

The induced dimension reduction (IDR) method, which has been introduced as a transpose-free Krylov space method for solving nonsymmetric linear systems, can also be used to determine approximate eigenvalues of a matrix or operator. The IDR residual polynomials are the products of a residual polynomial constructed by successively appending linear smoothing factors and the residual polynomials of a two-sided (block) Lanczos process with one right-hand side and several left-hand sides. The Hessenberg matrix of the OrthoRes version of this Lanczos process is explicitly obtained in terms of the scalars defining IDR by deflating the smoothing factors. The eigenvalues of this Hessenberg matrix are approximations of eigenvalues of the given matrix or operator.

Reorthogonalized block classical Gram–Schmidt

A reorthogonalized block classical Gram---Schmidt algorithm is proposed that factors a full column rank matrix $$A$$ into $$A=QR$$ where $$Q$$ is left orthogonal (has orthonormal columns) and $$R$$ is upper triangular and nonsingular. This block Gram---Schmidt algorithm can be implemented using matrix---matrix operations making it more efficient on modern architectures than orthogonal factorization algorithms based upon matrix-vector operations and purely vector operations. Gram---Schmidt orthogonal factorizations are important in the stable implementation of Krylov space methods such as GMRES and in approaches to modifying orthogonal factorizations when columns and rows are added or deleted from a matrix. With appropriate assumptions about the diagonal blocks of $$R$$ , the algorithm, when implemented in floating point arithmetic with machine unit $$\varepsilon _M$$ , produces $$Q$$ and $$R$$ such that $$\Vert I- Q ^T\!~ Q \Vert =O(\varepsilon _M)$$ and $$\Vert A-QR \Vert =O(\varepsilon _M\Vert A \Vert )$$ . The first of these bounds has not been shown for a block Gram---Schmidt procedure before. As consequence of these results, we provide a different analysis, with a slightly different assumption, that re-establishes a bound of Giraud et al. (Num Math, 101(1):87---100, 2005) for the CGS2 algorithm.

A thin-sheet model for global electromagnetic induction

We consider the frequency domain electromagnetic (EM) induction problem in a layered spherical Earth overlain by an inhomogeneous thin sheet of finite conductance, and derive an integral equation for the secondary electric fields resulting from surface conductivity variations. The Green′s function kernel of the integral equation is singular, and must in practice be regularized for a numerical treatment. We compare two approaches: truncation of the spherical harmonic representation of the Green′s function, and projection onto a discrete spatial grid, with elements of the resulting discrete system matrix evaluated as contour integrals. In either case rotational symmetry and the fast Fourier transform (FFT) may be exploited to reduce computational costs, and Krylov space methods may be used to solve the resulting discretized integral equation efficiently. We illustrate the methods for several source field geometries. The two regularization schemes result in solutions that agree well at large scales, but exhibit substantial differences as the nominal resolution is approached. The agreement improves with increased resolution. Although the level of solution errors for the two approaches is likely comparable in an L2 sense, the truncation approach results in Gibb′s oscillations near sharp conductivity contrasts. These may be reduced by tapering the truncated kernel, essentially equivalent to low-pass filtering the solution.

Higher order Galerkin time discretizations and fast multigrid solvers for the heat equation

We discuss numerical properties of continuous Galerkin–Petrov and discontinuous Galerkin time discretizations applied to the heat equation as a prototypical example for scalar parabolic partial differential equations. For the space discretization, we use biquadratic quadrilateral finite elements on general two-dimensional meshes. We discuss implementation aspects of the time discretization as well as efficient methods for solving the resulting block systems. Here, we compare a preconditioned BiCGStab solver as a Krylov space method with an adapted geometrical multigrid solver. Only the convergence of the multigrid method is almost independent of the mesh size and the time step leading to an efficient solution process. By means of numerical experiments we compare the different time discretizations with respect to accuracy and computational costs.

Expository Research Papers

The two papers in this issue are concerned with fair division and with approximation of functions. “Divide-and-Conquer: A Proportional, Minimal-Envy Cake-Cutting Algorithm,” by Steven Brams, Michael Jones, and Christian Klamler, has nothing to do with wedding cakes or Food Network Challenges. Instead it deals with the mathematical problem of proportional division, which often arises in social science applications. Let us illustrate the cake-cutting paradigm with the divide-and-choose algorithm for two players. Suppose you and I are to cut a cake, so that each of us gets half of the cake, where “half” is interpreted according to our preferences. Suppose you like chocolate and prefer a piece with at least half of all the chocolate decoration. I, on the other hand, like marzipan and prefer a piece with at least half of all the marzipan decoration. In the divide-and-choose algorithm, I start by cutting the cake into two pieces, so that each piece contains the same amount of marzipan. Then you pick the piece with the most chocolate. This is a proportional and envy-free division of the cake. The division is proportional because each of us gets at least half of what we want (chocolate in your case, marzipan in mine), and it is envy-free because none of us would want the other's piece. Mathematically, one expresses the preference for chocolate or marzipan in terms of a probability density function. For n players the situation is much more complicated. The authors propose a divide-and-conquer algorithm that requires only the minimal number of $n-1$ cuts. Although the algorithm is not envy-free, it does minimize the maximal number of players that any player can envy. The algorithm is also simple and easy to implement. Now for the second paper. Back in 1901 the German mathematician Carl Runge showed that polynomial interpolation of a function at n equally spaced points in the interval $[-1,1]$ can produce interpolants that fail to converge as $n\rightarrow\infty$, because they exhibit large-amplitude oscillations near the interval endpoints. In their paper “Impossibility of Fast Stable Approximation of Analytic Functions from Equispaced Samples,” Rodrigo Platte, Lloyd Trefethen, and Arno Kuijlaars look at the situation in a more general setting, and with finite precision computation in the back of their minds. Now the approximants do not have to be interpolating polynomials. They can be much more general. All that is required is that an approximant $\phi_n(f)$ depend on f at the n equally spaced grid points in $[-1,1]$, but nowhere else. The question is: How can perturbations in f affect the convergence of the approximant $\phi_n(f)$ in the worst case? To answer this question, the authors introduce a condition number $\kappa_n$ that measures how sensitive $\phi_n(f)$ can be to changes in f. You can think of $\kappa_n$ as sort of a Lipschitz constant for $\phi_n$. The main result is: If the approximants converge exponentially fast, in the sense that $\|f-\phi_n(f)\|\sim \sigma^n$ for some $\sigma<1$, then $\kappa_n\sim C^n$ for some $C>1$. This means, in the worst case, exponentially fast convergence comes at the price of exponential ill-conditioning. The authors also identify connections with potential theory, matrix iterations in the form of Krylov space methods, and quadrature formulae.

Complex dispersion relation calculations with the symmetric interior penalty method

AbstractA high‐order discontinuous Galerkin method for calculations of complex dispersion relations of two‐dimensional photonic crystals is presented. The medium is characterized by a complex‐valued permittivity and we relate for this absorptive system the spectral parameter to the time frequency. We transform the non‐linear eigenvalue problem for a Lorentz material in air into a non‐Hermitian linear eigenvalue problem and uses a Krylov space method to compute approximate eigenvalues. Moreover, we study the impact of the penalty term numerically and illustrate the high convergence rate of the method. Copyright © 2010 John Wiley &amp; Sons, Ltd.