{'fraction_non_alphanumeric': 0.047236249957057955, 'fraction_numerical': 0.011783297262015184, 'mean_word_length': 4.235611510791367, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Preconditioned iterative methods for eigenvalue counts

We describe preconditioned iterative methods for estimating the number of eigenvalues of a Hermitian matrix within a given interval. Such estimation is useful in a number of applications.In particular, it can be used to develop an efficient spectrum-slicing strategy to compute many eigenpairs of a Hermitian matrix. Our method is based on the Lanczos- and Arnoldi-type of iterations. We show that with a properly defined preconditioner, only a few iterations may be needed to obtain a good estimate of the number of eigenvalues within a prescribed interval. We also demonstrate that the number of iterations required by the proposed preconditioned schemes is independent of the size and condition number of the matrix. The efficiency of the methods is illustrated on several problems arising from density functional theory based electronic structure calculations.


Introduction
The problem of estimating the number of eigenvalues of a large and sparse Hermitian matrix A within a given interval [ξ , η] has recently drawn a lot of attention, e.g., [13,12]. One particular use of this estimation is in the implementation of a "spectrum slicing" technique for computing many eigenpairs of a Hermitian matrix [1,11]. Approximate eigenvalue counts are used to determine how to divide the desired spectrum into several subintervals that can be examined in parallel. In largescale data analytics, efficient means of obtaining approximate eigenvalue counts is required for estimating the generalized rank of a given matrix; see, e.g., [21].
A traditional approach for counting the number of eigenvalues of A in [ξ , η] is based on the Sylevester's law of inertia [14]. The inertia of the shifted matrices A − ξ I and A −ηI are obtained by performing LDL T factorizations of these matrices [1]. This approach, however, is impractical if A is extremely large or not given explicitly.
Several techniques that avoid factoring A have recently been described in [13,12]. These methods only require multiplying A with a number of vectors. In [12], a survey that describes several approaches to approximating the so-called density of states (DOS), which measures the probability of finding eigenvalues near a given point on the real line is presented. The DOS approximation can then be used to obtain an estimate of the number of eigenvalues in [ξ , η]. The potential drawback of a DOS estimation based approach is that, instead of directly targeting the specific interval [ξ , η], it always tries to approximate the eigenvalue distribution on the entire spectrum first.
Conceptually, the approaches in [13,12] are based on constructing a least-squares polynomial approximation of a spectral filter. Such approximations, however, often yield polynomials of a very high degree if A is ill-conditioned or the eigenvalues to be filtered are tightly clustered. These are common issues in practical large-scale computations. In particular, matrices originating from the discretization of partial differential operators tend to become more ill-conditioned as the mesh is refined. As a result, the polynomial methods of [13,12] can become prohibitively expensive. The overall cost of the computation becomes even higher if the cost of multiplying A with a vector is relatively high.
In this work we explore the possibility of using preconditioned iterative methods to reduce the cost of estimating the number of eigenvalues within an interval. By applying the Lanczos or Arnoldi iteration to preconditioned matrices with properly constructed Hermitian positive definite (HPD) preconditioners, we can significantly reduce the number of matrix-vector multiplications required to obtain accurate eigenvalue counts. Furthermore, when a good preconditioner is available, we can keep the number of matrix-vector multiplications (roughly) constant even as the problem size and conditioning of A increase. The methods we present in this paper do not require the lower and upper bounds of the spectrum of A to be estimated a priori. This feature compares favorably with the methods of [13,12] since obtaining such bounds can by itself be a challenging task. This paper is organized as following. Section 2 outlines the main idea, followed by derivation of the preconditioned Lanczos-type estimator based on Gauss quadrature in Section 3. The preconditioned Arnoldi-type algorithm is presented in Section 4. In Section 5, we discuss the proposed methods from the polynomial perspective. The performance of the introduced schemes depends to a large extent on the quality of the HPD preconditioner associated with the matrix A − τI. While the development of such a preconditioner is outside the scope of this paper, we point to several available options in Section 6. Several numerical experiments are reported in Section 7.

Basic idea
To simplify our presentation, let us assume that the endpoints ξ and η are different from any eigenvalue of A. Then the number of eigenvalues c(ξ , η) of A in [ξ , η] is given by the difference c(ξ , η) = n (A − ηI) − n (A − ξ I), where n (A − τI) denotes the negative inertia (i.e., the number of negative eigenvalues) of A − τI. Hence, in order to approximate c(ξ , η), it is sufficient to estimate n (A − τI) for a given real number τ.
The problem of estimating n (A − τI) can be reformulated as that of approximating the trace of a matrix step function. Namely, let Then Now let us assume that T is an HPD preconditioner for the shifted matrix A − τI in the sense that the spectrum of TA is clustered around a few distinct points on the real line. Specific options for constructing such preconditioners will be discussed in Section 6.
If T is available in a factorized form If T = MM * is chosen in such a way that its spectrum has a favorable distribution, i.e., the eigenvalues of C is clustered in a few locations, then estimating trace {h(C)} can be considerably easier than estimating trace {h(A − τI)} If the multiplication of C with a vector can be performed efficiently, then the trace of C can be estimated as where the entries of each vector v j are i.i.d. random variables with zero mean and unit variance; see [10,2]. It follows that for a sufficiently large sample size m. The variance of the stochastic trace estimator is known to depend on the magnitude of off-diagonal entries of the considered matrix, which is h(C) in (5). Clearly, different choices of the preconditioned operator C yield different matrices h(C), and hence lead to different convergence rates of the estimator (5).

Preconditioned Lanczos
If A is large, then the exact evaluation of h(C) in (5) can be prohibitively expensive, because it requires a full eigendecomposition of the preconditioned matrix. A more practical approach in this situation would be to (approximately) compute v * h(C)v for a number of randomly sampled vectors v without explicitly evaluating the matrix function.

The Gauss quadrature rule
Let us assume that T = M * M is available in the factorized form and let C = M(A − τI)M * in (5). We also assume that the Hermitian matrix C has p ≤ n distinct eigenvalues µ 1 < µ 2 < . . . < µ p .
Consider the orthogonal expansion of v in terms of the eigenvectors of C, i.e., v = ∑ p i=1 α i u i , where u i is an normalized eigenvector associated with the eigenvalue µ i , and α i = u * i v. It is then easy to verify that where p denotes the number of negative eigenvalues. The right-hand side in (6) can be viewed as a Stieltjes integral of the step function h with respect to the measure defined by the piecewise constant function Therefore, using (7), we can write (6) as Computing the above integral directly is generally infeasible because the measure (7) is defined in terms of the unknown eigenvalues of C. Nevertheless, the right-hand side of (8) can be approximated by using the Gauss quadrature rule [6], where the k nodes θ 1 ≤ θ 2 ≤ . . . ≤ θ k and weights w 1 , w 2 , . . . , w k of the quadrature are determined from k steps of the Lanczos procedure (see Algorithm 1) applied to the preconditioned matrix C with the starting vector v. In (9), k denotes the number of negative nodes θ i . Tridiagonal matrix J k+1,k and the Lanczos basis Q k+1 = [q 1 , q 2 , . . . , q k+1 ]. 1:

7: end for
Specifically, given q 1 = v/ v , running k steps of the Lanczos procedure in Algorithm 1 yields the relation where J k+1,k is the tridiagonal matrix The eigenvalues of the leading k × k submatrix of J k+1,k , denoted by J k , are ordered so that Then the Gauss quadrature rule on the right-hand side of (9) is defined by eigenvalues and eigenvectors of J k , i.e., where z i is the eigenvector of J k associated with the eigenvalue θ i , z i (1) denotes its first component [6], and k denotes the number of negative Ritz values. If the preconditioner T = MM * is chosen in such a way that the spectrum of and will be constant elsewhere. Hence, the integral in (8) will be determined only by integration over [a, b] because h vanishes in [c, d]. Therefore, in order for quadrature rule (9) to be a good approximation to (8), its nodes should be chosen inside [a, b].
In the extreme case in which clustered eigenvalues of C coalesce into a few eigenvalues of higher multiplicities, the number of Lanczos steps required to obtain an accurate approximation in (12) is expected to be very small.
Then the Gauss quadrature (12) will be exact with at most k = p nodes.
where u i is an eigenvector of C associated with the eigenvalue µ i . Then p steps of Lanczos process with v as a starting vector produce a tridiagonal matrix J p and an orthonormal basis Q p , such that the first column of and, by comparing with (6), we see that the quadrature (12) gives the exact value Proposition 1 implies that in the case of an ideal preconditioner, where M(A−τI)M * has two distinct eigenvalues, the Gauss quadrature rule (12) is guaranteed to be exact after at most two Lanczos steps.

The algorithm
Let J ( j) k denote the k-by-k tridiagonal matrix resulting from the k-step Lanczos procedure applied to C = M(A − τI)M * with a random starting vector v j . Assume that k j is the number of its negative eigenvalues. Then, by (5) and (12), the quantity n (A − τI) can be approximated from the estimator where z ( j) i (1) denotes the first components of a normalized eigenvector z associated with the negative eigenvalues. It is expected that, for a sufficiently large m, L τ (k, m) ≈ n (A − τI). The expression (13) is what Algorithm 2 uses to estimate the number of eigenvalues of A that are to the left of τ.
In order to estimate the number of eigenvalues in a given interval [ξ , η], Algorithm 2 should be applied twice with τ = ξ and τ = η. The difference between the estimated n (A − ξ I) and n (A − ηI) yields the desired count. The two runs of Algorithm 2 generally require two different HPD preconditioners, one for A − ξ I and the other for A − ηI. In some cases, however, it can be possible to come up with a single preconditioner that works well for both runs.
The cost of Algorithm 2 is dominated by computational work required to perform the preconditioned matrix-vector multiplication of M(A − τI)M * v at each iteration of the Lanczos procedure. The eigenvalue decomposition of the tridiagonal matrix J k , as well as reorthogonalization of the Lanczos basis in step 6 of Algorithm 1, is Algorithm 2: The preconditioned Lanczos-type estimator for n (A − τI)

Input:
Matrix A, shift τ, HPD preconditioner T = M * M for A − τI, number of steps k, and parameter m. Output: approximate number C τ of eigenvalues of A that are less than τ; Generate v ∼ N (0, I).

4:
Run k steps of Lanczos process in Algorithm 1 with the starting vector v to obtain tridiagonal matrix J k .

5:
Find the eigendecomposition (Θ , Z) of J k . Let z 1 , . . . , z k be unit eigenvectors associated with negative eigenvalues. 6: negligibly small for small values of k, which can be ensured by a sufficiently high quality preconditioner. Note that, in exact arithmetic, the Lanczos basis Q i should be orthonormal [14]. However, in practice, the orthogonality may be lost; therefore, we reorthogonalize Q i at every iteration of Algorithm 1.

Bias of the estimator
A relation between the Gauss quadrature (12) and where ε k is the error of the quadrature rule. Thus, (13) can be written as where ε ( j) k denotes the error of the quadrature rule for v * j h(C)v j . As m increases, the first term in the right-hand side of (14) converges to trace {h(C)} = n (A − τI). Thus, L τ (k, m) is a biased estimate of n (A−τI), where the bias is determined by the (average) error of the quadrature rule, given by the second term in the right-hand side of (14). In other words, the accuracy of L τ (k, m) generally depends on how well the Gauss quadrature captures the value of the matrix functional Bounds on the quadrature error for a matrix functional v * f (C)v, where f is a sufficiently smooth function and C is a Hermitian matrix, are well known. In particular, the result of [3] gives the bound where the constant N k is such that | f (2k) (x)| ≤ N k for x in the interval containing spectrum of C, and β j are the off-diagonal entries of (11). (1) is discontinuous. Therefore, bound (15)  , which contains zero, we leth(x) to be a polynomial p(x) of degree 4k + 1, such that p(b) = 1, p(c) = 0, and p (l) (b) = p (l) (c) = 0 for l = 1, . . . , 2k. This choice of polynomial will ensure that the piecewise functioñ h(x) is 2k times continuously differentiable. (Note that p(x) can always be be constructed by (Hermite) interpolation with the nodes b and c; see, e.g., [15].) We then apply (15) to obtain the bound on the quadrature error for v * h (M(A − τI)M * )v. As discussed above, this yields the estimate of the error ε k of quadrature rule (12) for functional v * h(M(A − τI)M * )v. Thus, we can conclude that the latter is bounded by (15), where N k is the maximum of |p (2k) (x)| on the interval [b, c].
This finding shows that we can expect that (12)  Fortunately, a good choice of the preconditioner T = M * M can ensure that eigenvalues of M(A − τI)M * are clustered and away from zero. In this case, the Ritz values typically converge rapidly to these eigenvalues after a few Lanczos steps. Thus, with a good preconditioner, the Gauss quadrature (12) can effectively approximate the matrix functional v * h(M(A − τI)M * )v, yielding small errors ε k for a relatively small number of quadrature nodes. As a result, the bias of the estimator L τ (k, m) in (14) will be small and, as confirmed by numerical experiments in Section 7.

The generalized averaged Gauss quadrature rule
The Gauss quadrature rule (12) is exact for all polynomials of degree at most 2k − 1; e.g., [6].
In the recent work of [16] (and references therein), a so-called generalized averaged (GA) Gauss quadrature rules was introduced. This quadrature rule make use of the same information returned by a k-step Lanczos process, but gives an exact integral value for polynomials of degree 2k. Hence it is more accurate at essentially the same cost.
When applying the GA Gauss quadrature rule to the matrix functional v * h(C)v in (8), we still use the expression (12), except that we have (2k − 1) nodes θ 1 , θ 2 , . . . , θ 2k−1 which are the eigenvalues of the matrix obtained from J k+1,k in (11) by extending its tridiagonal part in a "reverse" order. The set (α i ) of numbers in (16) gives the diagonal entries of J k+1,k , whereas (β i ) define the upper and lower diagonals. Similarly, the associated weights w i are determined by squares of the first components of the properly normalized eigenvectors z i ofJ 2k−1 associated with the eigenvalues θ i ; see [16] for more details. Thus, we can expect to increase accuracy of the estimator by a minor modification of Algorithm 2. This modification will only affect step 5 of the algorithm, where J k must be replaced by the extended tridiagonal matrix (16).

Preconditioned Arnoldi
Sometimes, the preconditioner T is not available in a factored form T = MM * . In this case, it may be necessary to work with T (A − τI) or (A − τI)T directly. One possibility is to make use of the fact that (A − τI)T is self adjoint with respect to an inner product induced by T . This property allows us to carry out a T -inner product Lanczos procedure that produces Similarly, we can use a T −1 -inner product based Lanczos procedure to obtain where Y k = M * Q k . Even though it may appear that we do not need T in a factored form in either (17) or (18), the starting vectors we use to generate (17) and (18) are related to M. In particular, (17) must be generated from x 1 = M −1 q 1 and (18) must be generated from y 1 = M * q 1 , where q 1 is a random vector with i.i.d entries.
Another approach is to construct an estimator based on (5), where C = T (A−τI). This will require evaluating the bilinear form v * h(T (A − τI))v, where h is a function of a matrix T (A − τI) that has real spectrum but is non-Hermitian in standard inner product. Similar to the Hermitian case, the matrix functional v * h(T (A − τI))v can be viewed as an integral, such that (19) where Γ is a contour that encloses the spectrum of T (A − τI) and the bar denotes complex conjugation; see, e.g., [9]. This integral can be approximated by a quadrature rule based on a few steps of the Arnoldi process (Algorithm 3) applied to the preconditioned operator T (A − τI) with a starting vector v [4,6].

Algorithm 3: The Arnoldi procedure for T (A − τI)
Input: Matrix A − τI, HPD preconditioner T , starting vector v, and number of steps k.

8: end for
Given q 1 = v/ v , Algorithm 3 produces an orthonormal Arnoldi basis Q k+1 and an extended upper Hessenberg matrix such that T (A − τI)Q k = Q k+1 H k+1,k , Q * k+1 Q k+1 = I. An Arnoldi quadrature rule for the integral (19) is fully determined by the k-by-k leading submatrix H k of (20). Similar to (12) where w i are determined by the first components of the (right) eigenvectors z 1 , . . . , z k of H k associated with its k eigenvalues that have negative real parts, and t i is the ith entry of the first column of S = Z −1 . Similar to Proposition 1, it can be shown that if T (A − τI) has p distinct eigenvalues, then (21) is exact with at most p nodes. Let H ( j) k be the upper Hessenberg matrix produced by the Arnoldi process applied to C = T (A − τI) with the starting vector v j . Then (21) and (5) yield the estimator where z k . Similar to (13), we expect that, for a sufficiently large m, the real part of A τ (k, m) approximates n (A − τI). The computation of Re (A τ (k, m)) is described in Algorithm 4.

Input:
Matrix A, shift τ, HPD preconditioner T for A − τI, number of steps k, and parameter m. Output: approximate number A τ of eigenvalues of A that are less than τ; Generate v ∼ N (0, I).

4:
Run k steps of Arnoldi process in Algorithm 3 with the starting vector v to obtain upper Hessenberg matrix H k .

6:
The cost of Algorithm 4 is comparable to that of Algorithm 2, and is slightly higher mainly due to the need to invert the eigenvector matrix of H k . In contrast to Algorithm 2, the above described scheme assumes complex arithmetic, because the upper Hessenberg matrix H k is non-Hermitian and can have complex eigenpairs. However, for good choices of T , the imaginary parts tend to be small in practice as, for a sufficiently large k, the eigenpairs of H k converge rapidly to those of T (A−τI), which are real. Finally, note that the derivation of the estimator (22) assumes an extension of the definition of the step function (1), such that h(x) has the value of one on the left half of the complex plane, and is zero elsewhere. ). Then, we can replace h(C) in (5) by a polynomial approximation p l (C) of degree l. There are several ways to choose this polynomial. One option is to take p l (C) as formal truncated expansion of h(x) in the basis of Chebyshev polynomials. This choice is related the approach described in [13].

