The Adjoint Newton Algorithm for Large-Scale Unconstrained Optimization in Meteorology Applications

A new algorithm is presented for carrying out large-scale unconstrained optimization required in variational data assimilation using the Newton method. The algorithm is referred to as the adjoint Newton algorithm. The adjoint Newton algorithm is based on the first- and second-order adjoint techniques allowing us to obtain the Newton line search direction by integrating a tangent linear equations model backwards in time (starting from a final condition with negative time steps). The error present in approximating the Hessian (the matrix of second-order derivatives) of the cost function with respect to the control variables in the quasi-Newton type algorithm is thus completely eliminated, while the storage problem related to the Hessian no longer exists since the explicit Hessian is not required in this algorithm. The adjoint Newton algorithm is applied to three one-dimensional models and to a two-dimensional limited-area shallow water equations model with both model generated and First Global Geophysical Experiment data. We compare the performance of the adjoint Newton algorithm with that of truncated Newton, adjoint truncated Newton, and LBFGS methods. Our numerical tests indicate that the adjoint Newton algorithm is very efficient and could find the minima within three or four iterations for problems tested here. In the case of the two-dimensional shallow water equations model, the adjoint Newton algorithm improves upon the efficiencies of the truncated Newton and LBFGS methods by a factor of at least 14 in terms of the CPU time required to satisfy the same convergence criterion.The Newton, truncated Newton and LBFGS methods are general purpose unconstrained minimization methods. The adjoint Newton algorithm is only useful for optimal control problems where the model equations serve as strong constraints and their corresponding tangent linear model may be integrated backwards in time. When the backwards integration of the tangent linear model is ill-posed in the sense of Hadamard, the adjoint Newton algorithm may not work. Thus, the adjoint Newton algorithm must be used with some caution. A possible solution to avoid the current weakness of the adjoint Newton algorithm is proposed.     

Truncated Newton method for sparse unconstrained optimization using automatic differentiation

When solving large complex optimization problems, the user is faced with three major problems. These are (i) the cost in human time in obtaining accurate expressions for the derivatives involved; (ii) the need to store second derivative information; and (iii), of lessening importance, the time taken to solve the problem on the computer. For many problems, a significant part of the latter can be attributed to solving Newton-like equations. In the algorithm described, the equations are solved using a conjugate direction method that only needs the Hessian at the current point when it is multiplied by a trial vector. In this paper, we present a method that finds this product using automatic differentiation while only requiring vector storage. The method takes advantage of any sparsity in the Hessian matrix and computes exact derivatives. It avoids the complexity of symbolic differentiation, the inaccuracy of numerical differentiation, the labor of finding analytic derivatives, and the need for matrix store. When far from a minimum, an accurate solution to the Newton equations is not justified, so an approximate solution is obtained by using a version of Dembo and Steihaug's truncated Newton algorithm (Ref. 1).This paper was presented at the SIAM National Meeting, Boston, Massachusetts, 1986.     

Sparse quadratic programming in chemical process optimization

The quadratic programming aspects of a full space successive quadratic programming (SQP) method are described. In particular, fill-in, matrix factor and active set updating, numerical stability, and indefiniteness of the Hessian matrix are discussed in conjunction with a sparse modification of Bunch and Parlett factorization of symmetric indefinite (Kuhn-Tucker) matrices of the type often encountered in optimization. A new pivoting strategy, called constrained pivoting, is proposed to reduce fill-in and compared with complete, partial and threshold pivoting. It is shown that constrained pivoting often significantly reduces fill-in and thus the iterative computational burdens associated with the factorization and solution of Kuhn-Tucker conditions within the QP subproblem. A numerical algorithm for updating the lower triangular and diagonal factors is presented and shown to be very fast, usually requiring only about 5% of the cost of refactorization. Two active set strategies are also presented. These include the options of adding inequalities either one or several at a time. In either case, the effects on matrix factor updating is shown to be small. Finally, a simple test is used to maintain iterative descent directions in the quadratic program. Our sparse symmetric indefinite QP algorithm is tested in the context of a family of SQP algorithms that include a full space Newton method with analytical derivatives, a full space BFGS method and a Range and Null space Decomposition (RND) method in which the projected Hessian is calculated from either analytical second derivatives or the BFGS update. Several chemical process optimization problems, with small and large degrees of freedom, are used as test problems. These include minimum work calculations for multistage isothermal compression, minimum area targeting for heat exchanger networks, and distillation optimizations involving some azeotropic and extractive distillations. Numerical results show uniformly that both the proposed QP and SQP algorithms, particularly the full space Newton method, are reliable and efficient. No failures were experienced at either level.     

