Search

The implicit Euler method integrates systems of ordinary differential equations $$\frac{d x}{d t}=G(t,x(t))$$ with differentiable right-hand side $G : {\mathbb R} \times {\mathbb R}^n \rightarrow {\mathbb R}^n$ from an initial state $x=x(0) \in {\mathbb R}^n$ to a target time $t \in {\mathbb R}$ as $x(t)=E(t,m,x)$ using an equidistant discretization of the time interval $[0,t]$ yielding $m>0$ time steps. We present a method for efficiently computing the product of its inverse Jacobian $$(E')^{-1} \equiv \left (\frac{d E}{d x}\right )^{-1} \in {\mathbb R}^{n \times n} $$ with a given vector $v \in {\mathbb R}^n.$ We show that the differential inverse $(E')^{-1} \cdot v$ can be evaluated for given $v \in {\mathbb R}^n$ with a computational cost of $\mathcal{O}(m \cdot n^2)$ as opposed to the standard $\mathcal{O}(m \cdot n^3)$ or, naively, even $\mathcal{O}(m \cdot n^4).$ The theoretical results are supported by actual run times. A reference implementation is provided., Comment: 12 pages, 3 figures

The efficient computation of Jacobians represents a fundamental challenge in computational science and engineering. Large-scale modular numerical simulation programs can be regarded as sequences of evaluations of in our case differentiable subprograms with corresponding elemental Jacobians. The latter are typically not available. Tangent and adjoint versions of the individual subprograms are assumed to be given as results of algorithmic differentiation instead. The classical (Jacobian) Matrix Chain Product problem is reformulated in terms of matrix-free Jacobian-matrix (tangents) and matrix-Jacobian products (adjoints), subject to limited memory for storing information required by latter. All numerical results can be reproduced using an open-source reference implementation., Comment: arXiv admin note: substantial text overlap with arXiv:2003.05755

We use Algorithmic Differentiation (AD) to implement type-generic tangent and adjoint versions of $$ y=\sum_{i=0}^{n-1} x_{2 i} \cdot x_{2 i+1} $$ in C++. We run an instantiation for char-arithmetic and we print the gradient at $(101~77~114~114~32~121~109~88~115~97)^T$ to std::cout, yielding the output ``Merry Xmas''. Similar instantiations of type-generic second-order tangent and second-order adjoint versions of $$ y=\frac{1}{6} \cdot \sum_{i=0}^{n-1} x^3_{i} $$ yield ``Happy 2024'' at $(72~97~112~112~121~32~50~48~50~52)^T.$ Prepend a sufficiently large number of zeros to the input vector to explore the varying run times of the different derivative codes. The entire source code can be found on https://github.com/un110076/SeasonsGreetings., Comment: Code available on github

The increasing use of stochastic models for describing complex phenomena warrants surrogate models that capture the reference model characteristics at a fraction of the computational cost, foregoing potentially expensive Monte Carlo simulation. The predominant approach of fitting a large neural network and then pruning it to a reduced size has commonly neglected shortcomings. The produced surrogate models often will not capture the sensitivities and uncertainties inherent in the original model. In particular, (higher-order) derivative information of such surrogates could differ drastically. Given a large enough network, we expect this derivative information to match. However, the pruned model will almost certainly not share this behavior. In this paper, we propose to find surrogate models by using sensitivity information throughout the learning and pruning process. We build on work using Interval Adjoint Significance Analysis for pruning and combine it with the recent advancements in Sobolev Training to accurately model the original sensitivity information in the pruned neural network based surrogate model. We experimentally underpin the method on an example of pricing a multidimensional Basket option modelled through a stochastic differential equation with Brownian motion. The proposed method is, however, not limited to the domain of quantitative finance, which was chosen as a case study for intuitive interpretations of the sensitivities. It serves as a foundation for building further surrogate modelling techniques considering sensitivity information., Comment: 11 pages

We want to obtain derivatives in discontinuous program code, where default Algorithmic Differentiation may not perform well. Specifically, we consider discontinuities induced by control flow statements, where meaningful derivatives should ideally be capable of representing the resulting kinks in the trajectory. To achieve this, one can interpolate the trajectory at the control flow statements before taking the derivative. We formulate a method to efficiently interpolate between all boundaries induced by control flow in program code. Theoretically, code can be viewed as a series of piecewise continuous functions applied in succession. These functions are nested inside one another and result in a function composition with several cases. We interpret this function composition as a tree and devise a heuristic to identify paths that are relevant to the interpolation. This allows us to conceive a language that smoothly interpolates control-flow statements automatically and efficiently, making it fully differentiable., Comment: 15 pages, 4 figures

