Why not use the third derivative for numerical optimization?

If Hessians are so good for optimization (see e.g. Newton's method), why stop there? Let's use the third, fourth, fifth, and sixth derivatives? Why not?

• Once you find the optimum, why look further? Indeed, what are you really trying to ask? What is your statistical question? – whuber Dec 22 '17 at 16:39
• In many cases, the limiting distribution of estimates that solve optimal estimating equations or minimize objective functions are jointly normal, so they can be characterized entirely by their first and second moments. – AdamO Dec 22 '17 at 21:41
• If you can do something, does not mean you should do it. Higher order derivatives are increasingly susceptible to noise. – Vladislavs Dovgalecs Dec 22 '17 at 22:03
• I'm voting to close this question as off-topic because it's not about statistics. It's about numerical optimization – Aksakal Dec 22 '17 at 23:40
• You haven't made a scientific breakthrough. Halley beat you by about 3 1/4 centuries. Halley, E., 1694, "A new, exact, and easy method of finding the roots of any equations generally, and that without any previous reduction" Philos. Trans. Roy. Soc. London,18, 136–145. 3rd derivative methods for optimization have existed and been studied for many years, but have not achieved great popularity. If implemented well, their greatest advantage can be an increase in robustness vs. a well-implemented Newton's method. This can be of advantage for the nastiest problems. – Mark L. Stone Dec 23 '17 at 0:04

I am interpreting the question as being "Why does Newton's method only use first and second derivatives, not third or higher derivatives?"

Actually, in many cases, going to the third derivative does help; I've done it with custom stuff before. However, in general, going to higher derivatives adds computational complexity - you have to find and calculate all those derivatives, and for multivariate problems, there are a lot more third derivatives than there are first derivatives! - that far outweighs the savings in step count you get, if any. For example, if I have a 3-dimensional problem, I have 3 first derivatives, 6 second derivatives, and 10 third derivatives, so going to a third-order version more than doubles the number of evaluations I have to do (from 9 to 19), not to mention increased complexity of calculating the step direction / size once I've done those evaluations, but will almost certainly not cut the number of steps I have to take in half.

Now, in the general case with $k$ variables, the collection of $n^{th}$ partial derivatives will number${k+n-1} \choose {k-1}$, so for a problem with five variables, the total number of third, fourth, and fifth partial derivatives will equal 231, a more than 10-fold increase over the number of first and second partial derivatives (20). You would have to have a problem that is very, very close to a fifth-order polynomial in the variables to see a large enough reduction in iteration counts to make up for that extra computational burden.

• Could you explain just how you are making use of higher derivatives? – whuber Dec 22 '17 at 16:41
• @whuber What the OP is referring to, extremely unclearly I have to admit, is Newton's method in optimization. The question really is "Why does Newton's method only use first and second derivatives, not third or higher derivatives?". It's off-topic as well as being unclear what he/she's asking, but I thought I'd just give an answer rather than vote to close for one reason or another. – jbowman Dec 22 '17 at 17:49
• +1 I think this is a good answer, but it could be improved by showing what you're up to based on the taylor expansion. – Matthew Drury Dec 22 '17 at 18:26
• As one of my professors - a very successful consultant as well - said to us once, "Whenever you think you've figured out how to build a better mousetrap, try to figure out why the 1,000 people who came up with that exact same idea before you haven't put it on the market." The whole point of using Newton is to save computation - otherwise, we'd just do exhaustive search. I assure you, adding a third derivative to a 3 dimensional problem will very, very rarely pay for the doubling of computation at each step with greatly reduced iterations unless the function is ~ a cubic. – jbowman Dec 22 '17 at 23:19
• No, it's not - it's a little deeper comment than it may first appear. The point is twofold - most ideas that appear good at first, aren't, for reasons that may not be at all obvious, and the real key to a breakthough may not be the idea itself but something that overcomes or works around the flaw in the idea. This reasoning, in effect, points that out, and tells you to look for weaknesses in the idea. It's not about giving up, it's about thinking things through, and with a critical eye at that. – jbowman Dec 23 '17 at 1:58

I don't really see what the statistical aspect of this question is, so I'll answer the optimization part.

There are 2 parts to convergence: iteration cost & iteration count

Pretty much every answer here is focusing on just the iteration cost and ignoring the iteration count. But both of them matter. An method that iterates in 1 nanosecond but that takes $10^{20}$ iterations to converge won't do you any good. And a method that blows up won't help either, no matter how cheap its iteration cost.

Let's figure out what's going on.

So: Why not use >2nd-order derivatives?

Partly because (and this is true for 2nd-order too, but more on that in a bit):

Higher-order methods generally only converge faster when near the optimum.

On the other hand, they blow up more easily when they are farther from the optimum!

(Of course, this isn't always true; e.g. a quadratic will converge in 1 step with Newton's method. But for arbitrary functions in the real world that don't have nice properties, this is generally true.)

This means that when you are farther away from the optimum, you generally want a low-order (read: first-order) method. Only when you are close do you want to increase the order of the method.

So why stop at 2nd order when you are near the root?

Because "quadratic" convergence behavior really is "good enough"!

To see why, you first have to understand what "quadratic convergence" means.

Mathematically, quadratic convergence means that, if $\epsilon_k$ is your error at iteration $k$, then the following eventually holds true for some constant $c$:

$$\lvert\epsilon_{k+1}\rvert \leq c\ \lvert\epsilon_{k}\rvert^2$$

In plain English, this means that, once you are near the optimum (important!), every extra step doubles the number of digits of accuracy.

Why? It's easy to see with an example: for $c = 1$ and $\lvert\epsilon_1\rvert = 0.1$, you have $\lvert\epsilon_2\rvert \leq 0.01$, $\lvert\epsilon_3\rvert \leq 0.0001$, etc. which is ridiculously fast. (It's super-exponential!)

Why not stop at 1st order rather than 2nd-order?

Actually, people often do this when second-order derivatives become too expensive. But linear convergence can be very slow. e.g. if you got $\epsilon_k = 0.9999999$ then you'd need maybe 10,000,000 iterations with linear convergence to get $\lvert\epsilon\rvert < 0.5$, but only 23 iterations with quadratic convergence. So you can see why there's a drastic difference between linear and quadratic convergence. This is not true for 2nd and 3rd-order convergence, for example (see next paragraph).

At this point, if you know any computer science, you understand that with 2nd-order convergence, the problem is already solved. If you don't see why, here's why: there is nothing practical to gain from tripling the number of digits every iteration instead of doubling it—what's it going to buy you? After all, in a computer, even a double-precision number has 52 bits of precision, which is around 16 decimal digits. Maybe it will decrease the number of steps you require from 16 to 3... which sounds great, until you realize it comes at the price of having to compute third derivatives at each iteration, which is where the curse of dimensionality hits you hard. For a $6$-dimensional problem, you just paid a factor of $6$ to gain a factor of $\approx 5$, which is dumb. And in the real world problems have at least hundreds of dimensions (or even thousands or even millions), not merely $6$! So you gain a factor of maybe 20 by paying a factor of, say, 20,000... hardly a wise trade-off.

But again: remember the curse of dimensionality is half the story.

The other half is that you generally get worse behavior when you're far from the optimum, which generally adversely affects the number of iterations you have to do.

Conclusion

In a general setting, higher-order methods than 2 are a bad idea. Of course, if you can bring additional helpful assumptions to the table (e.g. perhaps your data does resemble a high-degree polynomial, or you have ways of bounding the location of the optimum, etc.), then maybe you can find that they are a good idea—but that will be a problem-specific decision, and not a general rule of thumb to live by.

• Great answer, but I think the Abel-Ruffini theorem is a red herring. First of all, we are speaking about multivariate problems, so computing zeros of univariate polynomials is at most an easy subproblem of limited interest. And, more importantly, it doesn't matter if there is a closed formula for the solution or not: in practice, as far as I know, people don't use closed formulas even for degree-4 polynomials. They're just too long and complicated and unstable. Zeros of polynomials are computed numerically, in practice (using QR on the companion matrix). – Federico Poloni Dec 23 '17 at 11:10
• @FedericoPoloni: Yeah, the same thoughts came to my mind when I was deciding to put it in. I didn't have it in originally... I thought maybe I should put it in as just another example of why higher degrees can have unexpected problems. But I guess I'll take it out again if it's unhelpful, thanks for the comment. – Mehrdad Dec 23 '17 at 11:15
• @FedericoPoloni: P.S. while we're on the topic of numerical computation, you might find Sturm functions interesting (if you haven't heard of them already). – Mehrdad Dec 23 '17 at 11:22