A priori estimate for convex solutions to special Lagrangian equations and its application

We derive a priori interior Hessian estimates for special Lagrangian equations when the potential is convex. When the phase is very large, we show that continuous viscosity solutions are smooth in the interior of the domain. © 2008 Wiley Periodicals, Inc.     

Validation of an Augmented Lagrangian Algorithm with a Gauss-Newton Hessian Approximation Using a Set of Hard-Spheres Problems

An Augmented Lagrangian algorithm that uses Gauss-Newton approximations of the Hessian at each inner iteration is introduced and tested using a family of Hard-Spheres problems. The Gauss-Newton model convexifies the quadratic approximations of the Augmented Lagrangian function thus increasing the efficiency of the iterative quadratic solver. The resulting method is considerably more efficient than the corresponding algorithm that uses true Hessians. A comparative study using the well-known package LANCELOT is presented.     

Global Convergence Technique for the Newton Method with Periodic Hessian Evaluation

The problem of globalizing the Newton method when the actual Hessian matrix is not used at every iteration is considered. A stabilization technique is studied that employs a new line search strategy for ensuring the global convergence under mild assumptions. Moreover, an implementable algorithmic scheme is proposed, where the evaluation of the second derivatives is conditioned to the behavior of the algorithm during the minimization process and the local convexity properties of the objective function. This is done in order to obtain a significant computational saving, while keeping acceptable the unavoidable degradation in convergence speed. The numerical results reported indicate that the method described may be employed advantageously in all applications where the computation of the Hessian matrix is highly time consuming.     

Using the KKT matrix in an augmented Lagrangian SQP method for sparse constrained optimization

The augmented Lagrangian SQP subroutine OPALQP was originally designed for small-to-medium sized constrained optimization problems in which the main calculation on each iteration, the solution of a quadratic program, involves dense, rather than sparse, matrices. In this paper, we consider some reformulations of OPALQP which are better able to take advantage of sparsity in the objective function and constraints.The modified versions of OPALQP differ from the original in using sparse data structures for the Jacobian matrix of constraints and in replacing the dense quasi-Newton estimate of the inverse Hessian of the Lagrangian by a sparse approximation to the Hessian. We consider a very simple sparse update for estimating ² L and also investigate the benefits of using exact second derivatives, noting in the latter case that safeguards are needed to ensure that a suitable search direction is obtained when ² L is not positive definite on the null space of the active constraints.The authors are grateful to John Reid and Nick Gould of the Rutherford Appleton Laboratory for a number of helpful and interesting discussions. Thanks are also due to Laurence Dixon for comments which led to the clarification of some parts of the paper.This work has been partly supported by a CAPES Research Studentship funded by the Brazilian Government.     

Optimization of unconstrained functions with sparse hessian matrices-newton-type methods

Newton-type methods for unconstrained optimization problems have been very successful when coupled with a modified Cholesky factorization to take into account the possible lack of positive-definiteness in the Hessian matrix. In this paper we discuss the application of these method to large problems that have a sparse Hessian matrix whose sparsity is known a priori. Quite often it is difficult, if not impossible, to obtain an analytic representation of the Hessian matrix. Determining the Hessian matrix by the standard method of finite-differences is costly in terms of gradient evaluations for large problems. Automatic procedures that reduce the number of gradient evaluations by exploiting sparsity are examined and a new procedure is suggested. Once a sparse approximation to the Hessian matrix has been obtained, there still remains the problem of solving a sparse linear system of equations at each iteration. A modified Cholesky factorization can be used. However, many additional nonzeros (fill-in) may be created in the factors, and storage problems may arise. One way of approaching this problem is to ignore fill-in in a systematic manner. Such technique are calledpartial factorization schemes. Various existing partial factorization are analyzed and three new ones are developed. The above algorithms were tested on a set of problems. The overall conclusions were that these methods perfom well in practice.     

An optimal algorithm for a class of equality constrained quadratic programming problems with bounded spectrum

