An update on Herbie, Summer 2022

Rival is Herbie's interval-arithmetic library. It is, I think, a pretty high-quality implementation:

• It supports arbitrary precision, which Herbie uses for computable reals.
• It has had some attention put toward performance; for example, multiplication is quite a bit faster than in a "textbook" implementation.
• It supports many many functions, and not just the easy ones like sqrt or log; Rival has implementations of fmod, for example, and of pow for negative inputs.
• It has clear correctness criteria, and is tested extensively against them
• It has rigorous treatment of invalid inputs, like sqrt([-1, 1]), via what we call error intervals.
• It also supports arbitrary boolean operations, like comparing intervals, and produces boolean intervals as a result.
• It also supports "movability flags", which can detect lots of cases where higher-precision re-evaluation is fruitless.

Herbie's been using Rival for many releases now, and it ensures a level of soundness that I've very proud of.

So over the last few weeks I've done a bunch of minor maintenance and cleanup in Rival. First, I wrote some documentation, in case there are any other potential users of the library. That also lead me to think a little bit about the API we want to present, which is now cleaner. And finally, I added support for the final two functions, lgamma and tgamma, which had been missing in Rival. Amusingly, I talked about an interval version of these two functions in my last update post in February 2020.

Let me say just a bit more about that last point. The two Gamma functions are the most mathematically complex functions supported by Rival, and took quite a bit of work. The reason is basically this: computing the interval version of some function basically means solving an optimization problem, finding the function's minimum and maximum over a range. For all of the other functions supported by Rival, those minima and maxima are directly computable. However, there's no analytical closed form (as far as I know) for either the x or y positions of the Gamma function's minima and maxima. Therefore, finding those minima and maxima involves basically a binary search (actually, a ternary search) and also some correction for rounding error.

This all makes the interval gamma functions dreadfully slow—about 200ms on my computer for a single evaluation at 80 bits of precision. At higher precision, it's slower—much slower. (I tried running it at 8000 bits of precision and it never completed). But that's OK, because frankly I don't expect it to see much use. However, having these implementations around (even if they take a long time or time out for many unusual uses) means we can now enforce the invariant that all of Herbie's supported functions have interval versions, which in turn allows us to deprecate a bunch of fallback algorithms that were basically never used.

Amusingly, working on this also lead me to discover a bug I introduced into math/bigfloat, which I've now fixed.

A few years ago, I implemented an interactive profile explorer tool for Herbie. You can see it live on the Herbie nightly reports by scrolling to the bottom of the page. Basically, it takes the standard gprof output and hyperlinks it, so that instead of scrolling up and down looking for a function by number, you can click links to browse the profile. That's made it a lot easier to find inefficiencies in Herbie.

One such inefficiency I found recently was in Herbie's regimes pass. Regimes basically takes several different ways ("alts") to evaluate the same mathematical expression, and tries to learn a simple rule for which way is most accurate on which input. This involves a lot of repetitive trial-and-error, but we haven't put a lot of effort into speeding this up because it only makes up about 20% of Herbie runtime.

But when I was browsing around a Herbie profile, I noticed something odd. Of that 20%, about half was spent evaluating those mathematical expressions over and over again. This is strange, because in theory you should only need to evaluate each expression once, which shouldn't take long. When I checked the code, I saw the reason: Herbie was evaluating the expressions once during every search step, not once total. This was a super-easy fix (just hoist the evaluation out of the search loop) and lead to a big speed-up to regimes. Along the way I also upgraded regimes to use Herbie's batched evaluation function, which should produce another small speed-up.

A related component of Herbie is bsearch, which runs right after regimes and cleans up its output a bit. I simplified it, but I'm going to need to explain a bit about bsearch before the simplification makes sense.

regimes outputs a couple of things: an expression (e), a set of split points (v1, v2, …), and a list of alts (a1, a2, …) to use between each split point. Basically, this data corresponds to the following program:

(if (< e v1)
  a1
  (if (< e v2)
    a2
    ...))))

Moreover, each of the split points actually corresponds to an interval [v1a, v1b], where regimes has no way of picking any particular point in that interval.

The job of bsearch is to refine that interval using binary search. To do so, it considers the interval [v1a, v1b] and the two alts a1 and a2 used on either side of that interval. It then samples a bunch of inputs where e is equal to the midway points (v1a + v1b) / 2, evaluates a1 and a2 on them, determines which is more accurate, and uses that to refine the interval further.

