Geometry of nonlinear least squares with applications to sloppy models and optimization

Parameter estimation by nonlinear least-squares minimization is a common problem that has an elegant geometric interpretation: the possible parameter values of a model induce a manifold within the space of data predictions. The minimization problem is then to find the point on the manifold closest to the experimental data. We show that the model manifolds of a large class of models, known as sloppy models, have many universal features; they are characterized by a geometric series of widths, extrinsic curvatures, and parameter-effect curvatures, which we describe as a hyper-ribbon. A number of common difficulties in optimizing least-squares problems are due to this common geometric structure. First, algorithms tend to run into the boundaries of the model manifold, causing parameters to diverge or become unphysical before they have been optimized. We introduce the model graph as an extension of the model manifold to remedy this problem. We argue that appropriate priors can remove the boundaries and further improve the convergence rates. We show that typical fits will have many evaporated parameters unless the data are very accurately known. Second, “bare” model parameters are usually ill-suited to describing model behavior; cost contours in parameter space tend to form hierarchies of plateaus and long narrow canyons. Geometrically, we understand this inconvenient parametrization as an extremely skewed coordinate basis and show that it induces a large parameter-effect curvature on the manifold. By constructing alternative coordinates based on geodesic motion, we show that these long narrow canyons are transformed in many cases into a single quadratic, isotropic basin. We interpret the modified Gauss-Newton and Levenberg-Marquardt fitting algorithms as an Euler approximation to geodesic motion in these natural coordinates on the model manifold and the model graph, respectively. By adding a geodesic acceleration adjustment to these algorithms, we alleviate the difficulties from parameter-effect curvature, improving both efficiency and success rates at finding good fits.

Figure 1

(a) Fitting a nonlinear function to data, in this case the sum of two exponentials to three data points. Fit A has rate constants that decay too quickly, resulting in a poor fit; B is an improvement over fit A, although the rates are too slow; the best fit minimizes the cost (the sum of the squares of the residuals, which are deviations of model from data points). (b) Contours of constant cost in parameter space. Note the “plateau” in the region of large rates where the model is essentially independent of parameter changes. Note also the long, narrow canyon at lower rates, characteristic of a sloppy model. The sloppy direction is parallel to the canyon and the stiff direction is against the canyon wall. (c) Model predictions in data space. The experimental data are represented by a single point. The set of all possible fitting parameters induces a manifold of predictions in data space. The best fit is the point on the manifold nearest to the data. The plateau in (b) here is the small region around the short cusp near the corner. To help visualize the three-dimensional structure, an animation of this manifold rotating in three dimensions in available in the online supplemental material [23].

Figure 2

Hessian eigenvalues for three sloppy models. Note the extraordinarily large range of eigenvalues ($15$–$17$ orders of magnitude, corresponding to valley aspect ratios of ${10}^7$–${10}^9$) in Fig. 1b. Notice also the roughly equal fractional spacing between eigenvalues—there is no clean separation between important (stiff) and irrelevant (sloppy) direction in parameter space. (a) The model formed by summing six exponential terms with rates and amplitudes. We use this model to investigate curvature in Sec. 7 and as a test problem to compare algorithms in Sec. 8e. (b) The linear problem of fitting polynomials is sloppy with the Hessian given by the Hilbert matrix. (c) A more practical model from systems biology of signaling the epidermal growth factor in rat pheochromocytoma (PC12) cells [2], which also has a sloppy eigenvalue spectrum. Many more examples can be found in [6, 8].

Figure 3

Skewed coordinates. A sloppy model is characterized by a skewed coordinate mesh on the manifold. The volume of the parallel-piped is given by the determinant of the metric, which is equal to the product of the eigenvalues. Because sloppy models have many tiny eigenvalues, these volumes can be very small with extremely skewed coordinates. Our toy model has extremely skewed coordinates where the parameters are nearly equal (near the fold line). Most of the manifold is covered by regions where the coordinates are less skewed, which corresponds to a very small region in parameter space.

Figure 4

