# Approximate implicit function on real line

**URL:** <https://discourse.julialang.org/t/approximate-implicit-function-on-real-line/125421>\
**Category:** Numerics\
**Tags:** question, approximation\
**Created:** [January 31, 2025, 8:19am UTC](https://discourse.julialang.org/t/approximate-implicit-function-on-real-line/125421 "2025-01-31T08:19:23Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [January 31, 2025, 8:19am UTC](https://discourse.julialang.org/t/approximate-implicit-function-on-real-line/125421/1 "2025-01-31T08:19:23Z")

</div>

I have a family of functions y(x) defined implicitly by

x = \log\left(\sum\_{i=1}^N \exp(a\_i + y b\_i) \right)

where 0 \< b\_1 \< \dots \< b\_N can be assumed.

It is easy to find the solution numerically (eg with Newton’s method). But I need to evaluate it a gazillion times for a given set of a\_i, b\_i, so I want to speed up the calculation. Currently it is taking up 70% of my runtime 🙁

It can be shown (I will spare you the algebra) that y(x) is increasing and strictly concave. It also converges to asymptotes on both ends, which can be characterized using a\_1, b\_1, a\_N, b\_N in a simple way. Here is a mockup illustration:

 ![asymp](https://global.discourse-cdn.com/julialang/original/3X/e/e/eed1848a1232ce276ea1c5caf9d5b484b20d34f7.png)

The idea is that for a particular set of a's and b's, I would solve for y(x) on some grid, approximate that, and from then on look up values from a quick approximation. The approximation should be AD friendly but that should not be difficult.

If I had \lim\_{x \to \pm\infty} y(x) = 0 or some constant I would know how to approximate it with Chebyshev polynomials (appropriately transformed to \mathbb{R}). But the asymptotes are giving me a hard time.

It would be great to have a nice (= cheap, smooth, …) function g(x) that fits the asymptotes, then I would approximate y(x) - g(x).

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [January 31, 2025, 10:34am UTC](https://discourse.julialang.org/t/approximate-implicit-function-on-real-line/125421/2 "2025-01-31T10:34:49Z")

</div>

I am leaning to approximating the intercepts with a modified `squareplus`. Specifically, consider

f(x; A, b, C) = \frac{C}{2}\left(x + A/C + \sqrt{(x + A/C)^2+b^2}\right) 

where C \ne 0. It is easy to show that \lim\_{x\to-\infty} f(x) = 0 (just like `squareplus`) and \lim\_{x\to\infty} f(x) - Cx - A = 0 (generalizing the slope and intercept).

Then, hoping that I did the algebra right,

-f(-x; A\_1, b, C\_1) + f(x; A\_2, b, C\_2)

will have intercept and slope A\_1, C\_1 at -\infty and A\_2, C\_2 at \infty. I just need to think may way thought the numerics to see if there are any problems, but I don’t think I can find anything cheaper than this.

(I also considerd rotating a hyperbola to have two oblique asymptotes but I kind of gave up on the algebra of doing that. 😉)

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [January 31, 2025, 12:29pm UTC](https://discourse.julialang.org/t/approximate-implicit-function-on-real-line/125421/3 "2025-01-31T12:29:45Z")

</div>

ApproxFun and `solve`?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 31, 2025, 1:10pm UTC](https://discourse.julialang.org/t/approximate-implicit-function-on-real-line/125421/4 "2025-01-31T13:10:11Z")

</div>

> [@Tamas\_Papp](#):
>
> It is easy to find the solution numerically (eg with Newton’s method). But I need to evaluate it a gazillion times for a given set of a\_i, b\_i, so I want to speed up the calculation. Currently it is taking up 70% of my runtime 🙁 […]  
> I would solve for y(x) on some grid, approximate that, and from then on look up values from a quick approximation.

See also: [Approximate inverse of a definite integral - #3 by stevengj](https://discourse.julialang.org/t/approximate-inverse-of-a-definite-integral/83181/3)

For a smooth function, you can get exponential convergence on a finite interval by evaluating on a Chebyshev grid and then interpolating with a polynomial, ala ApproxFun.jl or FastChebInterp.jl (as in the above post).

But I guess you have already thought about Chebyshev polynomials?

> [@Tamas\_Papp](#):
>
> If I had lim x→±∞ y(x)=0 or some constant I would know how to approximate it with Chebyshev polynomials (appropriately transformed to R) But the asymptotes are giving me a hard time.

Just subtract off the asymptotes and approximate the what’s left? Oh, I see that you already thought of this:

> [@Tamas\_Papp](#):
>
> It would be great to have a nice (= cheap, smooth, …) function g(x) that fits the asymptotes, then I would approximate y(x) - g(x).

> [@Tamas\_Papp](#):
>
> The approximation should be AD friendly but that should not be difficult.

FastChebInterp.jl should be AD-friendly (it works with ForwardDiff, and it also hooks into ChainRules.jl to expose its own optimized derivatives to AD).

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [January 31, 2025, 1:30pm UTC](https://discourse.julialang.org/t/approximate-implicit-function-on-real-line/125421/5 "2025-01-31T13:30:47Z")

</div>

> [@stevengj](#):
>
> get exponential convergence on a finite interval

My problem is that I do not know the region I will be evaluating at in advance so I need the **whole real line**. Does ApproxFun.jl / FastChebInterp.jl allow infinite domains? I am only familiar with the latter and AFAIK it works on hypercubes.

(If I had a finite interval the problem would be trivial.)

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 31, 2025, 1:47pm UTC](https://discourse.julialang.org/t/approximate-implicit-function-on-real-line/125421/6 "2025-01-31T13:47:22Z")

</div>

> [@Tamas\_Papp](#):
>
> My problem is that I do not know the region I will be evaluating at in advance so I need the **whole real line**. Does [ApproxFun.jl](https://juliaregistries.github.io/General/packages/redirect_to_repo/ApproxFun) / [FastChebInterp.jl](https://juliaregistries.github.io/General/packages/redirect_to_repo/FastChebInterp) allow infinite domains? I am only familiar with the latter and AFAIK it works on hypercubes.

If your function goes exponentially quickly to zero (once you have subtracted off the asymptotics), as I think is the case here, you can can perform a change of variables to map it to a finite interval and it will still converge quickly.

For example, if you let x = t / (1-t^2), mapping t \in (-1,1) to x \in (-\infty, +\infty), then for something that decays exponentially fast you are mapping it to a [kind of bump function](https://en.wikipedia.org/wiki/Bump_function) with [essential singularities](https://en.wikipedia.org/wiki/Essential_singularity) at the endpoints, and Chebyshev approximation (similar to Fourier series/transforms) should still converge [faster than polynomially](https://arxiv.org/abs/1508.04376).

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 31, 2025, 6:52pm UTC](https://discourse.julialang.org/t/approximate-implicit-function-on-real-line/125421/7 "2025-01-31T18:52:26Z")

</div>

> [@stevengj](#):
>
> If your function goes exponentially quickly to zero (once you have subtracted off the asymptotics), as I think is the case here

> [@Tamas\_Papp](#):
>
> I am leaning to approximating the intercepts with a modified `squareplus`.

Note that your suggested asymptotic approximation is not good enough here, because it only approaches the asymptotic lines as O(b/x^2), not exponentially.

A simple exponential approximation to the asymptotes is:

```julia
function y_approx(x, a, b)
    f = (tanh(x) + 1)/2
    y1 = (x - a[begin]) / b[begin]
    y2 = (x - a[end]) / b[end]
    return y1 * (1 - f) + y2 * f
end

```

which converges exponentially fast to the exact y(x) for large |x|. However, you should be be able to do even better (get a steeper slope of exponential) by examining your x(y) function more closely and pulling out the next-order term.

With this approximation, I checked (on your function with random a and b) that x = t/(1-t^2) does indeed give a Chebyshev polynomial approximation on t = (-1,1) that converges faster than polynomially (the -\log |\text{error}| is proportional to \sqrt{\text{degree}}), and gives you 3 digits of accuracy with about degree 25. If you improve `y_approx` you should be able to do even better.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [February 3, 2025, 2:31pm UTC](https://discourse.julialang.org/t/approximate-implicit-function-on-real-line/125421/8 "2025-02-03T14:31:38Z")

</div>

Thank you very much for the concrete suggestion. I found that it is not (generally) concave, so I looked for another form, and found that

g(x) = (p - n) \log( \exp(kx) + 1) + knx

works to approximate asymptotes

\min(px, nx)

with p \< n, crossing at 0. This has

- g(x) \> 0,
- g''(x) \< 0,
- g'''(0) = 0,
- \lim\_{x \to \infty} g(x) - px = 0,
- \lim\_{x \to -\infty} g(x) - nx = 0,
- g'(0) can be parametrized to match c'(b) at some point.

These I can shift on the x and y axis to match the original function, and pick a k so that I match the derivatives at the crossing point. These provide a good approximation for my purposes with 10-12 Chebyshev polynomials.

So a costly Newton step is eliminated which was eating up 70% of my runtime. The replacement uses \<1% of the _reduced_ runtime. Took me about half a day to code with a lot of unit tests. Julia is amazing for a smooth transition from prototype to heavily optimized code.

Plot of g with p = 1/2, n = 2, k = 1:

 ![plot](https://global.discourse-cdn.com/julialang/original/3X/e/0/e0751d3fa9b41ac302f734752f65aa42174bfd5c.png)