Automatic, high-order, and adaptive algorithms for Brillouin zone integration

We thank the referee for a thoughtful and helpful report. Below, we give a detailed point by point response to the referee's comments, and also indicate our intended modifications to the manuscript. In particular, although the referee makes several interesting and useful points on the PTR which we will certainly mention in our revised manuscript, we believe that we have presented a highly efficient implementation of the PTR scheme and have thus provided a fair comparison with the IAI proposal.

(1) For 1D and 2D model functions, we feel that the benefit of adaptive integration over equispaced integration for functions with localized features is conceptually clear, and standard material. Since the degree of advantage is example-specific, we feel that additional comparisons in model cases would distract from the main points of the manuscript while increasing its length unnecessarily. An exception is Fig. 4, which, using a simple model, demonstrates a central and non-trivial conceptual point of the manuscript: that for functions with a specific structure often encountered in BZ integration, IAI is asymptotically superior to TAI. The PTR and IAI are compared for a real example in Figs. 5 and 7.

However, we will certainly add a remark comparing the number of integrand evaluations required by the PTR and adaptive integration to achieve a given accuracy, for our 1D model problem with a small value of eta.

(2) Algorithm 2 relies on a fixed 1D quadrature rule for smooth functions on an interval. The referee proposes to replace the Gauss rule (the optimal rule for non-periodic functions on an interval, exponentially convergent in p as recalled by Thm. 2) by a p-node nonperiodic trapezoidal rule (which is inferior, being merely $O(1/p^2)$ convergent). The poor accuracy of the latter scheme on each panel would result in a very large number of panels needed by Alg. 2 to achieve a small error tolerance. This would not be a fair comparison against the PTR applied globally to the periodic integrand, which in contrast to the above second-order accurate scheme, is exponentially convergent for analytic functions.

Our manuscript focuses only on high-order accurate methods, a necessary ingredient (from the point of view of efficiency) in algorithms which allow the user to specify a possibly small error tolerance. We will include a remark which explains this.

The other suggestion is to consider a global Gauss quadrature rule rather than the PTR. This would be inferior to the PTR, because the Gauss rule does not exploit the known global periodicity, whereas the PTR does. While both are exponentially convergent schemes for analytic functions, the Gauss rule is known to be the gold standard for an interval without a periodicity assumption, whereas the PTR is the gold standard for a periodic interval. This is standard material in numerical analysis, covered in Refs. [26,28,39]. Indeed, it can be shown that the convergence rate of the PTR is a factor $\pi/2$ faster than that of Gauss quadrature for a periodic function that is known only to be analytic in a given strip $|\text{Im}\, z| < a$. See, e.g., "New quadrature formulas from conformal maps," N. Hale and L. N. Trefethen, SIAM J. Numer. Anal. 46(2), 930--948 (2008). Thus there would be no point in testing the Gauss rule in the global periodic setting in which it is provably inferior.

We will make a remark addressing this point, and include a reference.

(3) The referee makes the interesting observation that, if one were starting from an $N_{trial}$ and $\Delta N$ that were smaller by a large factor than the $N_{max}$ at which the desired accuracy is achieved, then an arithmetic progression would require $O(N_{max}^2)$ evaluations and would therefore be inefficient. However, this is not the situation that we describe. Rather, we choose $\Delta N$ via a rigorous estimate that scales as $1/\eta$, as detailed in App. A, so that only $O(1)$ adaptation steps are taken. Thus we stand by the decision to use an arithmetic rather than a geometric progression of grid sizes in our efficient implementation of the automatic PTR.

To elaborate, since the PTR converges exponentially, once one is in the convergence regime, adding a fixed number $\Delta N$ of points decreases the error by approximately a fixed factor. Thus, a good strategy is to estimate a number $N_{trial}$ sufficient to be in the convergence regime (but not overkill), and then to repeatedly add $\Delta N$ points in order to drive the error down exponentially below the tolerance. As a consequence of the exponential convergence, this is expected to take at most a fixed, constant number of steps, giving a total number $O(N_{max})$ of points calculated, not $O(N_{max}^2)$. We will add a remark in the manuscript clarifying how $N_{trial}$ is chosen.

If a geometric progression of grid sizes were used, this would in certain cases result in significant over-convergence, and therefore a significantly larger cost. For example, if one chose $\alpha = 2$, one would often encounter situations like the following: $N = 1100$ is sufficient for convergence below $\epsilon$, but one has $N_{trial} = 1000$ and $N_{test} = 1000 \times 2 = 2000$. In this case, one does $2000^3 / 1100^3 \approx 2^3 = 8$ times more work than is needed. Using an arithmetic progression of grid sizes, with spacing chosen using $\eta$, avoids this common scenario.

The benefit of reusing previously computed values is only valid for discrete choices of $\alpha$, the smallest of which is $\alpha = 2$. Therefore, in 3D, the final calculation will be so much more expensive than the previous ones that the savings gained by reusing the computed values will be insignificant.

(4) This is a very interesting suggestion and has useful advantages: it avoids storing $H(k)$ (which also increases the flops-to-DRAM-access ratio), and it avoids using over-resolved grids for frequencies which may only require a coarse grid. However, this algorithm has the disadvantage that it requires specifying all frequency evaluation points in advance. In our case, we use adaptive interpolation in the frequency domain, and therefore our frequency evaluation points are determined on the fly.

