1. Introduction and motivation
In this paper, we describe a simple fast algorithm for evaluating expressions of the form
where are real numbers, and are points in a compact interval of . This expression can be viewed as representing the electrostatic potential generated by charges on a line in . We remark that fast algorithms for computing the electrostatic potential generated by general distributions of charges in exist, see for example the Fast Multipole Method that was introduced by , and which has been extended by several authors including [7, 13]. However, in a number of situations in computational physics it is useful to have a simple and extremely fast method for evaluating the potential of charges on a line; we present such a method in this paper. Under mild assumptions the presented method involves operations and has a small constant. The method is based on writing the potential as
We show that there exists a small set of quadrature nodes and weights such that for a large range of values of we have
see Lemma 4.5, which is a quantitative version of (2). Numerically the nodes and weights can be computed using a procedure for constructing generalized Gaussian quadratures, see §5.2. An advantage of representing as a sum of exponentials is that the translation operator
can be computed by taking an inner product of the weights
with a diagonal transformation of the vector. Indeed, we have
The algorithm described in §3 leverages the existence of this diagonal translation operator to efficiently evaluate (1). We remark that previous works have used the diagonal form (4) of the translation operator (3) to evaluate expressions of the form (1), see Dutt, Gu and Rokhlin , and Yarvin and Rokhlin . The current paper improves upon these past works by taking advantage of robust generalized Gaussian quadrature codes  that were not previously available; in particular, the presented algorithm is both simpler and more practical for large numbers of points than past methods, see the discussion in §5.1.
Expressions of the form (1) appear in a number of situations in computational physics. In particular, such expressions arise in connection with the Hilbert Transform
For example, the computation of the projection of a function onto the first functions in a family of orthogonal polynomials can be reduced to an expression of the form (1) by using the Christoffel–Darboux formula, which is related to the Hilbert transform; we detail the reduction of to an expression of the form (1) in the following.
Let be a family of monic polynomials that are orthogonal with respect to the weight on . Consider the projection operator
where . Let and be the point Gaussian quadrature nodes and weights associated with , and set
By construction the polynomial that interpolates the valuesat the points will accurately approximate on when is sufficiently smooth, see for example §7.4.6 of Dahlquist and Björck . Directly evaluating (5) would require operations. In contrast, the algorithm of this paper together with the Christoffel–Darboux Formula can be used to evaluate (5) in operations. The Christoffel-Darboux formula states that
where we have used the fact that the diagonal term of the double summation is equal to . The summation in (7) can be rearranged into two expressions of the form (1), and thus the method of this paper can be used to compute a representation of in operations.
Analogs of the Christoffel–Darboux formula hold for many other families of functions; for example, if is a Bessel function of the first kind, then we have
A simple modification of the algorithm presented in this paper can be used to evaluate more general expressions of the form
where are source points, and are target points. For simplicity, this paper focuses on the case where the source and target points are the same, which is the case in the projection application described above.
2. Main result
2.1. Main result
Our principle analytical result is the following theorem, which provides precise accuracy and computational complexity guarantees for the algorithm presented in this paper, which is detailed in §3.
Let and be given. Set
Given and , the algorithm described in §3 computes values such
in operations, where
The proof of Theorem 2.1 is given in §4. Under typical conditions, the presented algorithm involves operations. The following corollary describes a case of interest, where the points are Chebyshev nodes for a compact interval (we define Chebyshev nodes in §4.2).
Fix , and let the points be Chebyshev nodes on . If , then the algorithm of §3 involves operations.
The proof of Corollary 2.2
is immediate from standard probabilistic estimates. The following remark describes an adversarial configuration of points.
Fix , and let be a collection of points such that and are evenly spaced in and , respectively, that is
We claim that Theorem 2.1 cannot guarantee a complexity better than for this configuration of points. Indeed, if , then , and if , then . In either case
This complexity is indicative of the performance of the algorithm for this point configuration; the reason that the algorithm performs poorly is that structures exist at two different scales. If such a configuration were encountered in practice, it would be possible to modify the algorithm of §3 to also involve two different scales to achieve evaluation in operations.
3.1. High level summary
The algorithm involves passing over the points twice. First, we pass over the points in ascending order and compute
and second, we pass over the points in descending order and compute
Finally, we define for such that
We call the computation of the forward pass of the algorithm, and the computation of the backward pass of the algorithm. The forward pass of the algorithm computes the potential generated by all points to the left of a given point, while the backward pass of the algorithm computes the potential generated by all points to the right of a given point. In §3.2 and §3.3 we give an informal and detailed description of the forward pass of the algorithm. The backward pass of the algorithm is identical except it considers the points in reverse order.
3.2. Informal description
In the following, we give an informal description of the forward pass of the algorithm that computes
Assume that a small set of nodes and weights such that
where By definition, the points are all distance at least from . Therefore, by (12)
If we define
then it is straightforward to verify that
Observe that we can update to using the following formula
We can now summarize the algorithm for computing . For each , we compute by the following three steps:
Update as necessary
Use to evaluate the potential from such that
Directly evaluate the potential from such that
By (16), each update of requires operations, and we must update at most times, so we conclude that the total cost of the first step of the algorithm is operations. For each , the second and third step of the algorithm involve and operations, respectively, see (15). It follows that the total cost of the second and third step of the algorithm is operations, where is defined in (9). We conclude that can be computed in operations. In §4, we complete the proof of the computational complexity guarantees of Theorem 2.1 by showing that there exist nodes and weights that satisfy (12), where is the approximation error in (12).
3.3. Detailed description
In the following, we give a detailed description of the forward pass of the algorithm that computes . Suppose that and are given and fixed. We describe the algorithm under the assumption that we are given quadrature nodes and weights such that
The existence of such weights and nodes is established in §4.4, and the computation of such nodes and weights is discussed in §5.2. To simplify the description of the algorithm, we assume that is a placeholder node that does not generate a potential.
Input: , . Output: .
update and :
for i = 1,…,m
compute potential from such that
compute potential from such that
In some applications, it may be necessary to evaluate an expression of the form (1) for many different weights associated with a fixed set of points . For example, in the projection application described in §1.2 the weights correspond to a function that is being projected, while the points are a fixed set of quadrature nodes. In such situations, pre-computing the exponentials used in the Algorithm 3.1 will significantly improve the runtime, see §5.1.
4. Proof of Main Result
In this section we complete the proof of Theorem 2.1; the section is organized as follows. In §4.2 we give mathematical preliminaries. In §4.3 we state and prove two technical lemmas. In §4.4 we prove Lemma 4.5, which together with the analysis in §3 establishes Theorem 2.1. In §4.5 we prove Corollary 2.1, and Corollary 2.2.
Let and be fixed, and suppose that , and are given. The interpolating polynomial of the function at is the unique polynomial of degree at most such that
This interpolating polynomial can be explicitly defined by
where is the nodal polynomial for , that is,
We say are Chebyshev nodes for the interval if
The following lemma is a classical result in approximation theory. It says that a smooth function on a compact interval is accurately approximated by the interpolating polynomial of the function at Chebyshev nodes, see for example §4.5.2 of Dahlquist and Björck .
Let , and be Chebyshev nodes for . If is the interpolating polynomial for at , then
In addition to Lemma 4.1, we require a result about the existence of generalized Gaussian quadratures for Chebyshev systems. In 1866, Gauss  established the existence of quadrature nodes and weights for an interval such that
whenever is a polynomial of degree at most . This result was generalized from polynomials to Chebyshev systems by Kreĭn . A collection of functions on is a Chebyshev system if every nonzero generalized polynomial
has at most distinct zeros in . The following result of Kreĭn says that any function in the span of a Chebyshev system of functions can be integrated exactly by a quadrature with nodes and weights.
Lemma 4.2 (Kreĭn ).
Let be a Chebyshev system of continuous functions on , and be a continuous positive weight function. Then, there exists unique nodes and weights such that
whenever is in the span of .
4.3. Technical Lemmas
Fix and , and let be Chebyshev nodes for . If is the interpolating polynomial for at , then
By writing the derivative of as
we can deduce that the maximum of occurs at , that is,
It remains to show that . Since is a increasing function, we have
Exponentiating both sides of this inequality gives , which is a classical inequality related to Stirling’s approximation. This completes the proof. ∎
Suppose that and are given. Then, there exists
values such that for all we have
for some choice of coefficients that depend on .
We construct an explicit set of points and coefficients such that (22) holds. Set . We define the points by
for and , and define the coefficients by
for and . We claim that
Indeed, fix , and let be the unique integer such that . By the definition of the coefficients, see (24), we have
We claim that the right hand side of this equation is the interpolating polynomial for at , that is,
Since the proof is complete. ∎
Indeed, in (24) the coefficients are either equal zero or equal to the nodal polynomial, see (19), for Chebyshev nodes on an interval that contains . The nodal polynomials for Chebyshev nodes on an interval are bounded by on , see for example . The fact that can be approximated as a linear combination of functions with small coefficients means that the approximation of Lemma 4.4 can be used in finite precision environments without any unexpected catastrophic cancellation.
4.4. Completing the proof of Theorem 2.1
points and weights that satisfy (17); we show the existence of such nodes and weights in the following lemma, and thus complete the proof of Theorem 2.1. The computation of such nodes and weights is described in §5.2.
Fix , and let and be given. Then, there exists nodes and weights such that
Fix , and let be given. By the possibility of rescaling , , and , we may assume that such that we want to establish (25) for . By Lemma 4.4 we can choose points , and coefficients depending on such that
whenever is in the span of . By the triangle inequality
Recall that we have assumed , in particular, so it follows that
By (26), the function can be approximated to error in the -norm on by functions in the span of . Since our quadrature is exact for these functions, we conclude that