Quantitative Methods · Methods
24Numerical Optimisation in Practice
A desk calibrates a two-exponential kernel to a measured autocorrelation every morning: a fast component with a time scale of about three days and a slow one of about two months. For most of a year the slow time scale reads between 53 and 70 days. One Monday it reads 88, where Friday’s read 57, and the fit to the data is, if anything, slightly better. Nothing is wrong with the optimiser and nothing happened in the data: the objective has a long, nearly flat valley along which the slow time scale and the weight of the fast component trade off, and a little noise moved the bottom of the valley. The same objective has a mirror-image valley, reached by swapping the labels of the two components, and dead ends at the boundary. This chapter is about the optimisers that turn such objectives into numbers (gradient methods and their acceleration, Newton and quasi-Newton methods, nonlinear least squares, splitting and proximal methods, stochastic gradients) and about the practice that makes a calibration trustworthy: transforms, multistart, diagnostics of identifiability, and a penalty toward yesterday’s answer.
24.1 First-order methods
Definition 24.1 (Gradient descent, line search)
Gradient descent iterates . A line search chooses the step along the descent direction, for instance to satisfy the Wolfe conditions: sufficient decrease, , and a curvature condition, , with .
Proposition 24.2 (Rates of gradient descent)
If is -strongly convex with -Lipschitz gradient, gradient descent with step satisfies : linear convergence at a rate set by the condition number . For convex without strong convexity the gap falls like .
Proof. With step the descent lemma gives , and strong convexity gives ; combining, . ∎
Definition 24.3 (Momentum method, Nesterov acceleration)
A momentum method adds a fraction of the previous step to the gradient step. Nesterov acceleration evaluates the gradient at an extrapolated point, , , which improves the convex rate to and the strongly convex rate to about per iteration.
On the quadratic , plain gradient descent from with step removes the stiff coordinate in one step and then shrinks the other by per step, so the objective falls from 50.5 to : 0.067 after 100 iterations and 0.0090 after 200 (the proposition’s bound, , is conservative by the square); Nesterov’s method reaches in 100 (Figure 24.1). Condition number is the enemy of first-order methods, and ill-conditioning is the normal state of calibrations whose parameters trade off.
24.2 Newton, quasi-Newton and nonlinear least squares
Definition 24.4 (Newton’s method, quasi-Newton method, BFGS)
Newton’s method steps to the minimum of the local quadratic model, . A quasi-Newton method replaces the Hessian by an approximation updated from gradient differences; the BFGS method (Broyden, Fletcher, Goldfarb and Shanno, 1970) updates the inverse approximation by , with the step, the gradient change and ; L-BFGS (Liu and Nocedal, 1989) keeps only the last few pairs.
Near a minimum with a positive-definite Hessian Newton’s method converges quadratically: the number of correct digits doubles at each step. The BFGS update keeps positive definite as long as , which a line search satisfying the Wolfe curvature condition guarantees, since then . On the quadratic above, Newton’s method takes one step and L-BFGS seven.
Definition 24.5 (Nonlinear least squares, Gauss–Newton, Levenberg–Marquardt)
Nonlinear least squares minimises for a residual vector . The Gauss–Newton method approximates the Hessian by , the Jacobian of , and solves . The Levenberg–Marquardt algorithm (Levenberg, 1944; Marquardt, 1963) damps it, , raising after a failed step (toward scaled gradient descent) and lowering it after a success (toward Gauss–Newton).
Least-squares calibration is the chapter’s case: the residuals are the model’s autocorrelation minus the measured one at lags 1 to 30. The weight of the fast component lives in and the time scales are positive, so the optimiser works on unconstrained transforms (logit and log) and never proposes an illegal parameter. at the solution also measures identifiability: on the first day, its condition number is 10.5 and the implied correlation between the fast weight and the slow time scale is 0.96, the numerical signature of the flat valley.
Definition 24.6 (Calibration, identifiability, multistart)
Calibration is the choice of a model’s parameters to reproduce observed quantities (prices, quotes, moments), usually by nonlinear least squares. A parameter is locally unidentifiable, and identifiability fails, when the objective is flat in some direction at the optimum, so that data cannot pin it down. Multistart runs a local optimiser from many starting points to explore distinct minima.
From 100 random starts on the first day’s curve, 28 runs end at the best fit (root-mean-squared error 0.0096), 35 end at its mirror image with the labels of the two components swapped (the same curve, the same error, different parameters), 26 run off to the boundary where the slow time scale reaches the 1 000-day cap and the fit error is at least twice as large, and 11 stop elsewhere (Figure 24.2). A calibration run from one fixed start never learns any of this.
The valley itself is visible in a profile (Figure 24.3): fixing the slow time scale and re-optimising the other two parameters, the fit error is at 62 days, at 42 and at 89, against a noise level of . Anything between about 45 and 80 days fits within the noise.
Refit every day for 250 days from the same start, the slow time scale has a day-to-day change with median 4.5 days and maximum 30.7 (day 189: from 56.8 to 87.6 days, fit error from 0.00926 to 0.00921). With a Tikhonov penalty toward yesterday’s parameters (weight 0.01 on the squared change of each transformed parameter, the fit starting from yesterday’s answer), the median change is 1.9 days and the maximum 9.5, the standard deviation of the slow time scale over the year halves from 5.5 to 2.7 days, and the average fit error rises from 0.00945 to 0.00952, by 0.7% (Figure 24.4).
24.3 Operator splitting
Definition 24.7 (Proximal operator, proximal gradient method)
The proximal operator of a convex function is ; for it is soft thresholding. The proximal gradient method minimises , smooth, by ; FISTA (Beck and Teboulle, 2009) adds Nesterov’s extrapolation.
The lasso of chapter 16 is the textbook case. On a regression of 200 observations on 50 predictors, two of them nearly collinear, the plain proximal gradient method (ISTA) is still from the optimal objective after 50 iterations and after 100; FISTA is about away after 50 and below after 65 (Figure 24.5).
Definition 24.8 (Alternating direction method of multipliers)
The alternating direction method of multipliers (ADMM) minimises subject to by alternating , , (Boyd, Parikh, Chu, Peleato and Eckstein, 2011).
ADMM converges for closed convex and whenever a solution exists, at a modest rate, and each step is only as hard as the two proximal operators. For the nearest correlation matrix of chapter 23, is the distance to the given matrix plus the semidefinite constraint (an eigenvalue clip) and the unit-diagonal constraint. On chapter 23’s broken matrix ADMM needs 91 iterations to the alternating projections’ 74, and the two answers agree to .
24.4 Stochastic gradients
Definition 24.9 (Stochastic gradient descent)
Stochastic gradient descent minimises by steps along the gradient of one randomly chosen term (or a small batch), .
Each step is an unbiased but noisy estimate of a gradient step. With a constant step the iterates hover in a noise ball around the minimum; Robbins and Monro (1951) showed that steps with and , such as , converge. On a least-squares problem with 1 000 observations, 20 000 steps land 0.009 from the least-squares solution with Robbins–Monro steps and 0.113 with a constant step of 0.05. Adaptive methods such as Adam (Kingma and Ba, 2014) rescale each coordinate’s step by running moments of its gradients; they are the default for the neural networks of One Quant Book 9, and a poor default for a three-parameter calibration, where Levenberg–Marquardt converges in a few dozen steps.
24.5 Tutorial: the Monday the fit moved
Goal. Calibrate a two-exponential kernel every day, explore its minima by multistart, diagnose the flat valley, and stabilise the daily fit with a penalty toward yesterday. End state: Figures 24.2, 24.3 and 24.4; the median jump of 4.5 and 1.9 days and the 0.7% fit cost.
Levenberg–Marquardt with its damping schedule.
def levenberg_marquardt(resid, x0, jac=None, tol: float = 1e-10, max_iter: int = 500, lam: float = 1e-3) -> dict: """Minimise |r(x)|^2 / 2: solve (J'J + lam diag(J'J)) dx = -J'r, accept if the cost falls (then lam / 10), otherwise lam * 10.""" x = np.asarray(x0, dtype=float).copy() r = resid(x) cost = 0.5 * float(r @ r) J = jac(x) if jac else numerical_jacobian(resid, x) for it in range(1, max_iter + 1): A = J.T @ J g = J.T @ r dx = np.linalg.solve(A + lam * np.diag(np.maximum(np.diag(A), 1e-12)), -g) x_new = x + dx r_new = resid(x_new) c_new = 0.5 * float(r_new @ r_new) if np.isfinite(c_new) and c_new < cost: x, r, lam = x_new, r_new, lam / 10 done = cost - c_new < tol * max(cost, 1e-300) cost = c_new J = jac(x) if jac else numerical_jacobian(resid, x) if done or np.linalg.norm(dx) < tol * (1 + np.linalg.norm(x)): return {"x": x, "cost": cost, "iters": it, "jac": J, "status": "optimal"} else: lam *= 10 if lam > 1e12: break return {"x": x, "cost": cost, "iters": max_iter, "jac": J, "status": "stalled"}Listing 24.1. Levenberg–Marquardt for nonlinear least squares. code/firm/optim/firm_optim.py The penalty toward yesterday and the identifiability diagnostic.
def identifiability(J) -> dict: """Condition number and smallest singular value of the Jacobian, and the parameter correlation implied by (J'J)^-1: near-collinear columns mean a flat valley in the objective.""" J = np.asarray(J, dtype=float) sv = np.linalg.svd(J, compute_uv=False) cov = np.linalg.pinv(J.T @ J) d = np.sqrt(np.maximum(np.diag(cov), 1e-300)) return {"cond": float(sv[0] / sv[-1]), "smallest_sv": float(sv[-1]), "corr": cov / np.outer(d, d)} def tikhonov(resid, prev, weight): """Append sqrt-weighted deviations from the previous parameters to the residuals: the fit then minimises |r(x)|^2 + sum_j weight_j (x_j - prev_j)^2.""" prev, w = np.asarray(prev, dtype=float), np.sqrt(np.asarray(weight, dtype=float)) def r(x): return np.concatenate([resid(x), w * (np.asarray(x) - prev)]) return rListing 24.2. Identifiability diagnostics and the Tikhonov penalty. code/firm/optim/firm_optim.py - Run
daily(0.0),daily(0.01),multistart(),valley()anddiagnostics()inqm_optim.py, thenfig_optim.py.
What to change next. Order the time scales by construction () to remove the mirror valley; choose the penalty weight by cross-validating tomorrow’s fit error.
24.6 Build: the calibration harness
Purpose. Every model the miniature firm calibrates daily (kernels, volatility surfaces, term structures) runs through one harness that transforms, multistarts, penalises and logs.
Interface. gradient_descent, newton, lbfgs, levenberg_marquardt, prox_l1, proximal_gradient, admm_nearest_correlation, sgd, Bounded, multistart, identifiability, tikhonov.
Rules. Parameters are optimised in unconstrained coordinates; every production calibration logs its start, iterations, fit error, the condition number of and the change from yesterday; a new minimum far from yesterday’s is flagged, not published silently; multistart runs weekly on each model.
Acceptance tests. code/firm/optim/tests/: the linear rate of gradient descent, acceleration, Newton on Rosenbrock; L-BFGS on Rosenbrock; Levenberg–Marquardt recovers an exponential fit and a heavy penalty pins it; soft thresholding and the lasso by FISTA; ADMM’s nearest correlation matrix; SGD with Robbins–Monro steps; transforms round-trip; multistart ordering.
Stretch. Box constraints by projection in L-BFGS (L-BFGS-B); automatic differentiation of the residuals (chapter 28); a trust-region Levenberg–Marquardt.
Sources and further reading
- K. Levenberg, “A method for the solution of certain non-linear problems in least squares”, Quarterly of Applied Mathematics 2, 1944; D. W. Marquardt, “An algorithm for least-squares estimation of nonlinear parameters”, Journal of the Society for Industrial and Applied Mathematics 11, 1963.
- D. C. Liu and J. Nocedal, “On the limited memory BFGS method for large scale optimization”, Mathematical Programming 45, 1989.
- H. Robbins and S. Monro, “A stochastic approximation method”, Annals of Mathematical Statistics 22, 1951.
- A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding algorithm for linear inverse problems”, SIAM Journal on Imaging Sciences 2, 2009.
- S. Boyd, N. Parikh, E. Chu, B. Peleato and J. Eckstein, “Distributed optimization and statistical learning via the alternating direction method of multipliers”, Foundations and Trends in Machine Learning 3, 2011.
- D. P. Kingma and J. Ba, “Adam: a method for stochastic optimization”, arXiv:1412.6980, 2014.
- J. Nocedal and S. J. Wright, Numerical Optimization, 2nd ed., Springer, 2006.
24.7 Exercises
Exercise 24.1 ★
How many gradient-descent iterations reduce the error of a quadratic with condition number 1 000 by a factor of ? And with Nesterov acceleration, roughly?
Solution
Solution of Exercise 24.1.
With the bound , iterations; with acceleration, about .
Exercise 24.2 ★
Compute for , .
Solution
Solution of Exercise 24.2.
Soft thresholding by 1: .
Exercise 24.3 ★
Why does the two-exponential kernel have two minima with exactly the same fit?
Solution
Solution of Exercise 24.3.
Swapping for gives the same function of the lag: the model is invariant under relabelling its components, so every minimum has a mirror image. Ordering the time scales in the parameterisation removes it.
Exercise 24.4 ★★
Show that one Newton step from any point minimises a strictly convex quadratic.
Solution
Solution of Exercise 24.4.
For with positive definite, , the unique minimiser, whatever .
Exercise 24.5 ★★
Show that the Wolfe curvature condition implies for a descent direction.
Solution
Solution of Exercise 24.5.
With , : , since and .
Exercise 24.6 ★★
Do the steps satisfy the Robbins–Monro conditions? And ?
Solution
Solution of Exercise 24.6.
: the sum diverges but so does the sum of squares (), so not in the classical form (though averaging the iterates rescues it). : the sum diverges and the sum of squares converges, so yes.
Exercise 24.7 ★★★
Coding. Run the daily calibration with penalty weights 0, 0.001, 0.01 and 0.1 and tabulate the median day-to-day change of the slow time scale against the average fit error.
Solution
Solution of Exercise 24.7.
daily(w) for , 0.001, 0.01, 0.1: median change 4.48, 3.93, 1.87 and 0.42 days; average RMSE () 9.450, 9.454, 9.520 and 9.643. Stability is cheap up to 0.01 and then starts to cost fit.
Exercise 24.8 ★★★
Find the flaw. “Our calibration always converges, with a gradient norm below , so its parameters are reliable.”
Solution
Solution of Exercise 24.8.
Convergence to a stationary point says nothing about which minimum was found (mirror images, boundary solutions), nor about how flat the objective is there: parameters in a flat valley move with noise while the gradient is zero at each day’s point. Check multistart, the condition number of and the day-to-day stability.
24.8 Problem: The Monday the Fit Moved
Problem 24.1
Weekend problem — a calibration’s flat valley
Each day a desk fits , , to a measured autocorrelation whose true parameters are , and days, measured with noise of standard deviation 0.01. The fit is by Levenberg–Marquardt in logit and log coordinates. Two hundred and fifty days are simulated.
Part I — One day.
- From the fixed start, what does the first day’s fit give, and with what error?
- What do 100 random starts find?
- What is the mirror-image solution, and why does it exist?
- What are the condition number of and the correlation of the parameters?
- Over what range of the slow time scale does the fit stay within the noise?
Part II — Every day.
- What are the median and largest day-to-day changes of the slow time scale?
- What happened on the day of the largest change?
- With a penalty toward yesterday of weight 0.01, what are they?
- What does the penalty cost in fit error?
- How much does the year’s spread of the slow time scale fall?
Part III — Methods.
- How do gradient descent and Nesterov’s method compare on a quadratic with condition number 100?
- How many steps do Newton and L-BFGS need on it?
- How do ISTA and FISTA compare on the lasso?
- How do ADMM and alternating projections compare on the nearest correlation matrix?
- What do Robbins–Monro steps buy over a constant step?
Part IV — Judgement.
- Which parameter would you publish with a warning, and what warning?
- How would you choose the penalty weight?
- How would you remove the mirror valley?
- State the named result: the day-to-day parameter jump with and without the penalty, and the fit error the penalty costs.
- In one sentence: what does a converged calibration not tell you?
Solution
Solution of Problem 24.1.
1. , and days, RMSE 0.0096. 2. 28 at the best fit, 35 at its mirror image, 26 at the 1 000-day boundary with at least twice the error, 11 elsewhere. 3. : the same curve with the components’ labels swapped. 4. Condition number 10.5; correlation 0.96 between the fast weight and the slow time scale. 5. From about 45 to 80 days (profile error at 62 days, 12.3 at 42, 12.9 at 89, against noise of 10). 6. Median 4.5 days, maximum 30.7. 7. Day 189: the slow time scale went from 56.8 to 87.6 days and the fit error from 0.00926 to 0.00921. 8. Median 1.9 days, maximum 9.5. 9. The average RMSE rises from 0.00945 to 0.00952, 0.7%. 10. Its standard deviation falls from 5.5 to 2.7 days. 11. Gradient descent reaches 0.067 after 100 iterations and 0.0090 after 200; Nesterov after 100. 12. One Newton step (two iterations including the check); L-BFGS seven. 13. ISTA is from the optimum after 50 iterations; FISTA about after 50 and below after 65. 14. 91 ADMM iterations against 74 projections, the same matrix to . 15. Convergence: 0.009 from the least-squares solution after 20 000 steps, against 0.113 with a constant step. 16. The slow time scale, with its interval along the valley (about 45 to 80 days) and its day-to-day change. 17. By cross-validation: the weight that minimises tomorrow’s fit error of today’s parameters, or a weight that caps the day-to-day change at a tolerance the users accept. 18. Parameterise (and ), so only one labelling is representable. 19. Named result: the Monday the fit moved: refit daily from a fixed start, the slow time scale jumps by a median of 4.5 days (maximum 30.7); with a penalty toward yesterday (weight 0.01) by 1.9 (maximum 9.5), at a cost of 0.7% in fit error. 20. Which of several equally good parameter sets it found, and how precisely the data determine them.
24.9 Interview questions
Interview question 24.1 ★ researcher, mle
Why does gradient descent slow down on ill-conditioned problems, and what do you do about it?
Solution
Solution of Interview question 24.1.
The step must be small enough for the stiffest direction (), so progress along the flattest direction is only per step: the rate is set by the condition number. Remedies: rescale the variables, precondition, use curvature (Newton, quasi-Newton, Gauss–Newton), or accelerate (Nesterov, momentum).
What the interviewer is looking for: Condition number as the rate, and at least two remedies.
Interview question 24.2 ★★ researcher, developer
Your model calibration’s parameters jump from day to day while the fit error stays the same. What is happening, and what do you do?
Solution
Solution of Interview question 24.2.
The objective is flat in some direction (an identifiability problem), or has several equivalent minima, and noise moves the answer. Diagnose with the condition number of , a profile of the objective and multistart; fix with a penalty toward yesterday, reparameterisation, fixing a poorly identified parameter, or more informative data.
What the interviewer is looking for: Identifiability diagnosis, and a stabilisation that costs little fit.
Interview question 24.3 ★★ researcher
Explain Levenberg–Marquardt. Why not plain Gauss–Newton?
Solution
Solution of Interview question 24.3.
It solves : small gives Gauss–Newton (fast near the solution), large a short scaled gradient step (safe far from it); adapts to success. Plain Gauss–Newton can take huge steps when is nearly singular, exactly in flat valleys.
What the interviewer is looking for: Interpolation between Gauss–Newton and gradient descent, and the singular .
Interview question 24.4 ★★ mle
Why does SGD with a constant learning rate not converge, and what are the Robbins–Monro conditions?
Solution
Solution of Interview question 24.4.
Each stochastic gradient has variance that does not vanish at the optimum, so a constant step keeps the iterate bouncing in a noise ball of radius proportional to the step. With (enough total movement) and (vanishing noise), the iterates converge.
What the interviewer is looking for: Noise at the optimum and the two conditions.
Interview question 24.5 ★★ researcher, mle
What is a proximal operator, and why is it useful for the lasso?
Solution
Solution of Interview question 24.5.
minimises ; for the norm it is soft thresholding. The lasso’s non-smooth penalty is handled exactly by the prox, while the smooth least-squares part uses its gradient: the proximal gradient method (ISTA, FISTA) and coordinate descent both exploit it.
What the interviewer is looking for: Soft thresholding and the smooth-plus-simple split.
Interview question 24.6 ★★★ researcher
Why does BFGS need a line search satisfying the Wolfe conditions?
Solution
Solution of Interview question 24.6.
The BFGS update keeps the inverse-Hessian approximation positive definite only if ; the Wolfe curvature condition guarantees it, while a mere decrease condition does not. Without it the next direction may not descend.
What the interviewer is looking for: Curvature condition implies positive definiteness.
Terms defined in this chapter
- Alternating direction method of multipliers
- Calibration, identifiability, multistart
- Gradient descent, line search
- Momentum method, Nesterov acceleration
- Newton’s method, quasi-Newton method, BFGS
- Nonlinear least squares, Gauss–Newton, Levenberg–Marquardt
- Proximal operator, proximal gradient method
- Stochastic gradient descent