If one wishes to evaluate the Green's function, for example, at a pre-specified frequency grid, the algorithm proposed by the referee may be preferable. However, if one wishes to compute the Green's function to a user-specified accuracy $\epsilon$, as is the goal in our manuscript, the frequency grid must resolve features at $O(\eta)$ scale which are, in general, at unknown locations. Therefore, an a priori-specified frequency grid would require $O(1/\eta)$ points. We avoid this using an adaptive algorithm, described in Sec. V, which determines an efficient frequency grid of $O(\log(1/\eta))$ points on the fly, and is therefore not compatible with the algorithm proposed by the referee. We do not feel the possible advantage gained by switching to this algorithm is sufficient to justify replacing our $O(\log(1/\eta))$ scaling with $O(1/\eta)$ scaling.

We do not believe that our goal of computing the Green's function automatically to a user-specified tolerance is contrived. Although in some cases one only wants, for example, to roughly plot the density of states, in other cases the calculation of the Green's function appears in the inner loop of a complicated procedure (as in dynamical mean-field theory). In that case it may be required to perform the calculation automatically within a required error tolerance. We therefore feel that both algorithms have their place, and will make a remark explaining this in our revised manuscript.

(5) The suggestion of using a pruned FFT algorithm is interesting, but we do not feel that it is likely to give a significant advantage over the method we have proposed, for the following reasons.

Let $n_R$ denote the number of Fourier coefficients of $H(k)$ per dimension, with $N_R = n_R^3$, and let $n_k$ denote the number of $k$-evaluation points per dimension, with $N_k = n_k^3$. Although it is stated in the given reference (Tsirkin, npj Comput. Mater. 7, 33 (2021)) that direct evaluation of the Fourier series has a cost scaling as $O(N_R N_k)$, this is only the case if one evaluates the full 3D Fourier series at each $k$-point. However, since the $k$-points fall on a Cartesian grid, this is unnecessary and inefficient. Rather, as described in Sec IV.B (subsec. "Rapid evaluation of $H(k)$ in iterated adaptive integration") and in Appendix C, one should split the Fourier series dimension-by-dimension, so that the dominant cost is that of evaluating a one-dimensional Fourier series of $n_R$ terms at each of $N_k$ points. Each such Fourier series evaluation is extremely efficient when the geometric series property (addition formula) is exploited, requiring a single complex exponential evaluation, $2n_R$ multiplications, and $n_R$ additions. The computational complexity is therefore $O(N_k n_R)$, not $O(N_k N_R)$. This method also allows us to take full advantage of all crystal symmetries, as explained in Appendix C.

The proposed pruned FFT algorithm has complexity $O(N_k \log n_R)$. While this is asymptotically less than $O(N_k n_R)$, it is difficult to state a priori which method will be faster in practice, since $n_R$ is small and the direct evaluation algorithm has an extremely small prefactor. Furthermore, it appears that the pruned FFT algorithm does not necessarily take advantage of all symmetries for typical grid densities, which could be a significant disadvantage in practice.

Therefore, while a direct comparison of these approaches would be an interesting topic of future research, we feel that we have implemented a competitive and likely optimal Fourier series evaluation algorithm, and do not wish to distract from the main points of the paper by including this. We will make a remark in the text referring to the pruned FFT algorithm, and a summary of the points made above.

(6) The TAI method was eliminated in Section IV on a variety of grounds, and is highly unlikely to be competitive with IAI in practice. Indeed, Fig. 4 shows for a model problem that TAI requires more integrand evaluations than IAI, even for large $\eta$, and each integrand evaluation is unlikely to be more efficient in TAI than in IAI. Furthermore, integration on the irreducible Brillouin zone would require designing a TAI scheme with tetrahedral cells, and high-order quadratures on tetrahedra. This would be a very significant undertaking, and makes it difficult to include TAI, which is unlikely to be competitive in any case, in the comparisons in Sec. VI. Neither would it be fair to run the cube TAI on the full Brillouin zone and compare with IAI and PTR on the irreducible Brillouin zone.

(7) The PTR is expected to be faster than IAI for large values of $\eta$ for a few reasons. First, for large $\eta$, one does not need much adaptivity. In this case, IAI will produce panels of roughly equal size; if they were equal, this would yield a method of order $2p$ (for us, $p = 4$), as opposed to the spectrally convergent PTR. For smooth, periodic functions which do not require $k$-space adaptivity, the PTR is an optimal, highly efficient method.

Second, if one evaluates at multiple frequency points, the PTR can reuse previously stored values of $H(k)$, whereas for IAI $H(k)$ needs to be evaluated on the fly. As we state in our manuscript, this is a very significant advantage. To make a comparison which is sufficiently generous to the PTR, the cost of evaluating $H(k)$ was not even included in the PTR timings in Fig. 5, whereas obviously for IAI this evaluation does contribute to the timings (we will state this in the revised manuscript).

Lastly, the IAI method comes with some overhead which the much simpler PTR does not. The PTR is significantly easier to optimize at a low level, and we put considerable work into optimizing our implementation of the PTR in order to provide a fair comparison.

The number of function evaluations used by IAI and the PTR is not a fair metric for comparison of performance, because the cost per function evaluation is significantly different for the two methods. The most significant reason is that in IAI, $H(k)$ must be evaluated on the fly. For this reason, we have opted to only show comparisons using wall clock timings.

(8) This is a good idea, and we will try to add a corner zoom to these figures.

We again thank the referee for these thoughtful comments, which will lead to significant improvements in the clarity of the manuscript. We believe our responses should alleviate any concerns that we have not presented a high-performance implementation of the PTR.