Skip to content

Kernel_Polynomial_Method_(KPM)

Jump to bottom
Dominique Delande edited this page Feb 20, 2023 · 23 revisions

The Kernel Polynomial Method (KPM in the following) is a numerical method used in the and-python package to compute the spectral function, the density of states and the temporal evolution of an arbitrary initial state. In a forthcoming version, it will also be used for fitering in energy an arbitrary state.

Useful references are:

[WWAF] is a good general reference on KPM. A slightly more detailed version is [WF]. It mainly concentrate on density of states and spectral funnctions.

For temporal evolution, see [RM], [HWMSD] and [HKM].

The thesis of Sanjib Ghosh [SG] contains relevant details on the implementation of the KPM method, also for energy filtering. Although Sanjib used a different software (a pure-C program with a slightly different implementation), most of what is written in [SG] is relevant for the and-python package.

The description which follows is essentially based on [WWAF], with contributions of [HWMSD], [HKM] and [SG] for the temporal propagation.

Basics of the KPM method

The starting point is that we want to compute the expectation value of an operator which depends on a disordered Hamiltonian H. The operator can be Tr(\delta(H-E)) if one wants to compute a density of states, a diagonal element of \delta(H-E) for a spectral function, or \exp(-iHt) for the evolution operator.

In general, computing such an operator is a difficult task. In our case, the operator H itself is rather simple in the spatial basis that we use: the disorder is diagonal and the kinetic energy operator (discretized Laplacian) has only few non-zero non-diagonal matrix elements. The basic idea is thus to express the operator which is a function of the Hamiltonian as a power series of the Hamiltonian. While, for e.g. the evolution operator \exp(-iHt), the Taylor series of the exponential could be used, this turns out to be not very efficient and hard to generalize to more complex operator.

The second key ingredient is the use of Chebyshev polynomials as approximations of arbitrary functions. In the simple case of a (sufficiently smooth) real function f(x) of a real variable x in the [-1,1] interval, this is known as the best approximation polynomial of the function f, and can be constructed explicitely, see [WWAF] section 2. One can write:

f(x) = \frac{1}{\pi\sqrt{1-x^2}} \left[ \mu_0 + 2 \sum_{n=1}^{\infty}{\mu_n\ T_n(x)}\right]

with coefficients:

\mu_n = \int_{-1}^{1}{f(x)\ T_n(x)\ \mathrm{d}x}

where T_n is the nth Chebyshev polynomial defined by:

T_n(x) = \cos[n \arccos(x)]

T_n is a polynomial of degree n, obeying the recursion relations:

T_0(x) = 1

T_{-1}(x) = T_1(x) = x

T_{n+1}(x) = 2xT_n(x) - T_{n-1}(x)

By truncating the series at a finite order, one obtains the best polynomial approximation (in some well defined sense which does not matter here) of the function f(x).

It is important to note that the expansion is convergent only in the [-1,1] interval. Outside this interval, the Chebyshev polynomials increase extremely rapidly, spoiling the convergence.

How can one use such an expansion for a disordered Hamiltonian H? Simply by using Chebyshev polynomials of the Hamiltonian itself.

Suppose for a moment that we have an Hamiltonian H' whose spectrum is entirely included in the [-1,1] interval, and that we want to compute the matrix element:

f(E) = \langle \beta | \delta(E-H') | \alpha \rangle

where \alpha and \beta are some known states in the Hilbert space [such a matrix element is the basis for computing a spectral function or a density of states). We can expand the function f(E) (non-zero only in the [-1,1] interval):

f(E) = \frac{1}{\pi\sqrt{1-E^2}} \left[ \mu_0 + 2 \sum_{n=1}^{\infty}{\mu_n\ T_n(E)}\right]

where the weights \mu_n are given by:

\mu_n = \int_{-1}^1 {\langle \beta | \delta(E-H') | \alpha\rangle\ T_n(E)\ dE} = \langle \beta | T_n(H') | \alpha \rangle

Thanks to the recursion relations between Chebyshev polynomials, the \mu_n can be computed recursively. Indeed, one has:

T_{n+1}(H')|\alpha\rangle = 2H' T_n(H')|\alpha\rangle - T_{n-1}(H') |\alpha \rangle

so that it is enough to keep in memory two consecutive vectors (i.e. states in the Hilbert space) T_{n-1}(H')|\alpha\rangle and T_n(H')|\alpha \rangle to generate the \mu_n.

The expensive part of the calculation is the product H' T_n(H')|\alpha\rangle, which involves to apply the Hamiltonian H' onto an arbitrary state T_n(H')|\alpha\rangle. In the spatial basis used in and-python, this is a sparse matrix-product operation which involves only few arithmetic operations per vector component. As it is also essentially local, it can be organized in a cache-friendly way.

Once the \mu_n are computed, the reconstruction of f(E) is straightforward.

More generally, independently of the states \alpha and \beta, one has the useful identity:

\delta(E-H') = \frac{1}{\pi\sqrt{1-E^2}} \left[ 1 + 2 \sum_{n=1}^{\infty}{T_n(H')\ T_n(E)}\right]

Gibbs oscillations

The expansion in Chebyshev polynmials shares common properties with the Fourier series. When truncating a Fourier series, it is well known that spurious oscillations - known as Gibbs oscillations - may take place. This is especially true when the functions involved are not sufficiently smooth. This is clearly the case for the spectral function which involves a singular \delta function. Thus, truncating abruptly the expansion of f(E) by setting all \mu_n for n>N to zero results in unwanted oscillations. A workaround is to "smoothly" damp the \mu_n coefficients, modifying the expansion of f(E) so that:

f(E) = \frac{1}{\pi\sqrt{1-E^2}} \left[ g_0 \mu_0 + 2 \sum_{n=1}^{N}{g_n\ \mu_n\ T_n(E)}\right]

