Large scale density estimation with polynomial expansions and Johnson-Lindenstrauss transformations

Abstract

The computation of the charge density is a ubiquitous and time-consuming task in large scale electronic structure computations, e.g in Nanodevice Simulations, Density Functional Theory and the solution of Self-Consistent Field equations. In this work, we present a method for computing the charge density in a system which is based on techniques related to the Kernel Polynomial Method (KPM) and Johnson-Lindenstrauss transformations (JLT), a fundamental concept in the dimensionality reduction literature. As it is known, the charge density can be directly derived by a spectral decomposition of the Hamiltonian matrix, which gives the wave functions and the corresponding energy levels. However, a direct diagonalization scales proportionally to the cube of the number of discretization points, that is, O(n3) where n is the number of discretization points in the domain. For large systems this cost can be prohibitive, and many different alternative methods have been proposed in the literature to overcome this limitation. The charge density can be derived from the diagonal elements of the matrix function f(H; EF), where f(x, EF) is the Fermi distribution with Fermi Level EF. At zero temperature, the Fermi distribution is simply a low-pass filter, which truncates all energies above the Fermi level. Similar to the KPM, one can approximate the Fermi distribution with a Chebyshev polynomial expansion. Therefore, the problem reduces to evaluating the diagonal of a matrix polynomial. However, if the matrix polynomial is not evaluated carefully, the total cost can be larger than that of the diagonalization, since one needs to compute several matrix powers, which also scale as O(n3). To reduce the overhead of evaluating the matrix polynomial we exploit the aforementioned JLT. Combining these techniques we derive a very simple approximation algorithm, with asymptotic complexity O(n2 × log(n) × M / ε2), where M is the degree of the polynomial and ε∈(0,½) is the desired relative error accuracy. The asymptotic complexity is near-optimal, up to constant and logarithmic factors. Numerical experiments are also presented, demonstrating the properties of the proposed method in practice.