Even calculating Hessians is quite a bit of work: $$H = \begin{bmatrix} \dfrac{\partial^2 f}{\partial x_1^2} & \dfrac{\partial^2 f}{\partial x_1\,\partial x_2} & \cdots & \dfrac{\partial^2 f}{\partial x_1\,\partial x_n} \\[2.2ex] \dfrac{\partial^2 f}{\partial x_2\,\partial x_1} & \dfrac{\partial^2 f}{\partial x_2^2} & \cdots & \dfrac{\partial^2 f}{\partial x_2\,\partial x_n} \\[2.2ex] \vdots & \vdots & \ddots & \vdots \\[2.2ex] \dfrac{\partial^2 f}{\partial x_n\,\partial x_1} & \dfrac{\partial^2 f}{\partial x_n\,\partial x_2} & \cdots & \dfrac{\partial^2 f}{\partial x_n^2} \end{bmatrix}.$$

Now see how the third derivative looks like: $$\partial H/\partial x=\begin{bmatrix} \dfrac{\partial H}{\partial x_1}\\ \dfrac{\partial H}{\partial x_2}\\ \vdots\\ \dfrac{\partial H}{\partial x_n} \end{bmatrix}$$ This is a three dimensional matrix. Here's how its elements look like: $$(\partial H/\partial x)_{ijk}=\dfrac{\partial^3 f}{\partial x_i\partial x_j\partial x_k}$$