Fixing a few data points greatly restricts the possible range of the model behavior between those data points (lower). This is a consequence of interpolation of analytic functions. In this case, $f(t)$ is a sum of three exponentials with six parameters (amplitudes and rates). Shown above is a three-dimensional slice of possible models plotted in data space, with the value of $f(0)$ fixed to 1 and the value of $f(1)$ fixed to $1/e$. With these constraints we are left with a four-dimensional surface, meaning that the manifold of possible data shown here is indeed a volume. However, from a carefully chosen perspective (upper right), this volume can be seen to be extremely thin—in fact, most of its apparent width is curvature of the nearly two-dimensional sheet, evidenced by being able to see both the top (green) and bottom (black) simultaneously. (An animation of points in this volume rotating in three-dimensional space is available in the online supplemental material [23].) Generic aspects of this picture illustrate the difficulty of fitting nonlinear problems. Geodesics in this volume are just straight lines in three dimensions. Although the manifold seems to be only slightly curved, its extreme thinness means that geodesics travel very short distances before running into model boundaries, necessitating the diagonal cutoff in Levenberg-Marquardt algorithms as well as the priors discussed in Sec. 5.

Figure 5

The possible values of a model at intermediate time points are restricted by interpolating theorems. Taking cross sections of the model manifold corresponds to fixing the model values at a few time points, restricting the possible values at the remaining times. Therefore, the model manifold will have a hierarchy of progressively thinner widths, much like a hyper-ribbon.

Figure 6

(a) Geodesic cross-sectional widths of an eight-dimensional model manifold along the eigendirections of the metric from some central point, together with the square root of the eigenvalues (singular values of the Jacobian) [22]. Notice the hierarchy of these data-space distances—the widths and singular values each spanning around four orders of magnitude. To a good approximation, the cross-sectional widths are given by singular values. In the limit of infinitely many exponential terms, this model becomes linear. (b) Geodesic cross-sectional widths of a feed-forward artificial neural network. Once again, the widths nicely track the singular values.

Figure 7

(a) Falling off the edge of the model manifold. The manifold in data space defines a “natural” direction, known as the Gauss-Newton direction, in which an algorithm will try to follow to the best fit. Often this direction will push parameters toward the edge of the manifold. (b) Gradient and Gauss-Newton directions in parameter space. The manifold edge corresponds to infinite values of the parameters. Following the Gauss-Newton direction to the edge of the manifold will cause parameters to evaporate while on the plateau. While in a canyon, however, the Gauss-Newton direction gives the most efficient direction to the best fit.

Figure 8

The effect of the damping parameter is to produce a new metric for the surface induced by the graph of the model vs the input parameters. (a) Model graph, $\lambda =0$. If the parameter is zero, then the resulting graph is simply the original model manifold, with no extent in the parameter directions. Here we see a flat two-dimensional cross section; the $z$ axis is a parameter value multiplied by $\sqrt{\lambda }=0$. (b) Model graph $\lambda \ne 0$. If the parameter is increased, the surface is “stretched” into a higher-dimensional embedding space. This is an effective technique for removing the boundaries, as no such boundary exists in the model graph. However, this comes at a cost of removing the geometric connection between the cost function and the structure of the surface. For very large damping parameters, the model graph metric becomes a multiple of the parameter space metric, which rotates the Gauss-Newton direction into the gradient direction. The damping term, therefore, interpolates between the parameter space metric and the data space metric. A three-dimensional animation of this figure is available in the online supplemental material [23].

Figure 9

(a) Gauss-Newton directions. The Gauss-Newton direction is prone to pointing parameters toward infinity, especially in regions where the metric has very small eigenvalues. (b) Rotated Gauss-Newton directions. By adding a small damping parameter to the metric, the Gauss-Newton direction is rotated into the gradient direction. The amount of rotation is determined by the eigenvalues of the metric at any given point. Here, only a few points are rotated significantly. (c) Gradient directions. For large values of the damping parameter, the natural direction is rotated everywhere into the gradient direction.

Figure 10