where the g_n are conveniently chosen coefficients. The possible choices for the g_n's are discussed in [WWAF], section II.C. In the package and-python, we use exclusively the so-called Jackson kernel, which is defined by:

g_n = \frac{(N+n-1)\cos(\pi n/(N+1)) + \sin(\pi n/(N+1)) \cot (\pi/(N+1))}{N+1}

where N is the truncation order.

Basically, the effect of the Jackson kernel is to dams Gibbs oscillations and transform a \delta peak in a Gaussian of width proportional to 1/\sqrt{N}.

The Jackson kernel is systematically used in the and-python package, although different kernels should be easily added.

Mapping of the Hamiltonian

In general, the spectrum of the Hamiltonian is not included in the [-1,1]. Thus, in order to use the KPM, it is necessary to rescale the original Hamiltonian H. Following [WWAF], section II.B.1, we consider an interval [Emin,Emax] containing the spectrum of H. Emin and Emax MUST contain all the spectrum of H, that is Emin (resp. Emax) must be smaller (resp. larger) than the lowest (resp. largest) eigenvalue of H. We then define the shifted and rescaled operator H':

H' = (H-b)/a

where:

a = (Emax-Emin)/2

b = (Emax+Emin)/2

such that the spectrum of H' is entirely in the [-1,1] interval and its eigenstates are the same than those of H. One can then strightforwardly use the Chebyshev expansion for the operator H'. The expensive operation of multiplying an arbitrary state by the operator H' boils down to multiplication by H combined with a shift and a global multiplicative factor, the expensive part being the multiplication by H. In other words, the operation:

|\phi_{n+1}\rangle = 2H' |\phi_n\rangle - |\phi_{n-1}\rangle

is replaced by:

|\phi_{n+1}\rangle = 2H/a |\phi_n\rangle - 2b/a |\phi_n\rangle - |\phi_{n-1}\rangle

Temporal evolution

The evolution operator over time t is:

U(t) = \exp{-iHt}

which can be rewritten using the scaled Hamiltonian H' as:

U(t) = \exp{-ibt} U'(t'H')

where t'=at and U' is the evolution operator for H':

U'(t') = \exp{-it'H'}

Using a simple Fourier transform, one can calculate the evolution operator from the spectral properties of H':

U'(t') = \int_{-1}^{1} dE \ \delta(E-H') \exp{-iEt'}

which uses the fact that the spectrum of H' is in [-1,1]. By inserting the last equation of the Basis of the KPM method section, one obtains:

U'(t') = \int_{-1}^{1}{ dE\ \exp{-iEt'}\ \frac{1}{\pi \sqrt{1-E^2}} \left[ 1 + 2 \sum_{n=1}^{\infty} T_n(H') T_n(E) \right]}

The integral over E is simple:

c_n(t') = \int_{-1}^{1}{ dE\ \exp{-iEt'}\ \frac{1}{\pi \sqrt{1-E^2}} T_n(E)} = (-i)^n\ J_n(t')

where J_n is the standard Bessel function. One thus obtains:

U'(t') = c_0(t') + 2 \sum_{n=1}^{\infty} c_n(t') T_n(H')

which again involves Chebyshev polynomials of the Hamiltonian which can be efficiently computed.

Usually, one wants to compute the action of the evolution operator on some initial state |\psi_0\rangle, which boils down to:

U(t) |\psi_0\rangle = \exp{-ibt} \left[ J_0(at) |\psi_0\rangle + 2 \sum_{n=1}^{\infty} -i)^n \ J_n(at) |\phi_n\rangle

where the |\phi_n\rangle are computed recursively:

|\phi_{n+1}\rangle = 2H/a |\phi_n\rangle - 2b/a |\phi_n\rangle - |\phi_{n-1}\rangle

In order to save some memory bandwidth, the sum is computed using the Clenshaw trick.

The Bessel functions J_n(at) decrease quickly when n is large, so that the sum is rapidly converging, without any trick to avoid Gibbs oscillations. The code automatically truncates the series depending on the requested accuracy. Roughly, J_n(at) becomes exponentially small when n>|at|, so that the numer of terms in the series is essentially proportional to time t. If one wants to compute the evolution over a long time t, it is possible to do it directly with a large order in the Chebyshev expansion (which is numerically stable) or do it in p steps of duration t/p (with p an integer); altogether, the CPU time is almost the same for the two methods (avoid nevertheless too small steps).

Lower and upper energy bounds

As discussed above, it is of utmost importance to make sure that the energy spectrum of the Hamiltonian H is entirely contained in the interval [Emin,Emax]. Otherwise, any quantity computed with the KPM method will be exponentially diverging with the method order. If, in your calculation, you observe that the spectral function or the density of states takes huge values near the edges of [Emin,Emax], it is almost certainly because you underestimated Emax or overestimated Emin. Similarly, if the temporal evolution is not unitary (the norm of the state increases with time), almost surely an eigenvalue of H is not in the [Emin,Emax] interval.

In the absence of disorder, the energy spectrum is Emin=0 and Emax=2.0one_over_mass\sum_{i=1}^d 1.0/\delta_i^2. Note that is for a standard Hamiltonian and is more complex when spin-orbit interaction is added. A safe estimate in the presence of disorder is Emin=min(potential) and Emax= 2.0one_over_mass\sum_{i=1}^d 1.0/\delta_i^2 + max(potential). Tighter bounds can often be found, but depend on the disorder properties. If the disorder distribution does not have a long tail, the safe estimate is often sufficient. There is a routine energy_range to determine bounds for Emin and Emax. It has an optional argument accurate. When False, it uses the naive estimate. When True, it uses more accurate estimates, but costs a little more CPU time.

Home

Clone this wiki locally