Successive parabolic interpolation is a well-known technique, but I rediscovered it today, working from the method of secants. I thought my derivation might make interesting reading.
The scenario is that you have a function whose derivative is not computable — perhaps it involves some kind of physical experiment — and you want to find its minimum.
If you could compute its derivative, you could find roots of the derivative function $f'$, whether algebraically or numerically. The method of secants is a numerical method for finding roots: you evaluate the derivative at a couple of guesses you think are close to the root, draw a line through the values at those guesses, find its intersection with the X-axis, and use that intersection as a new guess. This gives you a recurrence
$$g_n = g_{n-1} - f'(g_{n-1}) \frac{g_{n-1}-g_{n-2}} {f'(g_{n-1})-f'(g_{n-2})}$$which, under suitable circumstances, converges rapidly to a root of $f'$. And, as it happens, for functions that are expensive to evaluate, it converges more rapidly than Newton's method, requiring less than twice as many iterations as Newton's method, but only one new evaluation of $f'$ per iteration, while Newton's method requires evaluating both $f'$ and $f''$ in each iteration.
The root of the derivative $f'$ we find may or may not be a minimum of $f$, and if it is a minimum it may not be the global minimum, and moreover sometimes the algorithm fails to converge; but it works often enough to be worth trying.
Now, suppose we can't evaluate the derivative $f'$ at all; for whatever reason we are limited to evaluating $f$, the function itself. (This was common before the invention of automatic differentiation, but now it's rare.) In the same way that the method of secants is in some sense a cheap approximation of Newton's method by considering the secant line an approximation of the tangent line, we can treat two secants' slopes as the approximation of the derivative at their centers. Given three guesses at the minimum $g_{n-1},\ g_{n-2},\ g_{n-3}$ we can evaluate the function $f$ at those guesses and draw two secant lines; then we can construct a secant line between the thus-approximated values of $f'$ and find the root of that secant to get a new guess $g_n$ for the minimum of $f$.
Let $s_n$ be the slope of the secant line from $(g_{n-1}, f(g_{n-1}))$ to $(g_n, f(g_n))$:
$$s_n = \frac{f(g_n) - f(g_{n-1})}{g_n - g_{n-1}}$$Then our new guess at the location of the extremum can be calculated in this way, precisely analogous to the usual rule for the method of secants given above:
$$g_n = \frac 12 \left(g_{n-1} + g_{n-2} - s_{n-1} \frac{g_{n-1} - g_{n-3}}{s_{n-1} - s_{n-2}} \right)$$Now, it may be that this sequence diverges, or converges to a maximum or a "saddle point" where the second derivative is zero but the third derivative isn't. But it ought to sometimes close in on the true minimum rather fast.
Does it? Let's see.
import sympy, math
sympy.init_printing()
def secmin(f, g0, g1, g2):
g = [g0, g1, g2]
fs = [f(g0), f(g1), f(g2)]
s = [None]
while True:
while len(s) < len(g):
n = len(s)
s.append((fs[n] - fs[n-1])/(g[n] - g[n-1]))
if s[-1] == s[-2]:
return
gn = (g[-1] + g[-2] - s[-1] * (g[-1] - g[-3])/(s[-1] - s[-2]))/2
if gn == g[-1]:
return
g.append(gn)
fs.append(f(gn))
yield g[-1], fs[-1]
It should be able to minimize quadratics immediately from any starting point, unless it immediately blows up:
for a in [3, -3, 5]:
for b in [3, -1, -2, 6]:
print("##", a, b)
for gn, fn in secmin(lambda x: (x-a)*(x-b), 1, 2, 15):
print(gn, fn)
## 3 3 3.0 0.0 ## 3 -1 1.0 -4.0 ## 3 -2 0.5 -6.25 ## 3 6 4.5 -2.25 ## -3 3 0.0 -9.0 ## -3 -1 -2.0 -1.0 ## -3 -2 -2.5 -0.25 ## -3 6 1.5 -20.25 ## 5 3 4.0 -1.0 ## 5 -1 2.0 -9.0 ## 5 -2 1.5 -12.25 ## 5 6 5.5 -0.25
That seems to be working. (And it suggests another way to derive the same method: draw a parabola through the last three points and take the extremum of the parabola as your new guess. This may be a useful way to search for what this method is called in the literature.) How about a more interesting function?
k = lambda x: (x+1)*x*(x-1)
for gn, fn in secmin(k, .7, .8, .9):
print(gn, fn)
0.6062500000000002 -0.383429443359375 0.5963414634146341 -0.3842686390940352 0.5811552961647223 -0.384875047339243 0.5775551936438853 -0.3849001067153479 0.5773721869195749 -0.3849001786276855 0.5773505190716897 -0.3849001794596424 0.5773502705021032 -0.3849001794597505 0.5773502692057524 -0.3849001794597505 0.5773502698539278 -0.3849001794597505
.57735**2
Well, that does seem to have converged, to $\sqrt\frac 13$, which I think is the correct local minimum. Let's check.
x = sympy.symbols('x')
ks = k(x)
ks
ks.diff().simplify()
Yup, I guess $\sqrt\frac 13$ is right. We can just as easily find the local maximum, since this method doesn't distinguish between different ways the derivative can be zero:
for gn, fn in secmin(k, .7, -.8, -.9):
print(gn, fn)
-0.265 0.24639037500000002 -0.5522900763358779 0.3838281667089989 -0.5479324561598471 0.383426707716067 -0.5838522992425746 0.3848266797135901 -0.5774639810120341 0.3849001570622076 -0.5772935538340034 0.3849001738885636 -0.5773503739123432 0.38490017945973154 -0.5773502673298542 0.3849001794597505 -0.5773502691996591 0.38490017945975047 -0.5773502577930631 0.3849001794597503 -0.5773502650682871 0.38490017945975047 -0.5773502671339732 0.38490017945975047 -0.5773502661011302 0.38490017945975047
gn+(1/3)**(1/2)
Not very precise I guess. After 13 iterations it only has 9 decimal places right, though it had those after about 7 or 8 iterations.
We can also get a place where the derivative is zero that isn't an extremum, though this particular case has extremely slow convergence (just as for Newton's method):
for i, (gn, fn) in enumerate(secmin(lambda x: (x+1)**3, 1, 2, 3)):
print(i, gn, fn)
0 0.44444444444444464 3.01371742112483 1 0.3092105263157894 2.2440290015855076 2 -0.04453945232167994 0.8722445757650734 3 -0.3904085036790089 0.22652529265667123 4 -0.5422437632616212 0.09591859535888844 5 -0.6789480814212261 0.03309221282277016 6 -0.7761025270859593 0.011223997865658832 7 -0.8397651878779913 0.0041140600495487125 8 -0.887119864717784 0.0014383102111333276 9 -0.9202866260793142 0.0005065164728117844 10 -0.943516097226348 0.00018020800980435872 11 -0.960099775160477 6.352227284635844e-05 12 -0.9717860153086274 2.2459148105993188e-05 13 -0.9800434347299882 7.947991436617367e-06 14 -0.9858916308688045 2.8082155642968306e-06 15 -0.9900234128559526 9.92992575234892e-07 16 -0.9929452213022262 3.511156477678299e-07 17 -0.9950117203817831 1.2412303026650733e-07 18 -0.9964726991213056 4.388615368916267e-08 19 -0.9975058154407956 1.5516213933544667e-08 20 -0.9982363568945098 5.485700776359981e-09 21 -0.9987529113645043 1.9395097393716187e-09 22 -0.9991181755212277 6.857194229851004e-10 23 -0.9993764565046607 2.4243775814686385e-10 24 -0.9995590878341095 8.57148850096849e-11 25 -0.9996882280852705 3.030476848385423e-11 26 -0.9997795439855649 1.071435063725322e-11 27 -0.9998441140364727 3.788096509741164e-12 28 -0.9998897719847694 1.3392941217389454e-12 29 -0.9999220570230744 4.735119755431835e-13 30 -0.9999448859913433 1.6741177470668976e-13 31 -0.9999610285112626 5.918899819412864e-14 32 -0.9999724429959704 2.09264711577092e-14 33 -0.9999805142555317 7.398624887703915e-15 34 -0.9999862214979851 2.6158088947768817e-15 35 -0.9999902571277819 9.248281063944923e-16 36 -0.9999931107489849 3.2697611293785863e-16 37 -0.9999951285638919 1.1560351322817576e-16 38 -0.9999965553744931 4.08720140954964e-17 39 -0.9999975642819455 1.4450439161425947e-17 40 -0.9999982776872467 5.1090017604550546e-18 41 -0.9999987821409728 1.806304894931244e-18 42 -0.9999991388436233 6.386252200568818e-19 43 -0.9999993910704864 2.257881118664055e-19 44 -0.9999995694218117 7.982815250711023e-20 45 -0.9999996955352431 2.8223513998738135e-20 46 -0.9999997847109058 9.978519071107503e-21 47 -0.9999998477676216 3.527939245982905e-21 48 -0.999999892355453 1.2473148819587568e-21 49 -0.9999999238838109 4.409924047830225e-22 50 -0.9999999461777265 1.559143602448446e-22 51 -0.9999999619419054 5.512405059787782e-23 52 -0.9999999730888632 1.9489295030605575e-23 53 -0.9999999809709527 6.890506385037264e-24 54 -0.9999999865444316 2.436161908976965e-24 55 -0.9999999904854763 8.613133132052922e-25 56 -0.9999999932722158 3.0452023862212063e-25 57 -0.9999999952427382 1.0766416038175297e-25 58 -0.999999996636108 3.8065027943310834e-26 59 -0.999999997621369 1.345802098994627e-26 60 -0.999999998318054 4.758128492913854e-27 61 -0.9999999988106846 1.6822523881865023e-27 62 -0.9999999991590269 5.947661793926281e-28 63 -0.9999999994053423 2.102814896341257e-28 64 -0.9999999995795135 7.434577242407851e-29 65 -0.9999999997026712 2.6285186204265714e-29 66 -0.9999999997897567 9.293221553009814e-30 67 -0.9999999998513355 3.2856556366864105e-30 68 -0.9999999998948783 1.1616545344147633e-30 69 -0.9999999999256677 4.107078747313252e-31 70 -0.9999999999474392 1.4520681680184541e-31 71 -0.9999999999628338 5.133848434141565e-32 72 -0.9999999999737196 1.8150852100230677e-32 73 -0.999999999981417 6.417253033710556e-33 74 -0.9999999999868598 2.2688277581267398e-33 75 -0.9999999999907085 8.02142252101062e-34 76 -0.9999999999934299 2.8360346976584248e-34 77 -0.9999999999953543 1.0026778151263275e-34 78 -0.999999999996715 3.545043372073031e-35 79 -0.9999999999976771 1.2534371274731826e-35 80 -0.9999999999983574 4.431753515910376e-36 81 -0.9999999999988386 1.5665717682967035e-36 82 -0.9999999999991788 5.538568680800781e-37 83 -0.9999999999994192 1.958776339825499e-37 84 -0.9999999999995893 6.926019076028007e-38 85 -0.9999999999997097 2.4470664853523243e-38 86 -0.9999999999997947 8.650504231448578e-39 87 -0.9999999999998548 3.062343626363425e-39 88 -0.9999999999998974 1.0795595486854529e-39 89 -0.9999999999999274 3.827929532954281e-40 90 -0.9999999999999487 1.3494494358568161e-40 91 -0.9999999999999637 4.784911916192851e-41 92 -0.9999999999999744 1.6868117948210201e-41 93 -0.9999999999999819 5.926434687968075e-42 94 -0.9999999999999872 2.0812498065722716e-42 95 -0.999999999999991 7.272533761516066e-43 96 -0.9999999999999937 2.5342838525752524e-43 97 -0.9999999999999956 8.758115402030107e-44 98 -0.9999999999999969 3.0040335828963266e-44 99 -0.9999999999999978 1.0947644252537633e-44 100 -0.9999999999999984 3.755041978620408e-45 101 -0.9999999999999989 1.3684555315672042e-45 102 -0.9999999999999992 4.69380247327551e-46 103 -0.9999999999999996 8.758115402030107e-47
How about a more intractable kind of function?
for gn, fn in secmin(lambda x: math.sin(2*math.sin(x*5)), 1, 2, 3):
print(gn, fn)
1.4695450787225794 0.9841402917723512 2.2366730403515707 -0.9233648785586928 2.2368323692955316 -0.9234792748744269 2.392501302877895 -0.9068104171186121 2.3045549816727253 -0.9876003036974766 2.310199041471532 -0.9917474360135683 2.32566469457885 -0.999117573198154 2.3374102891979294 -0.9995444218579446 2.3326606991155168 -0.9999999432898783 2.332498722519485 -0.999999778440412 2.33260549707667 -0.9999999999876273 2.3326062822842495 -0.9999999999999934 2.332606300999506 -1.0 2.332606300778019 -1.0 2.3326063008887625 -1.0
math.sin(2*math.sin(2.3326063009*5))
So it found a point where the function reaches its minimum value, which is really the best we can hope for.
So this is far from a rigorous demonstration that this method will always work under certain conditions, but it does seem to have worked on the things we tried it on, so we can at least say it works some of the time. Some diagrams and maybe an animation would probably be helpful. (I did draw diagrams on paper while working this out...)