The implementation of the recently proposed semi-monotonic augmented Lagrangian algorithm for the solution of large convex equality constrained quadratic programming problems is considered. It is proved that if the auxiliary problems are approximately solved by the conjugate gradient method, then the algorithm finds an approximate solution of the class of problems with uniformly bounded spectrum of the Hessian matrix at O(1) matrix–vector multiplications. If applied to the class of problems with the Hessian matrices that are in addition either sufficiently sparse or can be expressed as a product of such sparse matrices, then the cost of the solution is proportional to the dimension of the problems. Theoretical results are illustrated by numerical experiments. This research is supported by grants of the Ministry of Education No. S3086102, ET400300415 and MSM 6198910027.     

The estimation of the hessian matrix in nonlinear least squares problems with non-zero residuals

Computational experiments by McKeown [11] have shown that specialised methods, based on the Gauss—Newton iteration, are not necessarily the best choice for minimising functions that are sums of squared terms. Difficulties arise when the Gauss—Newton approach does not yield a good approximation to the second derivative matrix of the function: and this is more likely to happen when the function value at the optimum is not near zero and the terms in the sum of squares are significantly nonlinear. This paper considers some specialised methods for the nonlinear least squares problem which seek to improve the Gauss—Newton estimate of the Hessian matrix without explicitly calculating second derivatives.     

Singular solutions of Hessian fully nonlinear elliptic equations

We study Hessian fully nonlinear uniformly elliptic equations and show that the second derivatives of viscosity solutions of those equations (in 12 or more dimensions) can blow up in an interior point of the domain. We prove that the optimal interior regularity of such solutions is no more than C1+?, showing the optimality of the known interior regularity result. The same is proven for Isaacs equations. We prove the existence of non-smooth solutions to fully nonlinear Hessian uniformly elliptic equations in 11 dimensions. We study also the possible singularity of solutions of Hessian equations defined in a neighborhood of a point and prove that a homogeneous order 0<α<1 solution of a Hessian uniformly elliptic equation in a punctured ball should be radial.     

Parameter estimation of ordinary differential equations

This paper addresses the development of a new algorithm forparameter estimation of ordinary differential equations. Here,we show that (1) the simultaneous approach combined with orthogonalcyclic reduction can be used to reduce the estimation problemto an optimization problem subject to a fixed number of equalityconstraints without the need for structural information to devisea stable embedding in the case of non-trivial dichotomy and(2) the Newton approximation of the Hessian information of theLagrangian function of the estimation problem should be usedin cases where hypothesized models are incorrect or only a limitedamount of sample data is available. A new algorithm is proposedwhich includes the use of the sequential quadratic programming(SQP) Gauss–Newton approximation but also encompassesthe SQP Newton approximation along with tests of when to usethis approximation. This composite approach relaxes the restrictionson the SQP Gauss–Newton approximation that the hypothesizedmodel should be correct and the sample data set large enough.This new algorithm has been tested on two standard problems.     

An incomplete Hessian Newton minimization method and its application in a chemical database problem

To efficiently solve a large scale unconstrained minimization problem with a dense Hessian matrix, this paper proposes to use an incomplete Hessian matrix to define a new modified Newton method, called the incomplete Hessian Newton method (IHN). A theoretical analysis shows that IHN is convergent globally, and has a linear rate of convergence with a properly selected symmetric, positive definite incomplete Hessian matrix. It also shows that the Wolfe conditions hold in IHN with a line search step length of one. As an important application, an effective IHN and a modified IHN, called the truncated-IHN method (T-IHN), are constructed for solving a large scale chemical database optimal projection mapping problem. T-IHN is shown to work well even with indefinite incomplete Hessian matrices. Numerical results confirm the theoretical results of IHN, and demonstrate the promising potential of T-IHN as an efficient minimization algorithm.     

Computation of sparse low degree interpolating polynomials and their application to derivative-free optimization

Interpolation-based trust-region methods are an important class of algorithms for Derivative-Free Optimization which rely on locally approximating an objective function by quadratic polynomial interpolation models, frequently built from less points than there are basis components. Often, in practical applications, the contribution of the problem variables to the objective function is such that many pairwise correlations between variables are negligible, implying, in the smooth case, a sparse structure in the Hessian matrix. To be able to exploit Hessian sparsity, existing optimization approaches require the knowledge of the sparsity structure. The goal of this paper is to develop and analyze a method where the sparse models are constructed automatically. The sparse recovery theory developed recently in the field of compressed sensing characterizes conditions under which a sparse vector can be accurately recovered from few random measurements. Such a recovery is achieved by minimizing the ℓ ₁-norm of a vector subject to the measurements constraints. We suggest an approach for building sparse quadratic polynomial interpolation models by minimizing the ℓ ₁-norm of the entries of the model Hessian subject to the interpolation conditions. We show that this procedure recovers accurate models when the function Hessian is sparse, using relatively few randomly selected sample points. Motivated by this result, we developed a practical interpolation-based trust-region method using deterministic sample sets and minimum ℓ ₁-norm quadratic models. Our computational results show that the new approach exhibits a promising numerical performance both in the general case and in the sparse one.     

