Optimized Sparse Matrix Operations for Reverse Mode Automatic Differentiation

The source transformation tool for automatic differentiation of Fortran programs ADIFOR uses a preaccumulation technique to speed up tangent-linear codes significantly compared to the standard forward mode. Reverse mode automatic differentiation is applied to all scalar assignments to generate efficient code for the computation of local gradients. It has been well known for some time that reverse mode is not necessarily the optimal choice for the computation of these statement-level gradients as it does not minimize the number of operations required. This article presents an efficient algorithm for the solution of this combinatorial optimization problem. The corresponding software is freely available for downloading on our website. Developers of software for automatic differentiation are invited to integrate the algorithm into their tools. Gradients of scalar multivariate functions can be computed by elimination methods on the linearized computational graph. The combinatorial optimization problem that aims to minimize the number of arithmetic operations performed by the elimination algorithm is known to be NP-complete. In this article we present a polynomial algorithm for solving a relevant subclass of this problem's instances. The proposed method relies on the ability to compute vertex covers in bipartite graphs in polynomial time. A simplified version of this graph algorithm is used in a research prototype of the differentiation-enabled NAGWare Fortran compiler for the preaccumulation of local gradients of scalar assignments in the context of automatic generation of efficient tangent-linear code for numerical programs.

Automatic differentiation (AD) has become an important tool for optimization problems in computational science, and it has been applied to the Hartree-Fock method. Although the reverse-mode AD is more efficient than the forward-mode, eigenvalue calculation in the self-consistent field (SCF) method has impeded the use of the reverse-mode AD. Here, we propose a method to directly minimize Hartree-Fock energy under the orthonormality constraint of the molecular orbitals using reverse-mode AD by avoiding eigenvalue calculation. According to our validation, the proposed method was more stable than the conventional SCF method and achieved comparable accuracy.

Automatic differentiation (AD) is conventionally understood as a family of distinct algorithms, rooted in two “modes”—forward and reverse—which are typically presented (and implemented) separately. Can there be only one? Following up on the AD systems developed in the JAX and Dex projects, we formalize a decomposition of reverse-mode AD into (i) forward-mode AD followed by (ii) unzipping the linear and non-linear parts and then (iii) transposition of the linear part. To that end, we define a (substructurally) linear type system that can prove a class of functions are (algebraically) linear. Our main results are that forward-mode AD produces such linear functions, and that we can unzip and transpose any such linear function, conserving cost, size, and linearity. Composing these three transformations recovers reverse-mode AD. This decomposition also sheds light on checkpointing, which emerges naturally from a free choice in unzipping let expressions. As a corollary, checkpointing techniques are applicable to general-purpose partial evaluation, not just AD. We hope that our formalization will lead to a deeper understanding of automatic differentiation and that it will simplify implementations, by separating the concerns of differentiation proper from the concerns of gaining efficiency (namely, separating the derivative computation from the act of running it backward).

This document summarizes the findings of a strategic planning exercise commissioned by the Weapons Simulation and Computing, Computational Physics (WSC/CP) program at the Lawrence Livermore National Laboratory (LLNL) in FY24. During the year, the committee met with multiple stakeholder communities to gather input, opinions, suggestions and concerns which have been incorporated throughout this document. The key findings from this exercise are summarized: • The development of multiphysics modelling and simulation (mod/sim) codes and software technologies, their deployment on exascale compute platforms, and their broad adoption across the NNSA is a major success of the Advanced Simulation and Computing (ASC) program and the Exascale Computing Project (ECP). Sustained investment in these core technologies is essential. • Today’s state of the art involves running ensembles of O(100K) simulations to perform uncertainty quantification (UQ) and design studies using multiple statistical methods such as Bayesian optimization to understand sensitivities of our models and explore parameterized design spaces. Even with exascale computing, we are practically limited to O(10) parameters in these studies since the number of simulations required to sample the space scales exponentially with the number of design parameters. • The data from these simulation ensembles is increasingly being used to train machine learned (ML) surrogates (or reduced order models, ROMs) which can then be used for optimization or real time design exploration. However, the trained surrogates are still limited in the number of parameters they can represent due to the sampling limitations previously noted. • Augmenting our suite of integrated multiphysics simulation codes, both current and emerging, with the ability to compute gradients (solution derivatives) of arbitrary simulation outputs with respect to (some or all) simulation inputs would be a breakthrough technology, opening the door to a new era of efficient and automated inverse design based on verified and validated mod/sim capabilities. • This capability, which we refer to as differentiable multiphysics codes (DMCs), would revolutionize both UQ and optimization studies by breaking the curse of dimensionality that presently limits our “gradient-free” ensemble based computing approach. A similar breakthrough occurred in the AI/ML community once the ability to compute gradients of arbitrary loss functions using back-propagation became commonplace. Gradient information from the multiphysics codes can also be used to dramatically improve the efficiency and scale of training of ML/ROM surrogates for rapid assessments. • Achieving this in our suite of codes will be a grand challenge, similar to the amount of effort that was required to transition from CPU to GPU computing. It will require buy-in from the entire WSC/CP program and beyond, including all integrated codes, physics and engineering models, third-party library dependencies and performance portability abstractions. It will also require investment in research and development of numerical methods for computing adjoints of coupled physics across multiple adaptively refined moving meshes and of stochastic (Monte Carlo) and mesh free (SPH) methods. • New software and numerical techniques, largely pioneered by the AI/ML community, make this feasible. Chief among these is automatic differentiation (AD), the ability to employ AD at point-wise locations in a physics calculation (instead of traditional black-box approaches) and the ability to perform “back-propagation in time” (or reverse mode AD) for non-linear partial differential equations (PDEs). Fundamentally, the conclusion of this strategic planning exercise is that the time is right to undertake a large scale effort in WSC, centered on the existing integrated codes, to continue the natural evolution of mod/sim in the age of AI/ML. Instead of attempting to replace mod/sim with purely data driven AI/ML models, we believe the key to success is to integrate AI/ML by building on top of the decades of hard-won knowledge and the verified/validated multiphysics modelling capability that is the hallmark of the ASC program.

