Math. Program., Ser. A 95: 555–571 (2003) Digital Object Identifier (DOI) 10.1007/s10107-002-0329-7 Andreas Griewank · Uwe Naumann Accumulating Jacobians as chained sparse matrix products Received: January 17, 2002 / Accepted: May 29, 2002 Published online: February 14, 2003 – c  Springer-Verlag 2003 Abstract. The chain rule – fundamental to any kind of analytical differentiation - can be applied in vari- ous ways to computational graphs representing vector functions. These variants result in different operations counts for the calculation of the corresponding Jacobian matrices. The minimization of the number of arithme- tic operations required for the calculation of the complete Jacobian leads to a hard combinatorial optimization problem. We will describe an approach to the solution of this problem that builds on the idea of optimizing chained matrix products using dynamic programming techniques. Reductions by a factor of 3 and more are possible regarding the operations count for the Jacobian accumulation. After discussing the mathematical basics of Automatic Differentiation we will show how to compute Jacobians by chained sparse matrix products. These matrix chains can be reordered, must be pruned, and are finally subject to a dynamic programming algorithm to reduce the number of scalar multiplications performed. Key words. chained matrix product – combinatorial optimization – dynamic programming – edge elimination in computational graphs 1. Introduction – derivatives of programs Many numerical computations require the evaluation of Jacobians F ′ (x) ∈ R n (1) for vector functions F : R n → R m . Automatic Differentiation (AD) provides a fast and convenient way for calculating de- rivatives with machine accuracy by simply exploiting the chain rule. This approach is completely different from the one of approximating derivatives numerically by differ- ence quotients. In most cases F is given as a computer program written in some high-level program- ming language such as C or Fortran which may be decomposed into a sequence of scalar assignments of the form (R ∋) v i = ϕ i (v j ≺i ) for i = 1,... ,l + m (2) Here j ≺ i means that v j is an argument of the elemental function ϕ i that computes v i . Normally we expect that only one or two scalar arguments occur, but this requirement A. Griewank: Technical University Dresden, Germany e-mail: griewank@math.tu-dresden.de U. Naumann: Argonne National Laboratory, MCS, Il e-mail: naumann@mcs.anl.gov