{'fraction_non_alphanumeric': 0.03494854202401372, 'fraction_numerical': 0.008104631217838765, 'mean_word_length': 4.659063334142198, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Geodesic acceleration and the small-curvature approximation for nonlinear least squares

It has been shown numerically that the performance of the Levenberg-Marquardt algorithm can be improved by including a second order correction known as the geodesic acceleration. In this paper we give the method a more sound theoretical foundation by deriving the geodesic acceleration correction without using differential geometry and showing that the traditional convergence proofs can be adapted to incorporate geodesic acceleration. Unlike other methods which include second derivative information, the geodesic acceleration does not attempt to improve the Gauss-Newton approximate Hessian, but rather is an extension of the small-residual approximation to cubic order. In deriving geodesic acceleration, we note that the small-residual approximation is complemented by a small-curvature approximation. This latter approximation provides a much broader justification for the Gauss-Newton approximate Hessian and Levenberg-Marquardt algorithm. In particular, it is justifiable even if the best fit residuals are large, is dependent only on the model and not on the data being fit, and is applicable for the entire course of the algorithm and not just the region near the minimum.


I. INTRODUCTION
In this paper we consider the problem of minimizing a scalar function whose form is the sum of squares where r : R N → R M is an M-dimensional vector function of N parameters, θ. We refer to r(θ) as the residuals and scalar function C(θ) as the cost. Functions of this form often arise in the context of data fitting and represent an important class of problems as evidenced by the large number of software packages dedicated to their optimization.
The structure of this particular problem lends itself to efficient optimization. In particular, consider the Hessian matrix of second derivatives, necessary to implement a quasi Newton method: where in the last two lines we have applied the so-called Gauss-Newton or small-residual approximation and introduced the Jacobian matrix of first derivatives J. The approximation in the second line is usually justified by the hope that at a minimum of C(θ), the individual residuals are small, so that the Hessian is dominated by the contributions from the first term.
Numerically, this approximation is advantageous since it allows one to implement a quasi-Newton method by calculating only the first derivatives of the residuals. This approximation comes at the cost of storing the derivative information for each residual individually, but this is rarely a bottleneck on modern computers. Consequently, the functional form of Eq. (1) effectively allows one to estimate the Hessian matrix with the same information used to calculate the gradient.
Applying a trust region method to the Gauss-Newton approximate Hessian results in the Levenberg-Marquardt algorithm [1][2][3][4][5][6] which iteratively updates the parameters according to where λ is an appropriately chosen Langrange multiplier for the step bound δθ T D T Dδθ ≤ ∆ 2 and g = J T r is the gradient. Because Levenberg-Marquardt is a quasi-Newton method, it usually has very good convergence properties. In particular, if the small-residual approximation is good, convergence can be super-linear to a local minimum. Furthermore, for well-constructed choices of λ and ∆, the method is globally convergent [3,5].
Although the Levenberg-Marquardt algorithm has many desirable properties, it is not always ideal. Often data fitting problems have a cost function characterized by narrow, winding canyons. Although, asymptotically the algorithm may converge super-linearly, it may nevertheless spend an unreasonable amount of time navigating the winding canyon before it finally zooms into the minimum. This problem is typically more severe on problems with many parameters, which in turn are often more computationally expensive to evaluate and lack good parameter estimates to use as starting points. An improved optimization method which can find minima with fewer function evaluations (and especially fewer Jacobian evaluations) would be a welcome improvement.
In order to help improve the Levenberg-Marquardt algorithm, the authors previously proposed the inclusion of a geodesic acceleration term in the algorithm [7,8]. This correction was derived using an information geometric interpretation of the least-square problem and was justified based on the empirically observed small-extrinsic curvatures of the relevant manifolds. In this paper, we will see that the geodesic acceleration correction can be understood as a generalization of the Gauss-Newton method extended to cubic order. Although other methods exist which utilize higher-order information, geodesic acceleration is complementary to these approaches as their primary motivation is to improve the estimate of the Hessian. By contrast, the geodesic acceleration assumes the Hessian estimate is adequate to proceed to higher-order.
In this paper, we derive the geodesic acceleration in a geometric independent way (section II) and prove that its inclusion in the Levenberg-Marquardt algorithm does not compromise its convergence properties (section III). From the explicit derivation in section II, we see that geodesic acceleration can be understood as a continuation of the small-residual approximation to higher-order terms. Indeed, the discarded terms are properly understood as the residuals coupled to the extrinsic curvatures of the Model Graph defined in references [7,8]. The small-curvature approximation, therefore, provides additional justification for the Gauss-Newton approximation and the geodesic acceleration correction. As we argue in section IV, the small-curvature approximation is more useful than the small-residual approximation for many reasons. In particular, it is valid not only near a local minimum of the cost, but for all parameter values; it is a property of the model and not the data being fit and so is valid even when the model cannot fit the data well; furthermore, numerical experiments on many models suggest the small-curvature approximation is nearly universally valid.

II. DERIVATION
In order to improve the efficiency of the Levenberg-Marquardt method, we propose modifying the step to include higher order corrections in a numerically efficient manner. To derive this correction, consider the minimization problem of finding the best residuals with a constrained step-size. We write the dependence of the residual on the shift δθ as where J and K are the arrays of first and second derivatives respectively. We wish to minimize min δθ with the constraint that δθ T D T Dδθ ≤ ∆ 2 . After introducing a Lagrange multiplier λ for the constraint in the step size, the minimization becomes min δθ r + Jδθ + 1/2 δθ T Kδθ 2 + λδθ T D T Dδθ.
By varying δθ we find the normal equations: where we have explicitly included all the indices to avoid any ambiguity. Since we constrain the step size, it is natural to assume that δθ is small, and we seek a solution of Eq. (7) as a perturbation series around the linearized equation: Let δθ 1 be a solution of the linearized equation: where in the second line we have made the usual Gauss-Newton approximation. We do not actually discard the term involving K, as it will help to cancel out a higher order correction later in the derivation. We therefore set which is the usual Levenberg-Marquardt step.
With this definition of δθ 1 , Eq. (7) becomes to second order, with the term r m K mµα δθ 1α the term neglected by the Gauss-Newton approximation at first order.
We now turn our attention to the second term in parentheses in Eq. (10). Using the Since this term is proportional to the residuals, r m , it can be ignored using the usual smallresidual arguments. However, we now make an appeal to geometric considerations by noting that δ mn − βν J mβ (J T J + λD T D) −1 βν J nν = P N mn is a matrix that projects vectors perpendicular to the tangent plane of the Model Graph as described in reference [8]. If the curvature of the model graph is small, then P N K ≈ 0 and this term can be neglected. We discuss the implications of this argument further in section IV.
Returning to Eq. (10), after ignoring the last term in parentheses, we find where we have introduced the directional second derivative r ′′ m = µν K mµν δθ 1µ δθ 1ν and in the second line made the usual Gauss-Newton approximation to the Hessian, giving the formula first presented in [7]. This formula was originally interpreted as the second order correction to geodesic flow on the model graph, and so we refer to it as the geodesic acceleration correction. By analogy, we refer to the first order term as the geodesic velocity.
The full step is therefore given by: As has been noted previously [7,8], although the geodesic acceleration correction includes second derivative information at each step of the algorithm, its calculation is not computationally intensive. In particular, it only requires the evaluation of a directional second derivative, which is computationally comparable to a single evaluation of the residuals. Indeed, in the absence of an analytic expression, a finite difference estimate of the relevant second derivative can be found by a single function evaluation. In contrast, the Jacobian evaluation at each step is comparable to N function evaluations. Particularly for large problems, the computational cost of including the geodesic acceleration is negligible compared the other elements of the algorithm.

III. CONVERGENCE
In order to show that geodesic acceleration does not impair the convergence guarantees of the Levenberg-Marquardt algorithm, we must make a few additional modifications. We first note that an algorithm that selects ∆ directly requires a solution of Eq. (5) given the step bound, i.e. find the value of λ corresponding to the step bound ∆. This so-called subproblem, can be solved accurately and efficiently for the case K = 0, using the methods described by Moré[9] and Nocedel and Wright [5]. When including geodesic acceleration however, such a simple solution does not exist. Indeed, the step size is no longer a monotonically decreasing function of λ. Furthermore, accounting for the contribution from the second term requires additional function evaluations, making an accurate solution computationally expensive.
Fortunately, convergence proofs for Levenberg-Marquardt do not require that this subproblem be solved accurately. We therefore content ourselves by approximately solving the problem as follows: We first require that |δθ 1 | ≤ ∆. This step can be satisfied easily using the algorithm in references [5,9]. We next require that the relative contribution from the second order step be bounded for some α > 0 [17]. In practice we implement the requirement in Eq. (13) by rejecting all proposed steps for which it is not satisfied and decreasing ∆ (or increasing λ) until an acceptable step is generated, recalculating δθ 1 for the new ∆. This may appear inefficient since the method will on occasion reject steps that would have decreased the cost. However, this requirement adds stability to the algorithm, helping it to avoid undesirable fixed points with infinite parameter values as argued and numerically justified in references [7,8,10].
We need to make an additional assumption about the behavior of the model. Specifically, it is necessary to assume that the directional second derivative of the model is bounded: | µν K mµν u µ u ν | < κ for any parameter-space unit vector u and some positive constant κ.
This assumption is necessary to guarantee that we can always find a step-size that satisfies Eq. (13). We do not anticipate this requirement to be a major restriction for the applicability of the algorithm as we discuss later in this section.
With this additional assumption, we can show that our modified Levenberg-Marquardt algorithm enjoys the same global convergence properties as original algorithm. Indeed, the proof is nearly identical to that of the unmodified Levenberg-Marquardt. Our proof here follows closely that of Theorem 4.5 in Nocedal and Wright [5]. First we define the model function m k (δθ) by and the reduction ratio ρ k by We now define an algorithm for which we will prove convergence. This algorithm is analogous to algorithm 4.1 in reference [5] which we recover if we were to set δθ 2 = 0 at each step.

Algorithm 1
Given∆ > 0, ∆ 0 ∈ (0,∆), and α > 0 for k = 0, 1, 2, . . . do Calculate δθ 1 and δθ 2 as described above Before presenting our proof, first notice that by definingθ = Dθ, the optimization problem inθ must at each step satisfy the bound |δθ 1 | < ∆. Without loss of generality, we therefore assume that D T D is the identity. This assumption essentially replaces the Jacobian matrix withJ = JD −1 . With this additional assumption, we now present a Lemma that will be useful in proving convergence.
With this Lemma, we are prepared to prove our main result: that including geodesic acceleration does not affect the convergence properties of the Levenberg-Marquardt algorithm.

Theorem 1
Suppose that |J T J| ≤ β for some positive constant β, that C is Lipschitz continuously differentiable in the neighborhood S(R 0 ) for some R 0 > 0, and that |δθ 1 | < ζ∆ k for some constant ζ > 1 at each iteration. Also assume that m k (0) − m k (δθ) ≥ c 1 |g k | min(∆ k , |g k |/|J T J|) for some positive constant c 1 ∈ (0, 1] (g k = J T k r k ), and that the second directional derivative of r k in any unit direction u is bounded | µν K mµν u µ u ν | < κ for some positive κ. Then using Algorithm 1, lim k→∞ inf |g k | = 0. (20)

Proof
The proof is nearly identical to that of Theorem 4.5 in reference [5]. We repeat the proof here in order to highlight the small-differences, leaving out the algebraic details.
With some algebraic manipulation, we obtain Using the same argument as in reference [5], we have where β 1 is the Lipschitz constant for g on the set S(R 0 ).
Using Eqs. (24), (23), and (21), we have We now derive a bound on the right-hand-side that holds for all sufficiently small-values of ∆ k , that is, for all ∆ k ≤∆, where∆ is defined as follows: As noted in reference [5], the R 0 /γ term in this definition ensures that the bound in Eq. (21) is valid because |δθ k | ≤ γ∆ k ≤ γ∆ ≤ R 0 . The third term is necessary to show convergence when including geodesic acceleration.
Since, we have that if ∆ k ∈ [0,∆], then |δθ 2 | < α|δθ 1 |/2 and that ρ k > 1 4 , by the workings of Algorithm 1, ∆ k+1 ≥ ∆ k whenever ∆ k falls below the threshold∆. Consequently, a reduction in ∆ k (by a factor of 1/4) can occur in the algorithm only if ∆ k ≥∆. We conclude that ∆ k+1 ≥ min(∆ k ,∆/4) for all k ≥ K. (27) Suppose now that there is an infinite subsequence K such that ρ k ≥ 1/4 for k ∈ K. For k ∈ K and k ≥ K, we have from Eq. (23) that Since C is bounded below, it follows from this inequality that lim k∈K,k→∞ contradicting Eq. (27). Hence no such infinite subsequence K can exist, and we must have ρ k < 1 4 for all k sufficiently large. In this case, ∆ k will eventually be multiplied by 1 4 at every iteration, and we have lim k→∞ ∆ k = 0, which again contradicts Eq. (27). Hence our original assertion, Eq. (22) must be false, proving the theorem.
Although Theorem 1 and its subsequent proof are nearly identical to Theorem 4.5 in Nocedal and Wright [5], we now briefly discuss its implication. In particular, the proof is that the algorithm converges to a point where cost has zero gradient. This does not automatically imply convergence in the parameters unless the Hessian matrix is bounded from below. In fact, the Hessian matrix can be very ill-conditioned, particularly for models with many parameters [11,12]. Although the inferred parameters in such problems may be ill-conditioned, the convergence properties of the algorithm are nevertheless robust.
The requirement that m k (0) − m k (δθ) ≥ c 1 |g k | min(∆ k , |g k |/|J T J|) at each step of the algorithm in this context is given without motivation. However, just like the analogous theorem in reference [5], this requirement is necessary to guarantee that the parameter iterates do not accumulate at points with nonzero gradient. We note that in practice this condition is never explicitly checked.
The only new assumption of Theorem 1 beyond those for the standard algorithms [5] is the bound on the second derivative. Since we are utilizing second derivative information, this addition is not unexpected and guarantees that the new algorithm is well-behaved. It is also clearly analogous to the bound on the first derivative used in the first algorithm (|J T J| <= β, equivalent to |Ju| <= √ β). The standard bound on the first derivative already excludes models where the cost diverges at interior points in the domain of the model as well as points like C(θ) ∼ |θ − θ 0 | where the cost would stay finite although the derivative diverges; the new bound on the second derivative additionally excludes cases like C(θ) ∼ |θ − θ 0 | 3/2 , which should not arise often in practice. In our experience, many models do have costs which become singular at unphysical points (i.e., positive-only parameters set less than zero), which can be addressed by an appropriate reparameterization of the model (i.e., shifting to log parameters). When this can be done, it is also likely to improve the convergence rate of the algorithm.
We have here shown convergence of geodesic acceleration for an algorithm that belongs to the broad class of trust region methods that operate by specifying a step bound ∆.
Originally Levenberg-Marquardt was proposed as an algorithm that heuristically selected λ directly rather than implicitly through ∆ [1,2]. For an algorithm that directly selects λ, one can similarly show convergence [3]. The proof that geodesic acceleration does not impair the convergence of this class of algorithm is almost identical to that of the original theorem, just as the proof of Theorem 1 above is almost identical to that of Theorem 4.5 in reference [5].
We do not give the rigorous proof here, but merely note that under the same additional assumptions, convergence can be shown for these methods as well.

METHODS
In deriving the geodesic acceleration in section II, we observed that the neglected terms were each proportional to the residuals and could be neglected in the small-residual approximation. We also noted that the neglected terms were also proportional to the extrinsic curvature on the model graph as described in reference [7,8]. Indeed, the geodesic acceleration was originally understood as a small-curvature approximation rather than a small-residual approximation. These two approximations are complementary; in general only one of the two needs to be valid in order for the approximations to hold.
In deriving the geodesic acceleration, the connection between the small-residual and smallcurvature approximation only became apparent when considering cubic order terms. How-ever, the equivalence of the two approximations can be seen in the neglected terms of the Gauss-Newton Hessian without considering the geodesic acceleration. At a fixed point of the cost, the residuals are perpendicular to the model manifold (a fact used elsewhere to justify a scale-free measure of convergence [13]), i.e. r n ≈ m r m P N mn . Thus, as the algorithm approaches a minimum, the neglected term can be written as As Eq. (28) makes clear, near the best fit the nonlinear contribution to the Hessian includes only components perpendicular to the tangent plane and is thus proportional to the extrinsic curvature.
The implication of Eq. (28) is that the Levenberg-Marquardt algorithm may attain very good convergence rates even with large residuals so long as the extrinsic curvature is sufficiently small. However, the approximation in Eq. (28) is only valid near a fixed point of the cost, just as is the small-residual approximation is only valid near the minimum. The advantage of identifying the equivalence of the small-residual and small-curvature approximations is that the latter is likely to have much wider applicability. In particular for data fitting, the small-curvature approximation is a feature of the model rather than the data. Thus, although the validity of the small-residual approximation cannot be identified without finding the best fit, the small-curvature approximation can be. Furthermore, the small-curvature approximation is likely to be an excellent approximation for most models. In particular, several examples in references [8,14] exhibit extrinsic curvatures many orders of magnitude smaller than the magnitude of the bare nonlinearities.
Although Eq. (28) is only valid near a minimum, the approximation that led to the geodesic acceleration is valid over a much larger parameter range, and it is for this reason it is likely to be a useful modification. The problem with Levenberg-Marquardt is not that its asymptotic convergence rate is poor (super-linear convergence is satisfactory in most cases), but that it may spend an unreasonable amount of time navigating a narrow canyon before that convergence rate is realized. As originally suggested, the geodesic acceleration can speed up this process since it describes the curvature of the canyon. The algorithm can find the minimum more quickly by following the canyon with a sequence of parabolic steps, rather than linear steps. A numerical demonstration of this will not be given here, as it has already been shown elsewhere [7,8,10].
There are many other methods which include second derivative information in order to improve the Levenberg-Marquardt algorithm. These approaches all try to improve the quality of the Hessian estimate in some way [5,15,16]. Thus, although geodesic acceleration may appear superficially similar to these approaches, it is actually quite different. The philosophy behind geodesic acceleration is to extend the small-residual/curvature approximation to higher order terms rather than estimate the neglected terms. Of course, the utility of such of an approach is ultimately measured by performance on real world problems. As has been shown elsewhere, geodesic acceleration can be very helpful on large problems for which Levenberg-Marquardt spends unreasonable time searching for a minimum before zooming into the best fit.

V. CONCLUSION
In this paper we have studied the geodesic acceleration correction to the Levenberg-Marquardt algorithm, originally suggested in references [7,8]. We have derived the correction without the use of differential geometry and shown that it can be interpreted as an extension of the small-residual approximation used to estimate the Hessian matrix of a sum of squares. We have also seen that the small-residual approximation is complemented by the small-curvature approximation, which is likely to be applicable under much more general circumstances.
We have seen that with just a few minor modifications, the geodesic acceleration correction can be incorporated into a standard Levenberg-Marquardt routine with minimal computational cost and without sacrificing its convergence properties. Numerical experiments given elsewhere [7,8,10], suggest that the benefit/cost ratio of this improvement can be substantial on many optimization problems.