Algorithmic Differentiation of Numerical Methods: Second-order Adjoint Solvers for Parameterized Systems of Nonlinear Equations

AbstractAdjoint mode algorithmic (also know as automatic) differentiation (AD) transforms implementations of multivariate vector functions as computer programs into first-order adjoint code. Its reapplication or combinations with tangent mode AD yields higher-order adjoint code. Second derivatives play an important role in nonlinear programming. For example, second-order (Newton-type) nonlinear optimization methods promise faster convergence in the neighborhood of the minimum through taking into account second derivative information. The adjoint mode is of particular interest in large-scale gradient-based nonlinear optimization due to the independence of its computational cost on the number of free variables. Part of the objective function may be given implicitly as the solution of a system of n parameterized nonlinear equations. If the system parameters depend on the free variables of the objective, then second derivatives of the nonlinear system's solution with respect to those parameters are required. The local computational overhead as well as the additional memory requirement for the computation of second-order adjoints of the solution vector with respect to the parameters by AD depends on the number of iterations performed by the nonlinear solver. This dependence can be eliminated by taking a symbolic approach to the differentiation of the nonlinear system

Computational Methods for Nonlinear Systems Analysis With Applications in Mathematics and Engineering

An investigation into current methods and new approaches for solving systems of nonlinear equations was performed. Nontraditional methods for implementing arc-length type solvers were developed in search of a more robust capability for solving general systems of nonlinear algebraic equations. Processes for construction of parameterized curves representing the many possible solutions to systems of equations versus finding single or point solutions were established. A procedure based on these methods was then developed to identify static equilibrium states for solutions to multi-body-dynamic systems. This methodology provided for a pictorial of the overall solution to a given system, which demonstrated the possibility of multiple candidate equilibrium states for which a procedure for selection of the proper state was proposed. Arc-length solvers were found to identify and more readily trace solution curves as compared to other solvers making such an approach practical. Comparison of proposed methods was made to existing methods found in the literature and commercial software with favorable results. Finally, means for parallel processing of the Jacobian matrix inherent to the arc-length and other nonlinear solvers were investigated, and an efficient approach for implementation was identified. Several case studies were performed to substantiate results. Commercial software was also used in some instances for additional results verification

Riemann solvers and undercompressive shocks of convex FPU chains

We consider FPU-type atomic chains with general convex potentials. The naive continuum limit in the hyperbolic space-time scaling is the p-system of mass and momentum conservation. We systematically compare Riemann solutions to the p-system with numerical solutions to discrete Riemann problems in FPU chains, and argue that the latter can be described by modified p-system Riemann solvers. We allow the flux to have a turning point, and observe a third type of elementary wave (conservative shocks) in the atomistic simulations. These waves are heteroclinic travelling waves and correspond to non-classical, undercompressive shocks of the p-system. We analyse such shocks for fluxes with one or more turning points. Depending on the convexity properties of the flux we propose FPU-Riemann solvers. Our numerical simulations confirm that Lax-shocks are replaced by so called dispersive shocks. For convex-concave flux we provide numerical evidence that convex FPU chains follow the p-system in generating conservative shocks that are supersonic. For concave-convex flux, however, the conservative shocks of the p-system are subsonic and do not appear in FPU-Riemann solutions

Symbolic computation of exact solutions expressible in hyperbolic and elliptic functions for nonlinear PDEs

Algorithms are presented for the tanh- and sech-methods, which lead to closed-form solutions of nonlinear ordinary and partial differential equations (ODEs and PDEs). New algorithms are given to find exact polynomial solutions of ODEs and PDEs in terms of Jacobi's elliptic functions. For systems with parameters, the algorithms determine the conditions on the parameters so that the differential equations admit polynomial solutions in tanh, sech, combinations thereof, Jacobi's sn or cn functions. Examples illustrate key steps of the algorithms. The new algorithms are implemented in Mathematica. The package DDESpecialSolutions.m can be used to automatically compute new special solutions of nonlinear PDEs. Use of the package, implementation issues, scope, limitations, and future extensions of the software are addressed. A survey is given of related algorithms and symbolic software to compute exact solutions of nonlinear differential equations.Comment: 39 pages. Software available from Willy Hereman's home page at http://www.mines.edu/fs_home/whereman

Interpolatory methods for H∞\mathcal{H}_\infty model reduction of multi-input/multi-output systems

We develop here a computationally effective approach for producing high-quality $\mathcal{H}_\infty$-approximations to large scale linear dynamical systems having multiple inputs and multiple outputs (MIMO). We extend an approach for $\mathcal{H}_\infty$ model reduction introduced by Flagg, Beattie, and Gugercin for the single-input/single-output (SISO) setting, which combined ideas originating in interpolatory $\mathcal{H}_2$-optimal model reduction with complex Chebyshev approximation. Retaining this framework, our approach to the MIMO problem has its principal computational cost dominated by (sparse) linear solves, and so it can remain an effective strategy in many large-scale settings. We are able to avoid computationally demanding $\mathcal{H}_\infty$ norm calculations that are normally required to monitor progress within each optimization cycle through the use of "data-driven" rational approximations that are built upon previously computed function samples. Numerical examples are included that illustrate our approach. We produce high fidelity reduced models having consistently better $\mathcal{H}_\infty$ performance than models produced via balanced truncation; these models often are as good as (and occasionally better than) models produced using optimal Hankel norm approximation as well. In all cases considered, the method described here produces reduced models at far lower cost than is possible with either balanced truncation or optimal Hankel norm approximation

Flexible Differentiable Optimization via Model Transformations

We introduce DiffOpt.jl, a Julia library to differentiate through the solution of optimization problems with respect to arbitrary parameters present in the objective and/or constraints. The library builds upon MathOptInterface, thus leveraging the rich ecosystem of solvers and composing well with modeling languages like JuMP. DiffOpt offers both forward and reverse differentiation modes, enabling multiple use cases from hyperparameter optimization to backpropagation and sensitivity analysis, bridging constrained optimization with end-to-end differentiable programming. DiffOpt is built on two known rules for differentiating quadratic programming and conic programming standard forms. However, thanks ability to differentiate through model transformation, the user is not limited to these forms and can differentiate with respect to the parameters of any model that can be reformulated into these standard forms. This notably includes programs mixing affine conic constraints and convex quadratic constraints or objective function