Sixth's derivative will be six dimensional matrix: $$\dfrac{\partial^6 f}{\partial x_i\partial x_j\partial x_k\partial x_l\partial x_m\partial x_n}$$

Usually, the trade-off is not favorable for going after higher than Hessian. I mean the trade-off between potential gain in speed through using higher order approximations vs. the noise amplification. You always have noise in inputs because we're talking about statistical applications. This noise will be amplified by the derivatives.

If you play golf then the analogy in optimization is to first swing trying to get to the green, not worry to much about a hole. Once, on the green, we'll putt aiming a hole.

Typically, when you analyze the effectiveness of such algorithms, you'll find results such as one step of a fourth order algorithm having roughly the same effectiveness as two steps of a second order algorithm.

So the choice of which algorithm to use is relatively simple: if one step of the fourth order algorithm takes twice as much work or more than one step of the second order algorithm, you should use the latter instead.

That is the typical situation for these sorts of methods: the classical algorithm has the optimal work-to-effectiveness ratio for general problems. While there are occasional problems where a higher order approach is unusually easy to compute and can outperform the classical variant, they are relatively uncommon.

You can think of the order of derivatives as the order of a polynomial approximation to the function. Most optimization routines rely on convexity. A quadratic polynomial will be convex/concave everywhere whereas a 3rd order or higher polynomial will not be convex everywhere. Most optimization routines rely on successive approximations of convex functions with quadratics for this reason. A quadratic approximation that is convex requires a positive definiteness condition to be imposed in order for the quadratic to be convex.

• No, quadratics are not necessarily convex or concave (think of $x^2-y^2$). – Dirk Dec 22 '17 at 21:32
• @Dirk $x^2-y^2$ equal to what? – Ovi Dec 23 '17 at 4:02
• It's a quadratic function but neither convex nor concave. – Dirk Dec 23 '17 at 9:17
• @Dirk yes you are right, I should've added a positive semi-definite caveat. I will add that to my answer. – Lucas Roberts Dec 24 '17 at 19:04

Let me be the only one here defending 3rd order methods for SGD convergence, but definitely not in the entire space what would need $$\approx dim^3/6$$ coefficients, but e.g. in just a single direction, which needs only a single additional coefficient if already having 2nd order model in this direction.

Why single direction 3rd order model can be beneficial? For example because close to zero second derivative in this direction basically means two alternative scenarios: plateau or inflection point - only the former requires larger step size, and 3rd derivative allows to distinguish them.

I believe we will go toward hybrid multi-order methods: 2nd order method in a low dimensional subspace e.g. from PCA of recent gradients, what still allows for free 1st order simultaneous gradient descent toward part of gradient orthogonal to this subspace ... and additionally I would add e.g. 3rd order model for a single most relevant direction.