[Paper Review] A recurrence scheme for least-square optimized polynomials

C. Gebert, I. Montvay|ArXiv.org|Feb 28, 2003

Particle physics theoretical and experimental studies2 references3 citations

TL;DR

This paper presents a recurrence scheme for computing high-degree least-square optimized polynomials used in fermionic lattice field theory simulations, particularly within the two-step multi-boson (TSMB) algorithm. By leveraging orthogonal polynomial expansions and high-precision arithmetic, the method enables accurate approximation of inverse powers near zero—critical for light fermion masses—while allowing efficient evaluation using single-precision arithmetic, with a C implementation supporting polynomial degrees up to several thousand.

ABSTRACT

A recurrence scheme is defined for the numerical determination of high degree polynomial approximations to functions as, for instance, inverse powers near zero. As an example, polynomials needed in the two-step multi-boson (TSMB) algorithm for fermion simulations are considered. For the polynomials needed in TSMB a code in C is provided which is easily applicable to polynomial degrees of several thousands.

Motivation & Objective

• To develop a numerically stable and efficient method for computing high-degree polynomial approximations to inverse power functions near zero, which are essential in fermionic lattice field theory simulations.
• To address the ill-conditioning of the least-square optimization problem for functions like $x^{-\alpha}\bar{P}(x)$ over intervals $[\epsilon, \lambda]$ with $\epsilon \ll \lambda$, where $\lambda/\epsilon$ can exceed $10^6$.
• To provide a practical C implementation that computes these polynomials with minimal computational cost, even for degrees in the thousands, while ensuring high accuracy.
• To enable accurate evaluation of the resulting polynomials using standard 32-bit floating-point arithmetic through orthogonal expansion techniques, avoiding the need for high-precision arithmetic in runtime.

Proposed method

• The method uses orthogonal polynomials $\Phi_\mu(x)$ defined via a three-term recurrence relation with coefficients $\beta_\mu = -p_\mu / q_\mu$ and $\gamma_{\mu-1} = -q_\mu / q_{\mu-1}$, where $p_\mu$, $q_\mu$, and $s_\mu$ are integrals over the interval $[\epsilon, \lambda]$ with a weight function $w(x)$.
• The expansion coefficients $d_{n\nu} = b_\nu / q_\nu$ are computed from inner products $b_\nu = \int_\epsilon^\lambda w(x)^2 f(x) \Phi_\nu(x) dx$, minimizing the least-square deviation norm $\delta_n^{\text{min}}$.
• High-precision arithmetic is used during the recurrence computation to ensure numerical stability in the face of ill-conditioned matrices arising from large $\lambda/\epsilon$ ratios.
• Polynomial evaluation is performed via the orthogonal expansion $P_n(x) = \sum_{\nu=0}^n d_{n\nu} \Phi_\nu(x)$, allowing accurate evaluation even with single-precision floating-point arithmetic.
• For the TSMB algorithm, the method computes three main polynomials: $P^{(1)}$ approximating $x^{-\alpha}$, $P^{(2)}$ approximating $P^{(1)}$ with a polynomial weight, and $P^{(3)}$ approximating the inverse square root of $P^{(2)}$, with iterative refinement via Newton-type iteration.
• A C code is provided that supports multiple polynomial types, including iterative refinement of $P^{(3)}$ using Neuberger's formula or a Newton iteration starting from a crude initial approximation, with adjustable interval bounds and weight functions.

Experimental results

Research questions

• RQ1 How can high-degree least-square optimized polynomials be computed stably for functions with large condition numbers, such as $x^{-\alpha}$ near zero, when $\epsilon \sim 10^{-6}$ and $\lambda \sim 1$?
• RQ2 What recurrence scheme enables accurate computation of such polynomials using high-precision arithmetic while allowing efficient evaluation with standard 32-bit floating-point operations?
• RQ3 How can the approximation of $P^{(3)} \approx (P^{(2)})^{-1/2}$ be improved, especially at the lower end of the interval, to reduce uncorrectable errors in the TSMB algorithm?
• RQ4 What iterative method can efficiently refine the polynomial approximation of $P^{(3)}$ without relying on expensive Neuberger’s formula for large $K$?
• RQ5 Can a low-order initial guess for $P^{(3)}$ be used effectively in a Newton-type iteration to achieve high-precision results with minimal computational cost?

Key findings

• The recurrence scheme enables stable computation of least-square optimized polynomials up to degree several thousand, even for condition numbers exceeding $10^6$, by using high-precision arithmetic during the recurrence phase.
• The use of orthogonal polynomial expansions allows accurate evaluation of the resulting polynomials using only 32-bit floating-point arithmetic, significantly reducing runtime cost.
• For the third polynomial $P^{(3)}$, a Newton-type iteration starting from a crude initial approximation (e.g., $P^{(2)'}$) achieves high-precision results in just two to three iterations, making Neuberger’s formula unnecessary for high-accuracy cases.
• The C implementation supports multiple polynomial types and configurations, with output files storing recurrence coefficients and polynomial parameters for $P^{(1)}$, $P^{(2)}$, and $P^{(3)}$ separately.
• The required number of digits for high-precision computation scales linearly with polynomial degree: approximately $40 + 1.6n_1$ for the first polynomial, and $70 + 2.2n_{2,4}$ for the second and fourth, reflecting increasing complexity due to nested integrals.
• The method achieves machine-precision accuracy for the approximations, with the final $P^{(5)}$ sequence (iterative refinement of $P^{(3)}$) showing convergence independent of the initial guess after a few iterations, validating the robustness of the Newton iteration approach.