HUMBOLDT–UNIVERSITÄT ZU BERLIN

Mathematisch-Naturwissenschaftliche Fakultät II Institut für Mathematik Prof. PhD. Andreas Griewank Kshitij Kulshreshtha

Humboldt-Universität zu Berlin, Institut für Mathematik, Unter den Linden 6, D-10099 Berlin

Introduction to Algorithmic Differentiation

Wintersemester 2004/05

Practice Excercise 2

To be finished by November 10, 2004

Excercise 2.1 Consider an evaluation procedure that contains no additions and subtractions, so that none of the elemental partials c_{ij} = ∂ϕ_{i}/∂v_{j} with j ≺ i has one of the special values in {−1, 0, +1}. Moreover, assume that there are no dead ends in that all intermediates v_{k} with 1 ≤ k ≤ l − m depend on at least one independent variable x_{j} and impact at least one dependent variable y_{i}. Finally, assume that no dependent is directly calculated from an independent, so that none of the c_{ij} itself is an element of the Jacobian. Then each value c_{ij} affects the value of the Jacobian element ∂y_{i}/∂x_{j} and must therefore enter somehow into any procedure for evaluating F′(x) at least once. Let the gradient ∇f(x) = \overline{F} (x, 1) of a scalarvalued function f = F with m = 1 be calculated by the adjoint procedure with \overline{y}_{1} = 1. Then each elemental partial c_{ij} enters exactly once, namely, as a factor in a multiplication. In the incremental form each one of these multiplications is followed by an addition.

a. Conclude that the number of arithmetic operations used to accumulate ∇f from the elemental partials in the reverse mode cannot exceed the minimal number needed to calculate ∇f by a factor of more than 4. For some d ≫ 1 and positive constants w_{i}, z_{i}, b_{i} for i = 1 … d, consider the function

f(x) = ln |∑_{i=1}^{d} w_{i}(exp(z_{i} ∗ sin(x)) − b_{i})^{2}| .

We might think of sin(x) as a parameter restricted to [−1, 1] that should be chosen to minimise the discrepancies between the values exp(z_{i} ∗ sin(x)) and the b_{i} as measured by the logarithm of a suitably weighted sum. To optimise such an exponential fit, one obviously needs the scalar-valued derivative f′(x). To satisfy the special assumptions made above, let us consider the ψ_{i}(u) ≡ (u − b_{i})^{2} and the weighted sum ψ_{0}(v) = ∑_{i=1}^{d} w_{i}v_{i} as elemental functions.

b. Write a procedure for evaluating f using elementaries ψ_{i} for i = 0, 1, 2, …, d and draw the computational graph.

c. Derive the corresponding adjoint procedure for accumulating ∇f from the elemental partials. Count the number of multiplications and additions used in the incremental form. Verify that the number of multiplications corresponds exactly to the number of edges, which is 4d + O(1).

d. Rewrite the formula for f(x) using ψ_{i}, and differentiate it by hand. Show that by keeping the common factors cos(x) and 1/ψ_{0} out of the internal sum, the number of multiplcations needed for the calculation can be kept to 3d + O(1). e. Show that the tangent procedure for calculating f′(x) = ḟ(x, ẋ) with ẋ = 1 involves

exactly as many multiplications as the reverse mode and is therefore also less economical than the procedure derived by hand in part d.

Excercise 2.2 Consider the function f(X) ≡ ln | det(X)|, where X ∈ ℝ^{n̂×n̂} is a square matrix whose n ≡ n̂×n̂ entries X_{ij} are considered as independent variables. This function f(X) has been used as a penalty term for maintaining positive definiteness of X when it is restricted to being symmetric.

a. Show that, provided det(X) ≠ 0 so that f(X) can be evaluated, the matrix G ∈ ℝ^{n̂×n̂} of corresponding gradient components G_{ij} = G_{ij}(X) ≡ ∂f(X)/∂X_{ij} is the transpose of the inverse X^{−1}. Hence the cheap gradient principle guearentees that any finite algorith for evaluating the determinant of a general square matrix can be transformed into an algorithm for computing the inverse with essentially the same temporal complexity.

b. Examine whether the last assertion is still true when X is restricted to having a certain sparsity pattern (for example, tridiagonal). Explain your conclusion.

c. Assume that f(X) can be evaluated on some neighbourhood of a particular matrix by performing LU factorisation without pivoting. Write this algorithm asd an evaluation procedure that overwrites the original matrix elements with intermediate quantities and finally with the non-trivial entries of the two triangular factors. Derive the corresponding adjoint procedure with tape and compare it to standard methods for computing matrix inverses.

d. For a constant matrix A ∈ ℝ^{n×n} and a vector b ∈ ℝ^{n}, determine the gradient of

f(x) ≡ ln | det(A + bx^{T})| at x = 0 .

e. Write an evaluation procedure that evaluates f(x) for n = 2 by first computing the entries A + bx^{T}, then applying the explicit determinant formula for 2 × 2 matricesm and finally taking the logarithm of the determinant’s modulus. Show that the resulting adjoint procedure for evaluating ∇f is not optimal.

f. Suppose that A + bx^{T} has for all x near the origin n simple non-zero real eigenvalues λ_{j}(x) so that f(x) = ∑_{j} ln |λ_{j}(x)|. Determine ∇_{x}λ_{j} at x = 0 by implicit differentiation of the identities (A+bx^{T})v_{j} = λ_{j}v_{j} and (A+bx^{T})^{T} w_{j} = λ_{j}w_{j}. Here v_{j}, w_{j} ∈ ℝ^{n} denote the corresponding left and right eigenvectors normalised such that v_{j}^{T} w_{j} = 1 for j = 1, 2, …, n. Write an evaluation procedure for f(x) treating the λ_{j} as elemental functions. Derive and interpret the corresponding adjoint procedure for evaluating ∇f(0).

Phone: (030) 2093-5414 Faxphone: (030) 2093-5848

e-mail: griewank@math.hu-berlin.de, kshitij@math.hu-berlin.de
WWW: http://www.mathematik.hu-berlin.de/∼gaggle/AD/ADuebung.html

2