Polynomial viewpoint
The quality of a polynomial approximation p l (x) of h(x) can be measured by the difference between p l (x) and h(x) on the set of eigenvalues of C. When the spectrum of C has an arbitrary distribution, constructing a polynomial that provides the best least squares fit on the entire interval containing all eigenvalues, as is done in [13], is well justified.
When a good preconditioner is used, the spectrum of C tends to cluster around several points on the real line. Thus, a natural approach would be to choose p l such that it is only close to h in regions that contain eigenvalue clusters. It can be quite different from h elsewhere. An example of such an approach is an interpolating polynomial, e.g., [15], that interpolates h at eigenvalue clusters. A practical construction of such a polynomial is given by the following theorem, which relates the the interpolation procedure to the Lanczos or Arnoldi process.
Theorem 1 (see [17,8]). Let Q k , T k be the orthonormal basis and the projection of the matrix C generated from a k-step Lanczos (Arnoldi) process, with the starting vector v. Then where p k−1,v is the unique polynomial of degree at most k − 1 that interpolates f in the Hermite sense on the spectrum of T k .
The subscript "v" in p k−1,v is used to emphasize the dependence of the polynomial on the staring vector v. Note that T k is a symmetric tridiagonal matrix if C is Hermitian. It is upper Hessenberg otherwise.
Using formula (23), it is easy to verify that if C = M(A − τI)M * , then the bilinear form v * p k−1,v (C)v is exactly the same as the Gauss quadrature rule on the right-hand side of (12). Similarly, if C = T (A − τI), then v * p k−1,v (C)v is given by the Arnoldi quadrature on the right-hand side of (21). Hence, both estimators (13) and (22) can be viewed as a stochastic approximation of trace p k−1,v (C) , where p k−1,v (x) is an interpolating polynomial of degree k − 1 for the step function h.

