# Verification¶

## Taylor remainder convergence test¶

The fundamental tool used in verification of gradients is the
*Taylor remainder convergence test*. Let \(\widehat{J}(m)\) be the functional, considered
as a pure function of the parameter of interest,
let \(\nabla \widehat{J}\) be its gradient, and let \(\delta m\) be a perturbation to
\(m\). This test is based on the observation that

but that

by Taylor’s theorem. The quantity \(\left|\widehat{J}(m + h\delta m) - \widehat{J}(m)\right|\) is called the *first-order
Taylor remainder* (because it’s supposed to be first-order), and the quantity \(\left|\widehat{J}(m + h\delta m) - \widehat{J}(m) - h\nabla \widehat{J} \cdot \delta m\right|\)
is called the *second-order Taylor remainder*.

Suppose someone gives you two functions \(\widehat{J}\) and \(d\widehat{J}\), and claims that \(d\widehat{J}\) is the gradient of \(\widehat{J}\). Then their claim can be rigorously verified by computing the second-order Taylor remainder for some choice of \(h\) and \(\delta m\), then repeatedly halving \(h\) and checking that the result decreases by a factor of 4.

## Applying this in dolfin-adjoint¶

In the case of PDE-constrained optimisation, computing \(\widehat{J}(m)\) involves solving the PDE
for that choice of \(m\) to compute the solution \(u\), and then evaluating the functional \(J\).
The main function in dolfin-adjoint for applying the Taylor remainder convergence test is `taylor_test`

.
To see how this works, let’s restructure our forward model so that we can run it as a pure function of
the diffusivity \(\nu\):

```
from dolfin import *
from dolfin_adjoint import *
n = 30
mesh = UnitSquareMesh(n, n)
V = VectorFunctionSpace(mesh, "CG", 2)
ic = project(Expression(("sin(2*pi*x[0])", "cos(2*pi*x[1])"), degree=2), V)
def main(nu):
u = ic.copy(deepcopy=True)
u_next = Function(V)
v = TestFunction(V)
timestep = Constant(0.01)
F = (inner((u_next - u)/timestep, v)
+ inner(grad(u_next)*u_next, v)
+ nu*inner(grad(u_next), grad(v)))*dx
bc = DirichletBC(V, (0.0, 0.0), "on_boundary")
t = 0.0
end = 0.1
while (t <= end):
solve(F == 0, u_next, bc)
u.assign(u_next)
t += float(timestep)
return u
if __name__ == "__main__":
nu = Constant(0.0001)
u = main(nu)
J = Functional(inner(u, u)*dx*dt[FINISH_TIME])
dJdnu = compute_gradient(J, Control(nu))
```

As you can see, we’ve taken the action part of the model into a function `main`

, so that we
can drive the forward model several times for verification. Now let’s see how to use `taylor_test`

:

```
from dolfin import *
from dolfin_adjoint import *
n = 30
mesh = UnitSquareMesh(n, n)
V = VectorFunctionSpace(mesh, "CG", 2)
ic = project(Expression(("sin(2*pi*x[0])", "cos(2*pi*x[1])"), degree=2), V)
def main(nu):
u = ic.copy(deepcopy=True)
u_next = Function(V)
v = TestFunction(V)
timestep = Constant(0.01)
F = (inner((u_next - u)/timestep, v)
+ inner(grad(u_next)*u_next, v)
+ nu*inner(grad(u_next), grad(v)))*dx
bc = DirichletBC(V, (0.0, 0.0), "on_boundary")
t = 0.0
end = 0.1
while (t <= end):
solve(F == 0, u_next, bc)
u.assign(u_next)
t += float(timestep)
return u
if __name__ == "__main__":
nu = Constant(0.0001)
u = main(nu)
J = Functional(inner(u, u)*dx*dt[FINISH_TIME])
dJdnu = compute_gradient(J, Control(nu))
Jnu = assemble(inner(u, u)*dx) # current value
parameters["adjoint"]["stop_annotating"] = True # stop registering equations
def Jhat(nu): # the functional as a pure function of nu
u = main(nu)
return assemble(inner(u, u)*dx)
conv_rate = taylor_test(Jhat, Control(nu), Jnu, dJdnu)
```

Download the adjoint code with verification.

We set `parameters["adjoint"]["stop_annotating"]`

to `True`

, so that the `solve`

calls
in the perturbed runs of `main`

are not registered by libadjoint.

Running this program yields the following output:

```
$ python tutorial4.py
...
Taylor remainder without gradient information:
[0.0023634768826859553, 0.001219686435181555, 0.0006197555788530762,
0.0003124116082189321, 0.0001568463925042951]
```

The first line gives the values computed for the first-order Taylor remainder. As you can see, each value is approximately half the previous one.

```
Convergence orders for Taylor remainder without gradient information (should all be 1):
[0.9544004555219242, 0.9767390399645689, 0.9882512926546192, 0.9940957131087177]
```

The second line shows the convergence orders of the first-order Taylor remainders: these should
always be 1. (If they are not, try passing the `seed`

argument to `taylor_test`

to
use a smaller perturbation.)

```
Taylor remainder with gradient information:
[0.00015639195264909554, 4.0247982485970384e-05, 1.0211629980686528e-05,
2.5719961979492594e-06, 6.454097041455739e-07]
```

The third line gives the values computed for the second-order Taylor remainder. These values should be much smaller than those on the first line.

```
Convergence orders for Taylor remainder with gradient information (should all be 2):
[1.9581779067535698, 1.9787032993204938, 1.9892527525050359, 1.9946013350664813]
```

The fourth line shows the convergence orders of the second-order Taylor remainders: if the gradient has been computed correctly with the adjoint, then these numbers should be 2.

As can be seen, the second-order Taylor remainders do indeed converge at second order, and so the gradient `dJdnu`

is correct.

So, what if the Taylor remainders are not correct? Such a situation could occur if the model
manually modifies `Function`

values, or if the model modifies the entries of assembled matrices and
vectors, or if the model is not differentiable, or if there is a bug in dolfin-adjoint. dolfin-adjoint offers many ways to pinpoint
precisely where an error might lie; these are discussed in the section on debugging.

Once the adjoint model development is completed, you may wish to run your model on much bigger
production runs. **By default, dolfin-adjoint stores all variables computed in memory**, as they
will be necessary to linearise the model about the forward trajectory. If you wish to solve much
bigger problems, or if the model must run for many timesteps, it may not be feasible to store
everything: some balance between storage and recomputation must be achieved, in the form of a
**checkpointing scheme**. Checkpointing is the topic of the next section.