Ternary search

At what angle should you throw a tennis ball for it to fly furthest? Without air any schoolchild knows the answer: 45°. With air there is no formula: the range is computed by simulating the flight step by step, and one such trial is thousands of steps. But the main thing is known: as the angle rises from zero, the range first grows and then falls, and the curve has a single peak. So you can search the way binary search does, throwing away pieces of the interval. Choose a projectile and a method and see how many throws it takes to find the angle to a hundredth of a degree.

A tennis ball thrown at 30 m/s from a height of one metre flies furthest at 38.02°, 40.93 metres. Ternary search finds this angle in 46 flights, the golden section in 21. A feathered shuttlecock, which is braked much more strongly, has to be thrown flatter still, at 24.39°. And even the shot is slightly hindered by air: at 42.01° it flies 13 centimetres shorter than it would in a vacuum at 42.09°. Exactly 45° suits none of them: the shot is released from a height of 2.1 metres, and the ball and the shuttlecock are braked by air.

The idea

A function on an interval is called unimodal if it strictly increases up to some point and strictly decreases after it. The range of a throw as a function of angle is exactly like that. Neither a formula nor a derivative is needed for the search, only the ability to evaluate the function at any point.

Take two points m₁ < m₂ inside the interval and compare the values. If the value at m₂ is larger, the peak cannot lie left of m₁: left of the peak the function increases, so the value at m₁ could be smaller than at m₂ only if both points are on the ascent. Throw away the piece left of m₁. If the value at m₁ is larger, throw away the piece right of m₂. The interval shrinks, and the peak stays inside.

Classic ternary search puts the points at the thirds of the interval and throws away a third for two evaluations. Hence the name: the interval is divided into three parts. In the scene this shows in the lower strip: two points, a red line under the discarded part and the remaining part highlighted.

Why the peak is not lost

Note the word “strictly”. It will matter when we come to plateaus: if the function is constant somewhere, equal values say nothing, and the theorem does not work.

Termination follows from the theorem too. Each step keeps 2/3 of the interval, after k steps the length is (2/3)^k·L, and reaching precision ε takes ⌈log_1.5(L/ε)⌉ steps. For an angle from 0 to 90° to a precision of 0.01° that is 23 steps, that is, 46 flights.

The golden section: one evaluation per step

Ternary search is wasteful. After a step one of the two points ends up inside the new interval and its value is already known, but at the next step its position is not at a third, and the value is thrown away. What if the points are placed so that the old point lands exactly where a new one is needed?

Let the interval be [0, 1] and the points stand symmetrically: a = 1 − r and b = r. Throw away the right piece: the new interval is [0, r], and in it the points must stand in the same places at the same scale, that is, at r·(1 − r) and r·r. The old point a = 1 − r will do as the new right one if

This is the golden ratio, hence the method’s name: golden-section search. Each step keeps 0.618 of the interval but costs one evaluation, not two. Switch the scene to “Golden section”: the green circle marks the point whose value was taken from the previous step.

The fairest comparison of the methods is by how much the interval shrinks per function evaluation:

To shrink the interval tenfold, ternary search needs 11.4 evaluations, dichotomy 6.6, the golden section 4.8. The 2.4-fold difference depends neither on the function nor on the precision: it is in everything that takes time. There is an even more economical variant, Fibonacci search: if the number of evaluations is known in advance, the points are placed in the ratios of adjacent Fibonacci numbers, and the last step comes out slightly better. The golden section is its limit when there are many evaluations.

Code

The listing searches for the angle for the shot, the ball and the shuttlecock by three methods, counts the flights, and checks itself against the formula for a throw in a vacuum: tan θ = v/√(v² + 2gh). The flight is the Runge–Kutta method with a 2 ms step and quadratic air drag. Each evaluation function is wrapped in a counter, so the numbers in the output are honest numbers of flights.

Ternary search has one convenience: it is hard to write wrongly. The golden section has more variables, and it is easy to mix up which point goes where. But if evaluating the function is expensive, as here, the extra 25 flights pay for any amount of care.

If the slope is known

Two points are compared to find out whether the function is rising or falling. If the slope can be obtained otherwise, the peak is simply where it changes sign, and an ordinary binary search on the sign of the derivative finds it: one evaluation per step, the interval halved. This beats the golden section if the derivative costs no more than the function itself.

The most common case is a function on integers. Instead of the derivative, compare f(m) and f(m + 1): if the next value is larger, the peak is to the right of m. This is a binary search for the first m where f(m) ≥ f(m + 1), and it is shorter and more reliable than ternary search on integers with its ending “three points left, check them all”. The dichotomy in the listing is the same thing for real numbers, where a tiny shift plays the role of “plus one”.

If the derivative is smooth and can be computed, Newton’s method finds its zero fastest: the number of correct digits doubles at each step. That is how second-order optimisation methods work. Ternary search and the golden section are chosen when there is no derivative, when it lies, or when the function is a black box such as a simulation.

A peak can be found only halfway

Finding a peak has a limitation that finding a root does not. Near the peak the function is flat: f(x* + d) ≈ f(x*) − c·d². A double holds about sixteen significant digits, and as soon as c·d² falls below the last of them, the two points are indistinguishable. So the position of the peak is known only to a precision of about √ε, eight digits out of sixteen, however much you ask for.

   asked   golden: |x - 1/2|   bisection of f′
   1e-08             0.0e+00           3.7e-09
   1e-10             6.4e-10           2.9e-11
   1e-12             6.8e-10           4.5e-13
   1e-14             6.8e-10           3.6e-15