Calibration of the SABR Model in Illiquid Markets

Recently the SABR model has been developed to manage the option smile which is observed in derivatives markets. Typically, calibration of such models is straightforward as there is adequate data available for robust extraction of the parameters required asinputs to the model. The paper considers calibration of the model in situations where input data is very sparse. Although this will require some creative decision making, the algorithms developed here are remarkably robust and can be used confidently for mark to market and hedging of option portfolios.     

Convergence of quasi-Newton matrices generated by the symmetric rank one update

Quasi-Newton algorithms for unconstrained nonlinear minimization generate a sequence of matrices that can be considered as approximations of the objective function second derivatives. This paper gives conditions under which these approximations can be proved to converge globally to the true Hessian matrix, in the case where the Symmetric Rank One update formula is used. The rate of convergence is also examined and proven to be improving with the rate of convergence of the underlying iterates. The theory is confirmed by some numerical experiments that also show the convergence of the Hessian approximations to be substantially slower for other known quasi-Newton formulae.The work of this author was supported by the National Sciences and Engineering Research Council of Canada, and by the Information Technology Research Centre, which is funded by the Province of Ontario.     

Inexact trust region PGC method for large sparse unconstrained optimization

We present an algorithm, partitioning group correction (PGC) algorithm based on trust region and conjugate gradient method, for large-scale sparse unconstrained optimization. In large sparse optimization, computing the whole Hessian matrix and solving the Newton-like equations at each iteration can be considerably expensive when a trust region method is adopted. The method depends on a symmetric consistent partition of the columns of the Hessian matrix and an inaccurate solution to the Newton-like equations by conjugate gradient method. And we allow that the current direction exceeds the trust region bound if it is a good descent direction. Besides, we studies a method dealing with some sparse matrices having a dense structure part. Some good convergence properties are kept and we contrast the computational behavior of our method with that of other algorithms. Our numerical tests show that the algorithm is promising and quite effective, and that its performance is comparable to or better than that of other algorithms available.     

The Hessian and Jacobi morphisms for higher order calculus of variations

We formulate higher order variations of a Lagrangian in the geometric framework of jet prolongations of fibered manifolds. Our formalism applies to Lagrangians which depend on an arbitrary number of independent and dependent variables, together with higher order derivatives. In particular, we show that the second variation is equal (up to horizontal differentials) to the vertical differential of the Euler-Lagrange morphism which turns out to be self-adjoint along solutions of the Euler-Lagrange equations. These two objects, respectively, generalize in an invariant way the Hessian morphism and the Jacobi morphism (which is then self-adjoint along critical sections) of a given Lagrangian to the case of higher order Lagrangians. Some examples of classical Lagrangians are provided to illustrate our method.     

Fast automatic differentiation as applied to the computation of second derivatives of composite functions

A technique for deriving formulas for the second derivatives of a composite function with constrained variables is proposed. The original system of constraint equations is associated with a linear system of equations, whose solution is used to determine the Hessian of the function. The resulting formulas are applied to discrete problems obtained by approximating optimal control problems with the use of Runge-Kutta methods of various orders. For a particular optimal control problem, the numerical results obtained by the gradient method and Newton’s method with the resulting formulas are described and analyzed in detail.     

Sample size selection in optimization methods for machine learning

This paper presents a methodology for using varying sample sizes in batch-type optimization methods for large-scale machine learning problems. The first part of the paper deals with the delicate issue of dynamic sample selection in the evaluation of the function and gradient. We propose a criterion for increasing the sample size based on variance estimates obtained during the computation of a batch gradient. We establish an complexity bound on the total cost of a gradient method. The second part of the paper describes a practical Newton method that uses a smaller sample to compute Hessian vector-products than to evaluate the function and the gradient, and that also employs a dynamic sampling technique. The focus of the paper shifts in the third part of the paper to L ₁-regularized problems designed to produce sparse solutions. We propose a Newton-like method that consists of two phases: a (minimalistic) gradient projection phase that identifies zero variables, and subspace phase that applies a subsampled Hessian Newton iteration in the free variables. Numerical tests on speech recognition problems illustrate the performance of the algorithms.