The tough part is the part where it samples inputs where e is equal to a particular value. How? Well, here's where a trick comes in. Suppose the programs a1 and a2 both have a subexpression equal to e. And furthermore suppose that whatever variables are used in that subexpression, they're not used anywhere else in a1 or a2. Then we can treat e not as an expression, but as a symbolic variable, and then sampling points where e is any particular value is easy. This trick lets us run bsearch on many (though not all) results of regimes.

Anyway, the simplification I made to bsearch is that up until now, bsearch was sampling new points with e equal to a particular value. This makes sense, but unfortunately requires us to expose the sample-new-points function to the core of Herbie, and we've been trying for a while now to separate the two. If Herbie is separate from the sampling function, we would be able to run Herbie of points gathered by, say, dynamic analysis of some binary, or points generated by some sampling program other than Herbie. This resampling step in bsearch was the final thing standing in the way of complete separation.

So here's my simple fix that eliminates resampling. Instead of sampling new points, where the symbolic variable e is set to some specific value, just take the old points and add the variable e to them. This leads to basically the same results (hypothetically, there could be some over-fitting, but I think the chances and impact of that are small) and avoids the need to invoke the sampler!

For a long time, Herbie's series expansion algorithm has been a source of timeouts. Not only is this annoying (to users and to us developers), but it also limits us from calling the series expander as often as we'd like. Now, the series expander is pretty complicated—complicated enough that I've already written two posts about how it works, plus another post promising an e-graph version (it failed). So it might look like this is just a problem we're stuck with.

However, my suspicions where raised when Brett discovered that a simple cache in the series simplifier sped it up by a factor of two. Perhaps it wasn't series expansion, but the simplifier, that was the problem? Indeed, the simplifier took up something like 95% of series expansion runtime, and importantly, the simplifier is (frankly) crap. I wrote the simplifier way way back, in September 2014 or so, and at the time it was intended to be Herbie's main simplifier (this was before we discovered e-graphs). Therefore, it tries really hard to do a good job simplifying.

For series expansion, this is a waste. The only reason we have a special simplifier is because the series expander produces exponentially large terms if left to its own devices, and we can't feed those terms into our nice e-graph based simplifier for technical reasons. So the series simplifier should run fast and do very little; instead it seems to run forever.

Now, I had long totally forgotten how the series expansion simplifier works, but reading over the (mercifully short) code I noticed that it called cartesian-product, a sure sign of something bad happening. It turns out that the series expander would attempt to distribute terms in a multiplication, and to do so it would expand any product of sums into a sum of products. This is definitely super slow, and there's not really any reason to do it: even if this simplification were profitable, it would almost certainly be handled later by our full e-graph simplifier.

Turning this distribution handling off was very easy (just deleting a few lines of code) but it immediately reduced the time taken by series simplification, and series as a whole, to almost nothing, with no timeouts at all. And Herbie's results actually improved, presumably because we were now generating shorter expressions that later Herbie phases could handle better. This was an easy win and totally eliminates an annoying cause of timeouts.

One of Herbie's most important parameters is how many total iterations Herbie performs, controlled by the --num-iters flag. The default is four iterations, and I remember that way back when it used to be three, but we bumped it for reasons I don't remember. In the intervening time, though, Herbie has gotten dramatically faster, and in particular each iteration has gotten way, way faster (sampling, which happens before any of the iterations, has not been sped up as dramatically). So it makes sense to think about whether this default is still optimal.

To test this, I ran Herbie on its full test suite at a varying number of iterations and looked at its runtime and its accuracy. I expected accuracy to improve but runtime to worsen as the number of iterations crept up. Here are the results:

In this table, the first column shows the number of iterations, with the default (4) labeled in bold. The second column shows Herbie's "bits of error" metric, where higher means more accurate. Accuracy indeed increases steadily as the number of iterations goes up, but hits a clear knee at the default of four iterations. The first two and last one row are italicized, because with 0 or 1 iteration two fewer benchmarks timed out, while with 50 four more timed out; timed-out benchmarks don't contribute to the total. (Timeouts and other failures are shown in the fourth column.)

The third column shows total runtime, with the italicised rows run under less load and so artificially lower. Runtime increases pretty steadily throughout, with additional iterations adding a few minutes each. The last column shows how many benchmarks saturate and thus are not affected by additional iterations; there are 528 benchmarks in total, so with 50 iterations most benchmarks are saturated but there are still plenty that aren't. It looks like truly absurd iteration limits would be necessary to get everything to saturate.

All of the above suggests that our default of 4 is actually pretty much optimal, getting all of the large accuracy improvements without incurring the most extreme costs in running time. So I'm planning to leave the default in place. Of course, future changes to Herbie might recommend reducing the number—if each iteration gets better—or increasing it—if we find new things for additional iterations to do.