We are considering differential-algebraic equations with embedded optimization criteria (DAEOs) in which the embedded optimization problem is solved by global optimization. This actually leads to differential inclusions for cases in which there are multiple global optimizer at the same time. Jump events from one global optimum to another result in nonsmooth DAEs and thus reduction of the convergence order of the numerical integrator to first-order. Implementation of event detection and location as introduced in this work preserves the higher-order convergence behavior of the integrator. This allows to compute discrete tangents and adjoint sensitivities for optimal control problems., Comment: 16 pages, 4 figures

A modification of Newton's method for solving systems of $n$ nonlinear equations is presented. The new matrix-free method relies on a given decomposition of the invertible Jacobian of the residual into invertible sparse local Jacobians according to the chain rule of differentiation. It is motivated in the context of local Jacobians with bandwidth $2m+1$ for $m\ll n$. A reduction of the computational cost by $\mathcal{O}(\frac{n}{m})$ can be observed. Supporting run time measurements are presented for the tridiagonal case showing a reduction of the computational cost by $\mathcal{O}(n).$ Generalization yields the combinatorial Matrix-Free Newton Step problem. We prove NP-completeness and we present algorithmic components for building methods for the approximate solution. Inspired by adjoint Algorithmic Differentiation, the new method shares several challenges for the latter including the DAG Reversal problem. Further challenges are due to combinatorial problems in sparse linear algebra such as Bandwidth or Directed Elimination Ordering., Comment: 12 pages, 1 figure

All known elimination techniques for (first-order) algorithmic differentiation (AD) rely on Jacobians to be given for a set of relevant elemental functions. Realistically, elemental tangents and adjoints are given instead. They can be obtained by applying software tools for AD to the parts of a given modular numerical simulation. The novel generalized face elimination rule proposed in this article facilitates the rigorous exploitation of associativity of the chain rule of differentiation at arbitrary levels of granularity ranging from elemental scalar (state of the art) to multivariate vector functions with given elemental tangents and adjoints. The implied combinatorial Generalized Face Elimination problem asks for a face elimination sequence of minimal computational cost. Simple branch and bound and greedy heuristic methods are employed as a baseline for further research into more powerful algorithms motivated by promising first test results. The latter can be reproduced with the help of an open-source reference implementation.

We present and discuss seven different open problems in applied combinatorics. The application areas relevant to this compilation include quantum computing, algorithmic differentiation, topological data analysis, iterative methods, hypergraph cut algorithms, and power systems., Comment: 43 pages, 5 figures

Adjoint algorithmic differentiation by operator and function overloading is based on the interpretation of directed acyclic graphs resulting from evaluations of numerical simulation programs. The size of the computer system memory required to store the graph grows proportional to the number of floating-point operations executed by the underlying program. It quickly exceeds the available memory resources. Naive adjoint algorithmic differentiation often becomes infeasible except for relatively simple numerical simulations. Access to the data associated with the graph can be classified as sequential and random. The latter refers to memory access patterns defined by the adjacency relationship between vertices within the graph. Sequentially accessed data can be decomposed into blocks. The blocks can be streamed across the system memory hierarchy thus extending the amount of available memory, for example, to hard discs. Asynchronous i/o can help to mitigate the increased cost due to accesses to slower memory. Much larger problem instances can thus be solved without resorting to technically challenging user intervention such as checkpointing. Randomly accessed data should not have to be decomposed. Its block-wise streaming is likely to yield a substantial overhead in computational cost due to data accesses across blocks. Consequently, the size of the randomly accessed memory required by an adjoint should be kept minimal in order to eliminate the need for decomposition. We propose a combination of dedicated memory for adjoint $L$-values with the exploitation of remainder bandwidth as a possible solution. Test results indicate significant savings in random access memory size while preserving overall computational efficiency.

In Computational Science, Engineering and Finance (CSEF) scripts typically serve as the "glue" between potentially highly complex and computationally expensive external subprograms. Differentiability of the resulting programs turns out to be essential in the context of derivative-based methods for error analysis, uncertainty quantification, optimization or training of surrogates. We argue that it should be enforced by the scripting language itself through exclusive support of differentiable (smoothed) external subprograms and differentiable intrinsics combined with prohibition of nondifferentiable branches in the data flow. Illustration is provided by a prototype adjoint code compiler for a simple Python-like scripting language., Comment: 7 pages