Preconditioning
The iterative scheme we presented earlier rely on the assumption that the operator T is HPD, as this property guarantees that the inertia of the original matrix A − τI is preserved after preconditioning. Furthermore, a good choice of T should cluster spectrum of the preconditioned matrix C around several points in the real axis.
An ideal HPD preconditioner will result in the preconditioned matrix with only two distinct eigenvalues. In this case, by Proposition 1, the Lanczos procedure should terminate in two steps. An example of such an ideal preconditioner is the matrix T = |A − τI| −1 , where the absolute value is understood in the matrix function sense.
Clearly, the choice T = |A − τI| −1 is prohibitively costly in practice. However, it is possible to construct HPD preconditioners that only approximate |A − τI| −1 . Such a preconditioning strategy was proposed in [19] and is referred to as the absolute value (AV) preconditioning. It was shown in [19] that, e.g., for discrete Laplacian operators, AV preconditioners can be efficiently constructed using multigrid (MG).
Another possible option is to employ the incomplete LDL T (ILDL) factorization. Given a matrix A − τI and a drop tolerance t, an ILDL(t) preconditioner is of the form T = L − * D −1 L −1 , where L is lower triangular and D is block-diagonal with diagonal blocks of size 1 and 2, such that T ≈ (A − τI) −1 .
Clearly, since A − τI is indefinite, the ILDL(t) procedure will generally result in an indefinite T , which cannot be applied within the preconditioned estimators of this paper. Therefore, we suggest to modify it by taking the absolute value of diagonal blocks of D, so that T = L − * |D| −1 L −1 . Such a preconditioner is HPD, and the cost of the proposed modification is marginal. This idea has been motivated by [5], where a similar approach was used in the context of full (complete) LDL T factorization.
Finally, in certain applications, HPD operators are readily available and traditionally used for preconditioning indefinite matrices. For example, this is the case in Density Functional Theory (DFT) based electronic structure calculations in which the solutions are expressed in terms of a linear combination of planewaves. A widely used preconditioner, often referred to as the Teter preconditioner [18], is diagonal in the planewave basis.

