@user21820 lmao that is an oof

Also came up with a cool bisection variant.

Previously I mentioned the idea of using the geometric mean in a numbering system similar to scientific notation which scales exponentially.

The big downside to it in a practical setting, most solutions tend to have an order of magnitude around zero. It's not every day that reaching $|x|\approx10^{\pm100}$ isn't just asymptotic behavior.

And so I came up with the idea of "log-log-scaling" things to accelerate convergence to smaller values while retaining both asymptotically optimal convergence (arithmetic mean) and fast convergence to the right order of magnitude (geometric mean).

The formula for $x,y>1$ is $\exp(\sqrt{ \ln(x)\ln(y)})$.

For example, with $x=10^{10}$ and $y=10^{1000}$, we have:

$AM(x,y)\approx 5\times10^{999}$

$GM(x,y)= 10^{505}$

$LM(x,y)=10^{100}$

I think I've also tried exponential searching in the past but in terms of the rate at which it converges to the right order of magnitude, it's only guaranteed to be about as the arithmetic mean, except it starts from the smaller side.

So such methods are not as good as far as I've seen.

In the context of double precision, my "log mean" tends to converge on an interval of roughly $(2,2^{1024})$

[cont.] in only say 4 or 5 iterations to the right order of magnitude, which is very satisfactory.

@SimplyBeautifulArt: I thought the optimal is simply to binary search on the exponent followed by the mantissa? This is the structure of the floating point type, so in some sense you don't have a choice if you want optimality.

@user21820 Yeah that's right, but practically speaking most problems do not have large exponents. Even if such a problem did involve large exponents, the suggested mean is still fast converging, especially considering that the double does not have very many exponents to search through.

@SimplyBeautifulArt What I'm asking is why you need to have a funny sort of mean (that is expensive to compute) instead of just using two binary searches.

@user21820 I mean whether or not it really is expensive in comparison depends on how expensive your function evaluations are in comparison. I also think that in a lot of problems, it tends to be more expensive to evaluate large inputs than smaller inputs, or it can be troublesome to deal with due to overflow.

@SimplyBeautifulArt Oh ok.

@SimplyBeautifulArt Is there a way to go around the x,y>1 condition?

The full formulation that I used was:

$GM(x,y)=0$ if $\operatorname{sign}(x) \ne\operatorname{sign}(y)$.

$LM(x,y)=0$ if $\operatorname{sign}(x) \ne\operatorname{sign}(y)$.

$LM(x,y)=\operatorname{sign}(x) \exp(GM(\ln|x|,\ln|y|))$ if $\operatorname{sign}(x) =\operatorname{sign}(y)$.

Ah so you did handle the sign problem.

As far as the bisection goes, it also involves swapping over to $AM$ once $x\in[y/4,4y]$, and it also involves minimum tolerance steps so that you don't just get stuck at something like $x=1$ on every step but instead move over to something like $1\pm10^{-13}$ or something.

Mm...

Feels a bit too ad-hoc haha..

But I guess that's the way floating point is.

Though your problem with getting stuck at 1 is not due to floating point structure, so maybe there is a better solution.

Yeah, but I haven't really much better ideas.

@SimplyBeautifulArt The curious thing is that taking AM on the log-scale doesn't have this issue. But you want something that skews towards small exponents...

I know. Why not just use that ( x,y ↦ exp(AM(ln(x),ln(y))) ) in the second stage? In the first stage, eliminate 0, then split into positive and negative cases, in each case doing a doubly-exponential search phase to get an interval that works. That is, try [2^(−2^k),2^(2^k)] for k from 1 upwards until finding one that works, and then apply AM on the log-scale.

@SimplyBeautifulArt: That way, you avoid large exponents in the first place, but still have optimality.

@user21820 That's just GM lol

Idk how I feel about the doubly-exponential search, it still feels a little wonky to me.

Worth noting that aside from quickly converging, bracketing both sides of the root is fairly important, which is the main reason I'm not so hot on the idea of searching outwards from the smaller side.

@SimplyBeautifulArt Oh lol yes it is, good observation! =D

Some methods, like Chandrupatla's method, will be able to decide when interpolation search seems feasible on their own, but they tend to require both sides of the bracket to be somewhat reasonable.

@SimplyBeautifulArt You don't solve the problem of bracketing in your solution either; you are starting from some humongous interval, right? That's actually conceptually inferior to having an initial search phase because it's not actually generic.

Well it pulls in values from the opposite side of the bracket a lot faster, usually in only 2 to 5 iterations I'd say. Even if the interval is still large, interpolation search might still be applicable.

Say for example our bracket is $(2,2^{1024})$ so that $LM(\dots)=2^32$. If this is insufficient, the next iteration or 2 should already give a sufficient bracket.

@user21820 I think this also has issues with solutions very close to $0$.

@SimplyBeautifulArt I don't think so because 2^(−2^k) very quickly cuts underneath (if the root is positive).

Hmm well it is only 10 iterations for double precision.

Still think mine tends to be faster though. And a bit more natural in the sense of bracketing.

@SimplyBeautifulArt That's not the point. I'm saying that you start with an ad-hoc initial bracket that is not conceptually correct, whereas an initial doubly-exponential search phase is conceptually watertight and additionally avoids even computing for unnecessary large points.

See, why even compute for 2^(2^10) if the solution isn't anywhere near that?

@user21820 The intended usage is that one may use $(-\infty,\infty)$ and similar as initial brackets, which gets "rounded down". Usually these points aren't computed but rather provided, either as asymptotic real value or $\pm\infty$. There isn't a clear guarantee that I see with your approach that any bracket will be found though, without such an assumption beforehand.

@SimplyBeautifulArt I don't understand why you think my approach can fail to find a bracket! It covers the entire real line!

Well I mean I can try implementing yours anyways and see how that'll fit in with everything else.

Maybe my description was not clear enough?

@user21820 I meant you still have to check the endpoints. You may not find any sign changes during your steps, and there's no other way to determine which side (towards $0$ or $\pm\infty$) should the root lie in. Besides, verification of a valid bracket is necessary for many methods (particularly those involving interpolation needing the actual values).

@SimplyBeautifulArt Yes, but the initial search only checks up to the smallest 2^(2^k) bigger than the true root, whereas yours just starts from something big regardless of where the root lies.

As a fun aside, $(\exp(x+20)-2) (1+10^{-10} + \sin(\ln(1+x^2)))$ with an initial bracket of $(-a,0)$ is quite a nasty test case for many bracketing method implementations. It tricks stuff like Brent's method into trying interpolation search repeatedly.

@user21820 I mean it depends on the provided bracket's size.

Honestly also not sure about how to fit your method into other brackets, like say $(5, 10^{100})$.

@SimplyBeautifulArt Are you saying you require a provided bracket? Mine doesn't.

Oh I was dumb. Mine totally fails.

For some reason, I assumed that the function was monotonic...

Am I correct that you require a bracket [a,b] such that f(a) < 0 < f(b)?

Yeah

Ok, and you expect the root to have small exponent regardless of the given bracket, right?