Remez and function approximations

I’ve been learning recently about how to approximate functions by low-degree polynomials.

This is useful in fully homomorphic encryption (FHE) in the context of “arithmetic FHE” (see my FHE overview article), where the computational model makes low-degree polynomials cheap to evaluate and non-polynomial functions expensive or impossible.

In browsing the state of the art I came across two interesting things. The first is the software package lolremez that implements polynomial (and rational polynomial $f(x) / g(x)$) function approximation using the so-called Remez algorithm. After building it it seems to work rather well. The author also wrote a Wiki on the GitHub project with a bunch of tips for using the API to do special things like round the coefficients to easily-representable numbers, or enforce odd-function symmetry in the approximation.

Lolremez outputs the approximated polynomial’s coefficients in Horner evaluation format. Horner’s method is bad for me, or at least non-optimal, because it does not optimize for the number of multiplication operations, which is a bottleneck in FHE. Instead, there’s a nice technique from the computational matrix algebra world called Paterson-Stockmeyer that reduces the number of (non-scalar) multiplications from $O(n)$ to $O(\sqrt{n})$. This is what most FHE people use by default, as far as I can tell.

So I wrote a little Python CLI to convert lolremez output to coefficient form. It uses python-fire and sympy to do the hard work.

import fire
import sympy


def convert_lolremez(
    name: str, poly_str: str, use_odd_trick: bool = True
) -> str:
  """Convert a lolremez output string to coefficient form.

  Example poly_str input:
    x*((((-2.6e+2*x**2+7.4e+2)*x**2-7.9e+2)*x**2+3.8e+2)*x**2-8.6e+1)*x**2+8.8

  Args:
    name: the name of the function being approximated
    poly_str: the polynomial expression to evaluate
    use_odd_trick: if True, map x -> x^2 and multiply whole poly by x

  Returns:
    The C++ header representation of the coefficient form of the polynomial.
  """
  x = sympy.symbols("x")
  expr = eval(poly_str)

  if use_odd_trick:
      expr = expr.replace(x, x**2)
      expr = x * expr

  poly = expr.as_poly()
  # change to coefficient list in degree-increasing order, with zeroes included
  coeffs = list(reversed(poly.as_list()))
  comma_sep_body = ",\n".join([f"  {v:.16f}" for v in coeffs])

  return f"""
// Polynomial approximation for {name}
//
//   {poly}
//
static constexpr double {name.upper()}_APPROX_COEFFICIENTS[{len(coeffs)}] = {{
  {comma_sep_body}
}};
"""


if __name__ == "__main__":
  fire.Fire(convert_lolremez)

The second interesting thing is how FHE researchers have adapted Remez to FHE specifically. This paper of Joon-Woo Lee, Eunsang Lee, Yongwoo Lee, Young-Sik Kim, and Jong-Seon No develops a “multi-interval” Remez method that jointly optimizes over multiple disjoint intervals of approximation. This is useful for approximating something like a sign function that has a discontinuity at zero (though you can do it in lolremez using odd-function symmetry tricks). The authors use the multi-interval method with $[-1, -\varepsilon] \cup [\varepsilon, 1]$ for small $\varepsilon > 0$ to approximate sign. There is an implementation of this algorithm in the lattigo library, which is also the project that tipped me off about Paterson-Stockmeyer. Kudos to Jean-Philippe Bossuat, its primary author, for making so much useful knowledge public and easy to find.