Many modern numerical methods in computational science and engineering rely on derivatives of mathematical models for the phenomena under investigation. The computation of these derivatives often represents the bottleneck in terms of overall runtime performance. First and higher derivative tensors need to be evaluated efficiently. The chain rule of differentiation is the fundamental prerequisite for computing accurate derivatives of composite functions which perform a potentially very large number of elemental function evaluations. Data flow dependences amongst the elemental functions give rise to a combinatorial optimization problem. We formulate {\sc Chain Rule Differentiation} and we prove it to be NP-complete. Pointers to research on its approximate solution are given., Comment: 8 pages

We investigate errors in tangents and adjoints of implicit functions resulting from errors in the primal solution due to approximations computed by a numerical solver. Adjoints of systems of linear equations turn out to be unconditionally numerically stable. Tangents of systems of linear equations can become instable as well as both tangents and adjoints of systems of nonlinear equations, which extends to optima of convex unconstrained objectives. Sufficient conditions for numerical stability are derived., Comment: 7 pages

Second derivatives of mathematical models for real-world phenomena are fundamental ingredients of a wide range of numerical simulation methods including parameter sensitivity analysis, uncertainty quantification, nonlinear optimization and model calibration. The evaluation of such Hessians often dominates the overall computational effort. The combinatorial {\sc Hessian Accumulation} problem aiming to minimize the number of floating-point operations required for the computation of a Hessian turns out to be NP-complete. We propose a dynamic programming formulation for the solution of {\sc Hessian Accumulation} over a sub-search space. This approach yields improvements by factors of ten and higher over the state of the art based on second-order tangent and adjoint algorithmic differentiation.

Validation is a major challenge in differentiable programming. The state of the art is based on algorithmic differentiation. Consistency of first-order tangent and adjoint programs is defined by a well-known first-order differential invariant. This paper generalizes the approach through derivation of corresponding differential invariants of arbitrary order.

We propose a generalization of separability in the context of global optimization. Our results apply to objective functions implemented as differentiable computer programs. They are presented in the context of a simple branch and bound method. The often significant search space reduction can be expected to yield an acceleration of any global optimization method. We show how to utilize interval derivatives calculated by adjoint algorithmic differentiation to examine the monotonicity of the objective with respect to so called structural separators and how to verify the latter automatically., Comment: 15 pages, 1 figure

We are interested in the development of an algorithmic differentiation framework for computing approximations to tangent vectors to scalar and systems of hyperbolic partial differential equations. The main difficulty of such a numerical method is the presence of shock waves that are resolved by proposing a numerical discretization of the calculus introduced in Bressan and Marson [Rend. Sem. Mat. Univ. Padova, 94:79-94, 1995]. Numerical results are presented for the one-dimensional Burgers equation and the Euler equations. Using the essential routines of a state-of-the-art code for computational fluid dynamics (CFD) as a starting point, three modifications are required to apply the introduced calculus. First, the CFD code is modified to solve an additional equation for the shock location. Second, we customize the computation of the corresponding tangent to the shock location. Finally, the modified method is enhanced by algorithmic differentiation. Applying the introduced calculus to problems of the Burgers equation and the Euler equations, it is found that correct sensitivities can be computed, whereas the application of black-box algorithmic differentiation fails.

Many problems in science and engineering require uncertainty quantification that accounts for observed data. For example, in computational neuroscience, Neural Population Models (NPMs) are mechanistic models that describe brain physiology in a range of different states. Within computational neuroscience there is growing interest in the inverse problem of inferring NPM parameters from recordings such as the EEG (Electroencephalogram). Uncertainty quantification is essential in this application area in order to infer the mechanistic effect of interventions such as anaesthesia. This paper presents C++ software for Bayesian uncertainty quantification in the parameters of NPMs from approximately stationary data using Markov Chain Monte Carlo (MCMC). Modern MCMC methods require first order (and in some cases higher order) derivatives of the posterior density. The software presented offers two distinct methods of evaluating derivatives: finite differences and exact derivatives obtained through Algorithmic Differentiation (AD). For AD, two different implementations are used: the open source Stan Math Library and the commercially licenced dco/c++ tool distributed by NAG (Numerical Algorithms Group). The use of derivative information in MCMC sampling is demonstrated through a simple example, the noise-driven harmonic oscillator. And different methods for computing derivatives are compared. The software is written in a modular object-oriented way such that it can be extended to derivative based MCMC for other scientific domains.