The final cost is plotted against the number of Jacobian evaluations for five strengths of priors. After minimizing with priors, the priors are removed and a maximum of 20 further Jacobian evaluations are performed. The prior strength is measured by $p$, with $p=0$ meaning no prior. The success rate is $R$. The strongest priors converge the fastest, with medium strength priors showing the highest success rate.

Figure 11

(a) Extended geodesic coordinates. The parameters of a model are not usually well suited to describing the behavior of a model. By considering the manifold induced in data space, one can construct more natural coordinates based on geodesic motion that are more well-suited to describing the behavior of a model (black grid). These coordinates remove all parameter-effect curvature and are known as extended geodesic coordinates. Note that we have moved the data point so that the best fit is not so near a boundary in this picture. (b) Cost contours in extended geodesic coordinates. Although the summing exponential model is nonlinear, that nonlinearity does not translate into large extrinsic curvature. This type of nonlinearity is known as parameter-effect curvature, which the geodesic coordinates remove. This is most dramatically illustrated by considering the contours of constant cost in geodesic coordinates. The contours are nearly circular all the way out to the fold line and the boundary, where the rates are infinite.

Figure 12

By changing the value of the Levenberg-Marquardt parameter, the course of the geodesics on the corresponding model graph are deformed, in turn distorting the shape of the cost contours in the geodesic coordinates. (a) $\lambda =0$ is equivalent to the model manifold. The cost contours for a relatively flat manifold, such as that produced by the sum of two exponentials, are nearly perfect, concentric circles. The geodesics can be evaluated up to the boundary of the manifold, at which point the coordinates are no longer defined. Here we can clearly see the stiff, long manifold direction (vertical) and the sloppy, thin manifold direction (horizontal). (b) Small $\lambda$ ($\lambda$ much smaller than any of the eigenvalues of the metric) will produce cost contours that are still circular, but the manifold boundaries have been removed. In this case, the fold line has disappeared, and cost contours that ended where parameters evaporated now stretch to infinity. (c) Moderate $\lambda$ creates cost contours that begin to stretch in regions where the damping parameter significantly affects the eigenvalue structure of the metric. The deformed cost contours begin to take the plateau and canyon structures of the contours in parameter space. (d) Large $\lambda$ effectively washes out the information from the model manifold metric, leaving just a multiple of the parameter space metric. In this case, the contours are those of parameter space—a long narrow curved canyon around the best fit. This figure analogous to Fig. 1b, although the model here is a more sloppy (and more realistic) example. An animation of the transition from small to large damping parameter is available in the online supplemental material [23].

Figure 13

Intrinsic and extrinsic curvature. Intrinsic curvature is inherent to the manifold and cannot be removed by an alternative embedding. A model that is the sum of two exponential terms has all types of curvature. This is the same model manifold as in Fig. 1c, viewed from an alternative angle to highlight the curvature. From this viewing angle, the extrinsic curvature becomes apparent. This is also an example of intrinsic curvature. An animation of this surface rotating in three dimensions is available in the online supplemental material [23].

Figure 14

A ruled surface has no intrinsic curvature; however, it may have extrinsic curvature. The model manifold formed from a single exponential rate and amplitude is an example of a ruled surface. This model could be isometrically embedded in another space to remove the curvature. An animation of this surface rotating in three-dimensional space is available in the online supplemental material [23].

Figure 15

Geodesic curvature. A direction on a curved surface defines a geodesic. The deviation of the geodesic from a straight line in the embedding space is measured by the geodesic curvature. It is the inverse radius of the circle fit to the geodesic path at the point. A three-dimensional animation of this surface is available in the online supplemental material [23].

Figure 16

Shape operator. Specifying a direction normal to a curved surface, $\mathrm{n̂}$, defines a shape operator. The eigenvalues of the shape operator are the principal curvatures and the corresponding eigenvectors are the directions of principal curvature. A three-dimensional animation of this surface is available in the online supplemental material [23].

Figure 17