For dynamic multileaf collimator-based intensity modulatedradiotherapy in which small beam elements are used to generatecontinuous modulation, the sheer size of the dose calculationmatrix could pose serious computational challenges. In order tocircumvent this problem, the dose calculation matrix was reducedto a sparse matrix by truncating the weakly contributing entriesbelow a certain cutoff to zero. Subsequently, the sparse matrix wascompressed and matrix indexing vectors were generated to facilitatematrix-vector and matrix-matrix operations used in inverse planning.The application of sparsity permitted the reduction of overallmemory requirement by an order of magnitude.In addition, the effect of disregarding the small scattercomponents on the quality of optimization was investigated byrepeating the inverse planning using the dense dose calculationmatrix. Comparison of dense and sparse matrix-based plansrevealed an insignificant difference in optimization outcome, thusdemonstrating the feasibility and usefulness of the sparse methodin inverse planning. Furthermore, two additional methods of memoryminimization are suggested, namely hexagonal dose sampling andlimited normal tissue sampling.

This article presents our work toward correct and efficient automatic differentiation of OpenMP parallel worksharing loops in forward and reverse mode. Automatic differentiation is a method to obtain gradients of numerical programs, which are crucial in optimization, uncertainty quantification, and machine learning. The computational cost to compute gradients is a common bottleneck in practice. For applications that are parallelized for multicore CPUs or GPUs using OpenMP, one also wishes to compute the gradients in parallel. We propose a framework to reason about the correctness of the generated derivative code, from which we justify our OpenMP extension to the differentiation model. We implement this model in the automatic differentiation tool Tapenade and present test cases that are differentiated following our extended differentiation procedure. Performance of the generated derivative programs in forward and reverse mode is better than sequential, although our reverse mode often scales worse than the input programs.

The matrix computation language and environment MATLAB is extended to include sparse matrix storage and operations. The only change to the outward appearance of the MATLAB language is a pair of commands to create full or sparse matrices. Nearly all the operations of MATLAB now apply equally to full or sparse matrices, without any explicit action by the user. The sparse data structure represents a matrix in space proportional to the number of nonzero entries, and most of the operations compute sparse results in time proportional to the number of arithmetic operations on nonzeros.

The linear and nonlinear harmonic methods are efficient frequency domain methods for analyzing time periodic unsteady flow fields. They have been widely used in both academia and industry. But the cost and complexity of developing a linear harmonic solver has been limiting its wider applications. On the other hand, the automatic differentiation (AD) has long been used in the CFD community with a focus on generating adjoint codes in a reverse mode. All those AD tools can do a much better job in generating linearized codes in a tangent mode, but so far very little, if any, attention is paid to using AD for developing linear harmonic solvers. The linear harmonic method, in comparison with the harmonic balance method, has its own advantages. For example, it can capture small disturbances very effectively, and avoids aliasing errors which can lead to solution instability since each wave component is solved for separately. This paper presents the effort of using an AD tool to generate major source codes for the development of a linear harmonic solver for analyzing time periodic unsteady flows. It includes the procedures and advice of using AD for such a purpose. A case study is also presented to validate the developed linear harmonic solver.

