An adaptive parallel algorithm for computing zeros of analytic functions

Technical Summary: An Adaptive Parallel Algorithm for Computing Zeros of Analytic Functions

Problem Statement
The paper addresses the challenge of robustly and efficiently computing the zeros of analytic functions within a specified rectangular region of the complex plane. While root-finding for non-linear equations is a fundamental problem in physics, engineering, and mathematics, standard methods often struggle with global convergence or fail to provide certified results for functions with complex behaviors, such as exponential polynomials and Dirichlet polynomials. Specifically, the authors note that finite values of a function and its derivatives are insufficient to determine the number of zeros in a bounded region without global data, and discretizing contours for argument counting requires careful handling to avoid missing zeros or accumulating round-off errors.

Methodology
The proposed algorithm integrates four distinct components to ensure robustness, convergence, and speed:

1. Adaptive Argument Counting (Global Information):
   The core of the zero-counting procedure relies on the Argument Principle. To compute the variation of the argument along a path, the authors adapt a criterion derived from Ying and Katz. This method discretizes the contour into segments where the variation of the argument is guaranteed to be less than $\pi$.

   • Mechanism: For a segment $[s_1, s_2]$, the function $f(s)$ is approximated by an affine function $P(s)$. The algorithm ensures that the minimum modulus of $P(s)$ exceeds the maximum error bound of the remainder $R(s) = f(s) - P(s)$.
   • Error Bounds: The error bound is derived using a bound on the second derivative, $|f''(s)| \le M(s_1, s_2)$, via a corrected interpolation error lemma.
   • Adaptivity: An adaptive backtracking line search is employed. If the condition for the argument variation is not met (e.g., due to a zero nearby or a loose bound), the segment is subdivided. This avoids the need for a priori fixed discretization and handles regions where the argument varies rapidly.

2. Two-Dimensional Bisection Strategy:
   Once the number of zeros in a rectangle is determined, the algorithm recursively subdivides the region to isolate individual zeros.

   • The rectangle is split by bisecting its longest side.
   • To handle zeros lying on or near the subdivision boundary, the algorithm perturbs the subdivision segment and recomputes the argument variation, ensuring no zeros lie on the boundaries of the resulting sub-rectangles.
   • This process yields a collection of size-balanced rectangles, each containing exactly one simple zero, with diameters converging linearly to zero.

3. Exclusion Criterion:
   To improve efficiency, the algorithm employs an exclusion criterion based on Taylor's formula. For a rectangle $R$ with center $z_0$ and radius $r$, if the inequality $|f(z_0)| - |f'(z_0)|r - \frac{1}{2}N(z_1, z_2)r^2 > 0$ holds (where $N$ bounds the second derivative), the rectangle is discarded as containing no zeros. This step is optional but significantly reduces the search space.

4. Acceleration and Certification (Newton–Kantorovich):
   Since the bisection method converges linearly, the algorithm accelerates the final stage using Newton's method.

   • Once a rectangle is isolated to contain a single zero, Newton's method is applied starting from the rectangle's center.
   • Certification: The convergence is certified using the Kantorovich theorem. The algorithm checks conditions involving the inverse of the derivative and the bound on the second derivative to guarantee that a unique solution exists within a specific radius and that Newton's method will converge to it.
   • Global Convergence: If Newton's method fails (e.g., diverges or exceeds iteration limits), the algorithm reverts to the bisection method on the current rectangle, ensuring global convergence.

Implementation and Parallelization
The algorithm is implemented in C++ using the Boost libraries for multiprecision arithmetic.

• Parallel Architecture: The subdivision of the initial rectangle creates independent subproblems. The implementation uses a stack-based approach where available processors pull rectangles from a mutex-protected stack.
• Thread Safety: Mutexes protect the task stack and the "answers" stack (where certified zeros are stored) to prevent race conditions.
• Precision: The use of multiprecision libraries allows the algorithm to handle the high precision required for functions like the Riemann zeta function partial sums.

Key Results
The authors applied the algorithm to two specific classes of functions:

1. Exponential Polynomials: These arise in delay-differential equations and exhibit retarded, advanced, and neutral asymptotic types. The algorithm successfully computed zeros for a test function exhibiting all three types.
2. Dirichlet Polynomials (Riemann Zeta Partial Sums): The primary benchmark involved computing the zeros of the partial sums $\zeta_n(s) = \sum_{k=1}^n k^{-s}$ for $3 \le n \le 26$.
   • Scale: The study computed over 2.5 million certified zeros.
   • Tolerance: All results were computed within a guaranteed tolerance of $9.9 \times 10^{-21}$.
   • Performance: Benchmarks were conducted on a 24-thread AMD Ryzen 9 system. For example, computing zeros for $n=26$ with imaginary parts up to $100,000$ took approximately 6,711 seconds. The paper notes that the parameter $\rho$ (controlling step reduction in the adaptive search) significantly impacts execution time, with optimal values varying by problem instance.

Significance and Claims
The paper claims to provide a robust and certified method for zero-finding that overcomes the limitations of purely local methods (like standard Newton iteration) and purely global methods that may be computationally expensive.

• Robustness: By combining global argument counting with local refinement, the method ensures that no zeros are missed and that the results are mathematically certified via the Kantorovich theorem.
• Efficiency: The adaptive nature of the argument counting and the parallel bisection strategy allow the method to scale to millions of zeros.
• Applicability: The method is particularly effective for exponential and Dirichlet polynomials, where bounds on derivatives can be explicitly derived, though the authors acknowledge that the requirement for derivative bounds is a necessary condition for the robustness of their specific approach.

The authors emphasize that their contribution lies in the synthesis of these four components (counting, subdivision, exclusion, and refinement) into a single, parallelizable framework that delivers high-precision, certified results for complex analytic functions.