(a) Linear grid. A sloppy linear model may have a skewed coordinate grid, but the shape of the grid is constant, having no parameter-effect curvature. (b) Compressed grid. By reparametrizing the model, the grid may become stretched or compressed in regions of the manifold. (c) Rotating, compressed grid. Another parametrization may not only stretch the grid, but also cause the coordinates to rotate. The parameter-effect curvature describes the degree to which the coordinates are stretching and rotating on the manifold. With more than two parameters, there is also a torsion parameter-effect curvature (twisting).

Figure 18

(a) Cross sections of a summing exponential model projected into the plane spanned by the two principal components in data space. Notice the widths of successive cross sections are progressively more narrow, while the deviations from flatness are uniformly spread across the width. The magnitude of the deviation from flatness is approximately the same for each width, giving rise to the hierarchy of curvatures. (b) Cross sections of a feed-forward neural network has many of the same properties as the exponential model. In both cases, the curvature is much smaller than it appears due to the relative scale of the two axes. In fact, the sloppiest directions (narrowest widths) have an aspect ratio of about 1.

Figure 19

The extrinsic and parameter-effect curvature on the model manifold are strongly anisotropic, with the largest curvatures along the shortest widths (see Figs. 6 and 18). The slopes of the (inverse) curvature vs eigenvalue lines are roughly twice that of the singular values (which are equivalent to the widths). The magnitude of the extrinsic curvature is set by the most narrow cross sections, while the magnitude of the parameter-effect curvature is set by the widest cross section. Consequently, the parameter-effect curvature is much larger than the extrinsic curvature. Here we plot the widths and curvatures for a model of four exponentials (above) from Ref. [22] and a feed forward artificial neural network (below).

Figure 20

Curvature anisotropy. (a) Inverse metric eigenvalues. The (inverse) metric has eigenvalues spread over several orders of magnitude, producing a strong anisotropy in the way distances are measured on the model manifold. (b) Geodesic curvature in eigendirections of the metric. The geodesic curvatures also cover many decades. The shortened distance measurements from the metric eigenvalues magnify the anharmonicities in the sloppy directions. (c) Parameter-effect geodesic curvature. The parameter-effect curvature is much larger than the extrinsic curvature, but shares the anisotropy. (d) The eigenvalues of the shape operator. The strong curvature anisotropy described by the geodesic curvature is also illustrated in the eigenvalue spectrum of the shape operator. (e) Parameter-effect shape operator eigenvalues. Two measures (geodesic and shape operator curvatures) span similar ranges, but in both cases the parameter-effect curvature is a factor of about ${10}^5$ larger than the extrinsic curvature equivalent.

Figure 21

A caricature of the widths and curvatures of a typical sloppy model. (a) The manifold deviates by an amount ${\Delta }^N$ from a linear model for each width. As each width is smaller than the last by a factor of $\Delta$, the curvature is largest along the narrow widths. This summary agrees well with the two real models in Fig. 18. (b) The scales of the extrinsic and parameter-effect curvature are set by the narrowest and widest widths, respectively. The parameter-effect curvature is, therefore, smaller than the extrinsic curvature by a factor of ${\Delta }^N$. Both are strongly anisotropic. Compare this figure to the corresponding result for the two real models in Fig. 19.

Figure 22

Model graph curvature. As the Levenberg-Marquardt parameter $\lambda$ is increased, directions with highest curvature become less curved. For stiff directions with less extrinsic curvature, the parameter-effect curvature may be transformed into extrinsic curvature. The damping term reduces the large anisotropy in the curvature. For sufficiently large values of the Levenberg-Marquardt parameters, all curvatures vanish.

Figure 23

(a) Curvature and step size for $\kappa <0$. If $\kappa <0$, then the nonlinearities in the proposed direction are diverting the algorithm away from the desired path. Distances are limited by the size of the curvature. (b) Curvature and step size for $\kappa >0$. The nonlinearities may be helpful to an algorithm, allowing for larger than expected step sizes when $\kappa >0$. (c) Curvature and step size for $\kappa$ with alternating sign. For small $\lambda$, $\kappa <0$ and the nonlinearities are restricting the step size. However, if $\kappa$ becomes positive (the cusp indicates the change of sign), the possible step size suddenly increases. (d) Cost contours for positive and negative values of $\kappa$. One can understand the two different signs of $\kappa$ in terms of which side of the canyon the given point resides. The upper point has positive $\kappa$ and can step much larger distances in the Gauss-Newton direction than can the lower point with negative $\kappa$, which quickly runs up the canyon wall.