The efficient computation of Jacobians represents a fundamental challenge in computational science and engineering. Large-scale modular numerical simulation programs can be regarded as sequences of evaluations of in our case differentiable modules with corresponding local Jacobians. The latter are typically not available. Tangent and adjoint versions of the individual modules are assumed to be given as results of algorithmic differentiation instead. The classical (Jacobian) matrix chain product formulation is extended with the optional evaluation of matrix-free Jacobian-matrix and matrix-Jacobian products as tangents and adjoints. We propose a dynamic programming algorithm for the minimization of the computational cost of such generalized Jacobian chain products without considering constraints on the available persistent system memory. In other words, the naive evaluation of an adjoint of the entire simulation program is assumed to be a feasible option. No checkpointing is required. Under the given assumptions we obtain optimal solutions which improve the best state of the art methods by factors of up to seven on a set of randomly generated problem instances of growing size., Comment: 14 pages

In this work we present useful techniques and possible enhancements when applying an Algorithmic Differentiation (AD) tool to the linear algebra library Eigen using our in-house AD by overloading (AD-O) tool dco/c++ as a case study. After outlining performance and feasibility issues when calculating derivatives for the official Eigen release, we propose Eigen-AD, which enables different optimization options for an AD-O tool by providing add-on modules for Eigen. The range of features includes a better handling of expression templates for general performance improvements, as well as implementations of symbolically derived expressions for calculating derivatives of certain core operations. The software design allows an AD-O tool to provide specializations to automatically include symbolic operations and thereby keep the look and feel of plain AD by overloading. As a showcase, dco/c++ is provided with such a module and its significant performance improvements are validated by benchmarks., Comment: Updated with accepted version for ICCS 2020 conference proceedings. The final authenticated publication is available online at https://doi.org/10.1007/978-3-030-50371-0_51. See v1 for the original, extended preprint. 14 pages, 7 figures

A new proof for adjoint systems of linear equations is presented. The argument is built on the principles of Algorithmic Differentiation. Application to scalar multiplication sets the base line. Generalization yields adjoint inner vector, matrix-vector, and matrix-matrix products leading to an alternative proof for first- as well as higher-order adjoint linear systems.

With recent developments in cloud computing, a paradigm shift from rather static deployment of resources to more dynamic, on-demand practices means more flexibility and better utilization of resources. This demands new ways to efficiently configure networks. In this paper, we will characterize a class of competitive cloud services that telecom operators could provide based on the characteristics of telecom infrastructure through an applicable streaming service architecture. Then, we will model this architecture as a cost-based mathematic model. This model provides a tool to evaluate and compare the cost of software services for different telecom network topologies and deployment strategies. Additionally, with each topology it acts as a means to characterize the deployment solution that yields the lowest resource usage over the entire network. These applications are illustrated through numerical analysis. Finally, a proof-of-concept prototype is deployed to shows dynamic properties of the service in the architecture and the model above., Comment: 20 pages; International Journal on Cloud Computing: Services and Architecture (IJCCSA), 2011

The South Pole Acoustic Test Setup (SPATS) was built to evaluate the acoustic characteristics of the South Pole ice in the 10 kHz to 100 kHz frequency range, for the purpose of assessing the feasibility of an acoustic neutrino detection array at the South Pole. The SPATS hardware consists of four vertical strings deployed in the upper 500 m of the South Pole ice cap. The strings form a trapezoidal array with a maximum baseline of 543 m. Each string has 7 stages equipped with one transmitter and one sensor module. Sound is detected or generated by piezoelectric ceramic elements inside the modules. Analogue signals are sent to the surface on electric cables where they are digitized by a PC-based data acquisition system. The data from all strings are collected on a central computer in the IceCube Laboratory from where they are send to a central data storage facility via a satellite link or stored locally on tape. A technical overview of SPATS and its performance is presented.

We give a gentle introduction to using various software tools for automatic differentiation (AD). Ready-to-use examples are discussed, and links to further information are presented. Our target audience includes all those who are looking for a straightforward way to get started using the available AD technology. The document is dynamic in the sense that its content will be updated as the AD software evolves., Comment: 23 pages