Memory access is known to be the main bottleneck for shared-memory parallel graph applications especially for large and irregular graphs. Propagation blocking (PB) idea was proposed recently to improve the parallel performance of PageRank and sparse matrix and vector multiplication operations. The idea is based on separating parallel computation into two phases, binning and accumulation, such that random memory accesses are replaced with contiguous accesses. In this paper, we propose an algorithm that allows execution of these two phases concurrently. We propose several improvements to increase parallel throughput, reduce memory overhead, and improve work efficiency. Our experimental results show that our proposed algorithms improve shared-memory parallel throughput by a factor of up to 2× compared to the original PB algorithms. We also show that the memory overhead can be reduced significantly (from 170 percent down to less than 5 percent) without significant degradation of performance. Finally, we demonstrate that our concurrent execution model allows asynchronous parallel execution, leading to significant work efficiency in addition to throughput improvements.

SimRank is an important metric to measure the topological similarity between two nodes in a graph. In particular, single-source and top-k SimRank has numerous applications in recommendation systems, network analysis, and web mining, etc. Mathematically, given a vertex, the computation of single-machine and single-source SimRank mainly lies in matrix-matrix operations. However, it is almost impossible to directly compute on large graphs. Thus, existing works yield to two main operations: a series of random walks, and sparse matrix and dense vector multiplication operations. This brings about high computation cost for SimRank on large graphs. In real-world applications, there is always the query time and accuracy trade-off, which hinders the computation of high-precision SimRank on large-scale graphs. To handle this problem, this paper proposesClipSim, the first GPU-friendly parallel framework that accelerates the single-source SimRank on GPU with accuracy guarantee. We design a novel data structure and GPU-friendly parallel algorithms for efficient computation of all the operations of SimRank on GPU. Moreover, our theoretical derivation enables ClipSim to largely reduce the number of random walks required for each node, while maintaining the same theoretical accuracy as the state-of-the-art algorithm, ExactSim. We conduct extensive experiments on real-world and synthetic datasets to demonstrate the accuracy and efficiency of ClipSim. The results show that compared with ExactSim, ClipSim obtains single-source SimRank vectors with the same accuracy and up to 160× faster computation time.

This paper provides a brief introduction to automatic differentiation and relates it to the tangent linear model and adjoint approaches commonly used in meteorology. After a brief review of the forward and reverse mode of automatic differentiation, the ADIFOR automatic differentiation tool is introduced, and initial results of a sensitivity-enhanced version of the MM5 PSU/NCAR mesoscale weather model are presented. We also present, a novel approach to the computation of gradients that, uses a reverse mode approach at the time loop level and a forward mode approach at every time step. The resulting “pseudoadjoint” shares the characteristic of an adjoint code that the ratio of gradient to function evaluation does not depend on the number of independent variables. In contrast to a true adjoint approach, however, the nonlinearity of the model plays no role in the complexity of the derivative code.Key wordsAutomatic DifferentiationAdjointTangent Linear ModelMM5ADIFORSparsLinCData Assimilation

We study the correctness of automatic differentiation (AD) in the context of a higher-order, Turing-complete language (PCF with real numbers), both in forward and reverse mode. Our main result is that, under mild hypotheses on the primitive functions included in the language, AD is almost everywhere correct, that is, it computes the derivative or gradient of the program under consideration except for a set of Lebesgue measure zero. Stated otherwise, there are inputs on which AD is incorrect, but the probability of randomly choosing one such input is zero. Our result is in fact more precise, in that the set of failure points admits a more explicit description: for example, in case the primitive functions are just constants, addition and multiplication, the set of points where AD fails is contained in a countable union of zero sets of polynomials.

In mathematics and computer algebra, automatic differentiation (AD) is a set of techniques to evaluate the derivative of a function specified by a computer program. AD exploits the fact that every computer program, no matter how complicated, executes a sequence of elementary arithmetic operations (addition, subtraction, multiplication, division, etc.), elementary functions (exp, log, sin, cos, etc.) and control flow statements. AD takes source code of a function as input and produces source code of the derived function. By applying the chain rule repeatedly to these operations, derivatives of arbitrary order can be computed automatically, accurately to working precision, and using at most a small constant factor more arithmetic operations than the original program. This paper presents AD techniques available in ROOT, supported by Cling, to produce derivatives of arbitrary C/C++ functions through implementing source code transformation and employing the chain rule of differential calculus in both forward mode and reverse mode. We explain its current integration for gradient computation in TFormula. We demonstrate the correctness and performance improvements in ROOT’s fitting algorithms.

Several software systems are available for implementing automatic differentiation of computer programs. The forward mode of automatic differentiation is limited by computational intensity and computer memory. The reverse mode, or adjoint approach, is limited by computer memory and disk storage. A modular technique for derivative computation that can significantly reduce memory required to compute derivatives in a complex FORTRAN model using the reverse mode of automatic differentiation is discussed and demonstrated.

Automatic differentiation is applied for improving productivity when solving optimal control problems. A gradient operator is evaluated with an algorithm which uses automatic differentiation in reverse mode. This technique is combined with a dedicated compiler with the use of which a symbolic rather than a programming language can be applied for expanding the formal description of optimal control problems. Numerical results concerning the optimal ascent of an orbital stage are presented.