Numerical experiments
We now study the numerical behavior of the proposed methods for three test problems listed in Table 1. The matrix "Laplace" represents a standard five-point finite differences (FD) discretization of the 2D Laplacian on a unit square with mesh size h = 2 −7 . The problems "Benzene" and "H2" originates from the DFT based electronic structure calculations. The former is a FD discretization of a Hamiltonian operator associated with a ground state benzene molecule 1 , whereas the latter corresponds to a Hamiltonian associated with the hydrogen molecule generated by the KSSOLV package [20]. Throughout, our goal is to estimate the quantity n (A − τI) for a given value of the shift τ.  Table 1 Estimates of n (A − τI) produced by Algorithm 2 with different preconditioners and the corresponding numbers k of Lanczos iterations for three test problems. Table 1 presents the results of applying the Lanczos-type estimator given in Algorithm 2 to the test problems with different preconditioner choices. For the "Laplace" and "Benzene" matrices, we use the positive definite ILDL(t) based preconditioning with different drop tolerance t, discussed in the previous section. The ILDL factorizations of A − τI are obtained using the sym-ildl package [7]. In the "H2" test, we employ the diagonal Teter preconditioner available in KSSOLV. In both cases, the preconditioner is accessible in the factorized form T = M * M. The number of random samples m is set to 50 in all tests.
In the table, we report estimates of n (A − τI) produced by Algorithm 2 along with the corresponding numbers of Lanczos iterations (k) performed at each sampling step. The reported values of k correspond to the smallest numbers of Lanczos iterations that result in a sufficiently accurate estimate. The error associated with these approximations have been observed to be within 5%. Table 1 demonstrates that the use of preconditioning significantly reduces the number of Lanczos iterations. Furthermore, k becomes smaller as the quality of the preconditioner, which is controlled by the drop tolerance t in the ILDL(t) based preconditioners, improves for the "Laplace" and "Benzene" tests.  Figure 1 shows that the quality of the estimates can be further improved by using the GA Gauss quadrature rules discussed in Section 3.4. In both plots, the horizontal axis corresponds to the number of Lanczos iterations (k) per sampling step, and the vertical axis is the corresponding estimate of n (A − τI). It can be seen that the estimator based on the GA Gauss quadrature (referred to as "GA Lanczos") is generally more accurate for the two test problems, with the accuracy difference being especially evident for smaller values of k.
In the context of linear systems arising from discretizations of partial differential equations, an important property of preconditioning is that it allows maintaining the same number of iterations needed to obtain solution regardless of problem size. A similar phenomenon can be observed when estimating n (A − τI) using the preconditioned methods of this paper.  Table 2 Independence of preconditioned Arnoldi-and Lanczos-type estimators for n (A − τI) on the discretization parameter for the "Laplace" (left) and "H2" (right) problems.
In Table 2 (left) we consider a family of discrete Laplacians, whose size and condition numbers increase as the mesh parameter h is refined. For each of the matrices, we apply the Arnoldi-type estimator of Algorithm 4 with the MG AV preconditioner from [19] and, similar to above, report the smallest numbers k of Arnoldi iterations per sampling step needed to obtain a sufficiently accurate estimate (within 5% error) of n (A − τI). The results are compared against those of an unpreconditioned estimator based on (5), where C = A and the step function h(A) is replaced by its least-squares polynomial approximation of degree k constructed using the basis of Chebyshev polynomials. The latter (referred to as "Chebyshev") is essentially the approach proposed in [13].
It can be seen from the table, that Algorithm 4 with the AV preconditioner exhibits behavior that is independent of h. Regardless of the problem size and conditioning, the number of Arnoldi steps stays (roughly) the same (between 16 and 19).
In Table 2 (right) we report a similar test for a sequence of "H2" problems obtained by increasing the kinetic energy cutoff (ecut) from 25 to 125 Ry in the plane wave discretization. This gives Hamiltonian matrices with sizes ranging from 1, 024 to 23, 583. Again, we observe that the behavior of the Lanczos-type estimator in Algorithm 2 with the Teter preconditioner [18] is essentially independent of the discretization parameter, whereas the "Chebyshev" approach tends to require higher polynomial degrees as the problem size grows.