The golden section gets stuck at 6.8·10⁻¹⁰: beyond that the points are indistinguishable, and rounding decides the answer. Binary search on the sign of the derivative 1 − 2x gets to 3.6·10⁻¹⁵: near a root the function is not flat but steep. That is why library functions do not accept a precision below about √ε: the default tolerance of optimize in R is the fourth root of machine epsilon, about 1.2·10⁻⁴. The good news is in the same place: if the peak is flat, the value of the function at the point found hardly differs from the maximum. For the throw, an angle error of a thousandth of a degree costs less than a millimetre of range.

Where it breaks: noise, plateaus, several peaks

Everything said rests on one assumption: the function is unimodal. For the flight simulation that holds. For real data almost never. Here is the daily mean temperature in Almaty for every year since 1941. Summer is hot, winter is cold: there seems to be one peak. Try to find the hottest day of the year with ternary search.

In 1979 the hottest day was 3 August, 28.3 °C. Ternary search names 24 June, 22.2 °C. At the fourth step it compared two days in midsummer, the second turned out cooler, and the whole second half of summer, real heat included, went into the discard pile. After that the search honestly converges to the top of a small hill in June. Out of 85 years it is wrong in 53. Smoothing helps, but not completely: with weekly averages there are 36 misses, with monthly 13, with a two-month window still 10.

The cause lies not in the data but in the assumption. A day in midsummer can be colder than its neighbours because of rain, and then the curve has not one peak but dozens. Comparing two points answers the question “where is it higher at these two places”, not “where is the peak”. For finding a maximum in data there is a simple and reliable way: go through all of it. 365 comparisons take less than a millisecond.

Plateaus are a separate trouble. If the function is constant somewhere, the equality f(m₁) = f(m₂) says nothing about where the peak is, and no trick will help.

So on integers with repeats, say for a function that rounds time to the minute, ternary search may work or may not, depending on where the plateaus are. Strict unimodality or convexity is genuinely needed. Convex functions are the most reliable source of unimodal ones: a sum of convex functions is convex, and so is a maximum of convex functions. The article’s second problem is built on this.

And finally, there can be several peaks: in a neural network’s loss function, in the range of a throw with a spinning ball, in the cost of a route. Then ternary search finds one of the peaks, not necessarily the highest, just as gradient descent finds one of the pits. A grid scan with a search started in each cell helps.

What it costs

Time is the number of function evaluations: O(log(L/ε)), where L is the length of the interval and ε the required precision. The constant depends on the method: per evaluation ternary search shrinks the interval 1.22 times, the golden section 1.62 times. Memory is a few variables. On integers from 1 to n, O(log n) evaluations.

The measurement below is the time to find the angle for the tennis ball at different required precisions, from a tenth of a degree to 10⁻¹⁰. Every flight costs the same, so the lines run at a constant ratio: ternary search is about 2.3 times slower than the golden section at any precision.

Time grows linearly with the number of digits: each new digit costs the same number of flights. Note that beyond the eighth digit it grows for nothing: as shown above, the peak does not get any more precise.

Where it works

Optimisation libraries. In SciPy the function minimize_scalar can search by the golden section, but by default it takes Brent’s method: the golden section, sped up where possible by a parabola through three known points. The documentation says as much: Brent is usually better. The optimize function in R is the same Brent’s method, after his 1973 book “Algorithms for Minimization without Derivatives”.

Choosing a step in multidimensional optimisation. A method that has chosen a descent direction must decide how far to go. That is a search for the minimum of a function of one variable, the step length, and textbooks solve it with the golden section or Brent’s method.

Olympiad problems. The minimum of a convex function: the distance between moving points, the spread of a group, a cost made up of a growing and a shrinking part. If the function depends on two variables and is convex, the searches are nested: for each x an inner search finds the best y. This is legitimate, because the minimum of a convex function over one variable is a convex function of the other.

Where it does not work. Finding the maximum in an array of data, as in the temperature scene: that needs a full scan. And searching in a sorted array: for that there is binary search, while ternary search, sometimes proposed “because three is more than two”, makes more comparisons: two checks throw away only two thirds.

History

In 1953 the statistician Jack Kiefer published the note “Sequential minimax search for a maximum”. He was looking for a way to find the maximum of a unimodal function in a given number of evaluations so that in the worst case the remaining interval would be as short as possible, and he proved that the best way is Fibonacci search. The golden section comes out of it when there are many evaluations. In 1973 Richard Brent, in a book on optimisation without derivatives, combined the golden section with parabolic interpolation: the parabola quickly finds the peak of a smooth function, and the golden section is the safety net where the parabola goes wrong.

The throwing problem is older than computers themselves. That in a vacuum a projectile launched at 45° flies furthest has been known since Galileo. Air resistance, which lowers this angle, tormented ballistics experts for centuries, and one of the first tasks of electronic computers was precisely firing tables, the same step-by-step flights as in our scene. Real shot putters release the shot below the theoretical angle for another reason: the steeper the release, the lower the speed a person manages to give the shot. In 2001 Nick Linthorne showed on five athletes that each has their own optimum and that it coincides with the athlete’s habitual angle.

Problems

The first is the throw from the scene: the function is expensive, and a full scan of angles will not be in time. In the second the function must first be checked for unimodality: at first glance it could have any number of peaks.