samhocevar / lolremez

📈 Polynomial Approximations using the Remez Algorithm
Do What The F*ck You Want To Public License
396 stars 36 forks source link

Multi-interval remez algorithm #28

Open ramym1 opened 5 months ago

ramym1 commented 5 months ago

Hi, this project is very impressive, thanks for sharing this. I wanted to ask if this implementation supports the multi-interval Remez algorithm. For example, is it possible to approximate sing(x) function on the intervals [-1,-0.01], [0.01,1]?

Thanks in advance!

aadler commented 5 months ago

Why don't you just find separate minimax interpolating polynomials for each interval?

ramym1 commented 5 months ago

I want to approximate a polynomial for sign function in [-1,1]. My usecase is that I have an encrypted input x \in [-1,1] but I can't know if x is negative or positive. I am using "fully homomorphic encryption" which supports evaluating polynomials on encrypted data.

samhocevar commented 5 months ago

That’s a very interesting use case! I’ll have a look at what can be done.

ramym1 commented 5 months ago

Thanks a lot! actually this problem is solved in the state of the art by what is called "multi-interval Remez algorithm", which computes Remez on [-1,-eps], [eps, 1]: https://eprint.iacr.org/2020/834.pdf

samhocevar commented 5 months ago

So, this specific case is actually a lot less complex than the general case, because you can use the technique for odd functions explained in the lolremez manual, the function you want to approximate is simply $f(x) = 1$ and the range is $[\sqrt{0.01},1]$.

For instance:

$ lolremez -d 5 -r 0.1:1 '1/sqrt(x)' '1/sqrt(x)'
double f(double x)
{
    double u = -27.674387189032068;
    u = u * x + 89.61852558307946;
    u = u * x + -113.12718846240877;
    u = u * x + 70.763443739725503;
    u = u * x + -23.464237705608777;
    return u * x + 4.8738244889926898;
}
$

Replacing $x$ with $x^2$ and multiplying the final result by $x$:

image

ramym1 commented 5 months ago

Thank you! do you have an idea on how to solve the general case?

ramym1 commented 5 months ago

Or more specifically, is it possible to compute an approximation of a stepwise function? For example, approximating the following function: f(x) = {-3 if x < -1, 0 if -1 < x < 1, 2 if 1 < x < 3} In the domains: [-3+eps, -1-eps], [-1+eps, 1-eps],[1+eps, 3-eps]`.

One way to do this would be to express f(x) as the sum of two stepwise functions that have two branches, and then using your solution to approximate a stepwise function with two branches. However, I wonder if there is a way to compute an approximation of f(x) using Remez directly.

samhocevar commented 5 months ago

The Remez algorithm doesn’t work well with non-continuous functions. I would recommend helping it by making the functions continuous, or even smooth (C-infinity). You can for instance build a good such function using tanh(). In your case:

$$-3 + 3\frac{1+\tanh k(x+1)}{2} + 2\frac{1+\tanh k(x-1)}{2}$$

The greater $k$ is, the closer to the stepwise function you will get. For instance with $k = 5$ and $k = 30$:

image

Using the Remez algorithm on $[-3,3]$ will then give something like that:

image

However there is clearly a problem with LolRemez and that kind of function; most of the time it fails to converge and find a solution. I am going to investigate.

ramym1 commented 5 months ago

Great, thanks!

yellow123Nike commented 2 months ago

So, this specific case is actually a lot less complex than the general case, because you can use the technique for odd functions explained in the lolremez manual, the function you want to approximate is simply f(x)=1 and the range is [0.01,1].

For instance:

$ lolremez -d 5 -r 0.1:1 '1/sqrt(x)' '1/sqrt(x)'
double f(double x)
{
    double u = -27.674387189032068;
    u = u * x + 89.61852558307946;
    u = u * x + -113.12718846240877;
    u = u * x + 70.763443739725503;
    u = u * x + -23.464237705608777;
    return u * x + 4.8738244889926898;
}
$

Replacing x with x2 and multiplying the final result by x:

image

"How to approximate only use 7th-order polynomial?"

j2kun commented 2 months ago

Funny that I also saw this in search of a remez-like algorithm for FHE. Note that a Go implementation of the multi-interval method exists in Lattigo: https://github.com/tuneinsight/lattigo/blob/0d754047425ff9513bf6fbea717415a754911408/he/hefloat/minimax_composite_polynomial.go#L125