Figure 24

The curvature along the width of a manifold effects if the best fit lies on the boundary or on the interior. For a cross-sectional width (thick black line), consider three possibilities: (a) extrinsically flat, (b) constant curvature along width, and (c) curvature proportional to distance from the boundary. Gray regions correspond to data points with best fits on the interior of the manifold, while white regions correspond to data with evaporated parameters. If the curvature is larger near the boundaries, there is less data space available for evaporated best-fit parameters.

Figure 25

Extended geodesic coordinates for the Rosenbrock function. The residuals are one choice of extended geodesic coordinates if the number of parameters equals the number of data points, as is the case for the Rosenbrock function. Because the Rosenbrock function is a simple quadratic, the coordinate transformation can be expressed analytically. Lines of constant $\rho$ are equicost lines, while lines of constant $\varphi$ are the paths a geodesic algorithm should follow to the best fit. Because the geodesics follow the path of the narrow canyon, the radial geodesics are nearly parallel to the equicost lines in parameter space. This effect is actually much more extreme than it appears in this figure because of the relative scales of the two axes.

Figure 26

Greedy step and delayed gratification step criterion. In optimization problems for which there is a long narrow canyon, such as for the Rosenbrock function, choosing a delayed gratification step is important to optimize convergence. By varying the damping term, the algorithm may choose from several possible steps. A greedy step will lower the cost as much as possible, but by doing so it will limit the size of future steps. An algorithm that takes the largest allowable step size (without moving uphill) will not decrease the cost much initially, but will arrive at the best fit in fewer steps and more closely approximate the true geodesic path. What constitutes the largest tolerable step size should be optimized for specific problems so as to guarantee convergence.

Figure 27

Geodesic acceleration in the Rosenbrock Valley. The Gauss-Newton direction, or velocity vector, gives the correct direction that one should move to approach the best fit while navigating a canyon. However, that direction quickly rotates, requiring an algorithm to take very small steps to avoid uphill moves. The geodesic acceleration indicates the direction in which the velocity rotates. The geodesic acceleration determines a parabolic trajectory that can efficiently navigate the valley without running up the wall. The linear trajectory quickly runs up the side of the canyon wall.

Figure 28

(a) Deacceleration when overstepping. Typically the velocity vector greatly overestimates the proper step size. (We have rescaled both velocity and acceleration to fit in the figure.) Algebraically, this is due to the factor of the inverse metric in the velocity, which has very large eigenvalues. The acceleration compensates for this by pointing antiparallel to the velocity. However, the acceleration vector is also very large, as it is multiplied twice by the velocity vector and once by the inverse metric. To make effective use of the acceleration, it is necessary to regularize the metric by including a damping term. (b) Acceleration indicating the direction of the canyon. As the Levenberg-Marquardt parameter is raised, the velocity vector shortens and rotates from the natural gradient into the downhill direction. The acceleration vector also shortens, although much more rapidly, and also rotates. In this two-dimensional cross section, although the two velocity vectors rotate in opposite directions, the accelerations both rotate to indicate the direction that the canyon is turning. By considering the path that one would optimally like to take (along the canyon), it is clear that the acceleration vector is properly indicating the correction to the desired trajectory.

Figure 29

Generalized Rosenbrock results for Levenberg-Marquardt variants. If the canyon that an algorithm must follow is very narrow (measured by the condition number of the metric at the best fit) or turns sharply, the algorithm will require more steps to arrive at the best fit. Those that use the geodesic acceleration term converge more quickly as the canyon narrows. As the parameter-effect curvature increases, the canyon becomes more curved and the problem is more difficult. Notice that changing the canyon’s path from a cubic function in (a) to a quartic function in (b) slowed the convergence rate by a factor of 5. We have omitted the quadratic path since including the acceleration allows the algorithm to find the best fit in one step, regardless of how narrow the canyon becomes.