The standard approach for computing the trace of the inverse of a very large, sparse matrix A is to view the trace as the mean value of matrix quadratures, and use the Monte Carlo algorithm to estimate it. This approach is heavily used in our motivating application of Lattice QCD. Often, the elements of A −1 display certain decay properties away from the non zero structure of A, but random vectors cannot exploit this induced structure of A −1 . Probing is a technique that, given a sparsity pattern of A, discovers elements of A through matrix-vector multiplications with specially designed vectors. In the case of A −1 , the pattern is obtained by distance-k coloring of the graph of A. For sufficiently large k, the method produces accurate trace estimates but the cost of producing the colorings becomes prohibitively expensive. More importantly, it is difficult to search for an optimal k value, since none of the work for prior choices of k can be reused.First, we introduce the idea of hierarchical probing that produces distance-2 i colorings for a sequence of distances 2 0 , 2 1 , . . . , 2 m , up to the diameter of the graph. To achieve this, we do not color the entire graph, but at each level, i, we compute the distance-1 coloring independently for each of the node-groups associated with a color of the distance-(2 i−1 ) coloring. Second, based on this idea, we develop an algorithm for uniform, toroidal lattices that simply applies bit-arithmetic on local coordinates to produce the hierarchical permutation. Third, we provide an algorithm for choosing an appropriate sequence of Hadamard and Fourier vectors, so that earlier vectors in the sequence correspond to hierarchical probing vectors of smaller distances. This allows us to increase the number of probing vectors until the required accuracy is achieved.Several experiments show that when a decay structure exists in the matrix, our algorithm finds it and approximates the trace incrementally, starting with the most important contributions. We have observed up to an order of magnitude speedup over the standard Monte Carlo.
We present a new algorithm that computes eigenvalues and eigenvectors of a Hermitian positive definite matrix while solving a linear system of equations with conjugate gradient (CG). Traditionally, all the CG iteration vectors could be saved and recombined through the eigenvectors of the tridiagonal projection matrix, which is equivalent theoretically to unrestarted Lanczos. Our algorithm capitalizes on the iteration vectors produced by CG to update only a small window of vectors that approximate the eigenvectors. While this window is restarted in a locally optimal way, the CG algorithm for the linear system is unaffected. Yet, in all our experiments, if the window has more than a properly chosen but small number of vectors, it converges to the required eigenvectors at a rate identical to unrestarted Lanczos. After the solution of the linear system, eigenvectors that have not accurately converged can be improved in an incremental fashion by solving additional linear systems. In this case, eigenvectors identified in earlier systems can be used to deflate, and thus accelerate, the convergence of subsequent systems. We have used this algorithm with excellent results in lattice quantum chromodynamics applications, where hundreds of right-hand sides may be needed. Specifically, about 70 eigenvectors are obtained to full accuracy after solving 24 right-hand sides. Deflating these from the large number of subsequent right-hand sides removes the dreaded critical slowdown, where the conditioning of the matrix increases as the quark mass reaches a critical value. Our experiments show almost a constant number of iterations for our method, regardless of quark mass, and speedups of 8 over original CG for light quark masses. Introduction.The numerical solution of linear systems of equations of large, sparse matrices is central to many scientific and engineering applications. One of the most computationally demanding applications is lattice quantum chromodynamics (QCD) because not only does it involve very large matrix sizes but it also requires the solution of several linear systems with the same matrix but different right-hand sides. Direct methods, although attractive for multiple right-hand sides, cannot be used because of the size of the matrix. Iterative methods provide the only means for solving these problems.QCD is the theory of the fundamental force known as strong interaction, which describes the interactions among quarks, one of the constituents of matter. Lattice QCD is the tool for nonperturbative numerical calculations of these interactions on a Euclidean space-time lattice [58]. The heart of the computations is the solution of the lattice-Dirac equation, which translates to a linear system of equations M x = b, often for a large number of right-hand sides [20]. The Dirac operator M is
Abstract. Large, sparse, Hermitian eigenvalue problems are still some of the most computationally challenging tasks. Despite the need for a robust, nearly optimal preconditioned iterative method that can operate under severe memory limitations, no such method has surfaced as a clear winner. In this research we approach the eigenproblem from the nonlinear perspective that helps us develop two nearly optimal methods. The first extends the recent Jacobi-Davidson-Conjugate-Gradient (JDCG) method to JDQMR, improving robustness and efficiency. The second method, Generalized-Davidson+1 (GD+1), utilizes the locally optimal Conjugate Gradient recurrence as a restarting technique to achieve almost optimal convergence. We describe both methods within a unifying framework, and provide theoretical justification for their near optimality. A choice between the most efficient of the two can be made at runtime. Our extensive experiments confirm the robustness, the near optimality, and the efficiency of our multimethod over other state-of-the-art methods.1. Introduction. The numerical solution of large, sparse, Hermitian and real symmetric eigenvalue problems is central to many applications in science and engineering. It is also one of the most time consuming tasks. Recently, electronic structure calculations, with eigenproblems at their core, have displaced Quantum Chromodynamics as the top supercomputer cycle user. The symmetric eigenvalue problem seems deceptively simple to solve, given well conditioned eigenvalues and a wealth of theoretical knowledge. Yet, these advantages have enabled researchers to push modeling accuracy to unprecedented levels, routinely solving for a few extreme eigenvalues of matrices of size more than a million, while an order of a billion has also been tried [48].The sheer size of these problems can only be addressed by iterative methods. At the same time, the size imposes limits on the memory available to these methods. Moreover, preconditioning becomes imperative to reduce the total number of iterations. While many eigenvalue iterative methods have been proposed, there is little consensus on which method is the best and in what situations. The unrestarted Lanczos method is known to be optimal for solving Hermitian eigenvalue problems but, unlike the Conjugate Gradient (CG) method for linear systems, it requires unlimited storage of its iteration vectors. With preconditioning, or under limited storage, the question of optimality remains open. Furthermore, there is a noticeable scarcity of high quality, general purpose software for preconditioned eigensolvers. In this research we seek an optimal, or nearly optimal, method that can utilize preconditioning and that can be implemented in a robust software that is also flexible.In the particular case of seeking one eigenpair, if the eigenvalue were known, solving a linear system using CG to obtain the corresponding eigenvector would yield the optimal Lanczos convergence. In practice both the eigenvalue and the eigenvector are unknown, so the appropriate w...
The Davidson method is a popular preconditioned variant of the Arnoldi method for solving large eigenvalue problems. For theoretical, as well as practical reasons the two methods are often used with restarting. Frequently, information is saved through approximated eigenvectors to compensate for the convergence impairment caused by restarting. We call this scheme of retaining more eigenvectors than needed`thick restarting', and prove that thick restarted, non-preconditioned Davidson is equivalent to the implicitly restarted Arnoldi. We also establish a relation between thick restarted Davidson, and a Davidson method applied on a de ated system. The theory is used to address the question of which and how many eigenvectors to retain and motivates the development of a dynamic thick restarting scheme for the symmetric case, which can be used in both Davidson and implicit restarted Arnoldi. Several experiments demonstrate the e ciency and robustness of the scheme. and its JD modi cation can be regarded as two ways of improving convergence and robustness at the expense of a more complicated step.Often, eigenvalue problems are very large and ill-conditioned. As a result, eigenvalue methods require a large number of steps, and need to save all the vector iterates for extracting the eigenvectors. Such cases exhibit overwhelming storage requirements. In addition, the Lanczos and Arnoldi processes, which traditionally had been considered without restarting, are plagued by orthogonality problems and spurious solutions. For the above reasons many restarting variants are used in practice 7, 18, 24, 2]. The GD method improves convergence, and solves some of the aforementioned problems through orthogonalization and preconditioning. However, the number of iterations can still grow very large, and cause similar storage problems. The problem is actually aggravated in the symmetric case, where the better theoretical framework and software has led researchers to consider matrices of huge size that allow only a few vectors to be stored. The GD method also can be restarted every time the basis contains m vectors (GD(m)). If the l lowest eigenvalues are needed, the l lowest Ritz values are computed at the m th step, and their corresponding Ritz vectors are used as initial guesses for the restarted GD iteration.Truncating the Krylov sequence is expected to impair the convergence rate of the method. There are two main reasons: the new vectors entering the basis repeat some of the information that was discarded when restarting, and the Rayleigh-Ritz procedure does not minimize over the whole Krylov subspace. There has been much discussion about the problems caused by restarted methods for both linear systems and eigenvalue problems 27, 20, 24]. Some methods tend to save additional information at each restart 14, 3, 11]. For the Davidson method, Murray et al. 17], and Van Lenthe et al. 29], have proposed restarting with two vectors per required Ritz vector with some success. In an e ort to minimize execution time, Crouzeix et al. 6], have pr...
scite is a Brooklyn-based organization that helps researchers better discover and understand research articles through Smart Citations–citations that display the context of the citation and describe whether the article provides supporting or contrasting evidence. scite is used by students and researchers from around the world and is funded in part by the National Science Foundation and the National Institute on Drug Abuse of the National Institutes of Health.
hi@scite.ai
10624 S. Eastern Ave., Ste. A-614
Henderson, NV 89052, USA
Copyright © 2024 scite LLC. All rights reserved.
Made with 💙 for researchers
Part of the Research Solutions Family.