Tutorials

Clad is an open source clang plugin which supports automatic differentiation of mathematical functions in C++. This page walks through one small example per mode.

Every example below is a file in Clad’s test suite, included here verbatim. The // prints: comments are what the suite checks the program writes, so the code and the numbers beside it are kept true by the build.

The Forward mode

Clad supports forward mode automatic differentiation through the clad::differentiate API call. It differentiates with respect to one parameter, named in the second argument, and the generated function returns the derivative.

#include "clad/Differentiator/Differentiator.h"
#include <iostream>

double func(double x) { return x * x; }

int main() {
  // Ask clad for the derivative of func with respect to x.
  auto d_func = clad::differentiate(func, "x");
  // Call it the way func itself would be called.
  std::cout << d_func.execute(/*x=*/3) << "\n"; // prints: 6
  // And print the code clad generated for it.
  d_func.dump();
}

.dump() prints the derivative Clad generated, which is often the quickest way to see what it did.

The Reverse Mode

Clad also supports reverse mode automatic differentiation, through the clad::gradient API call. One call computes the derivatives with respect to every parameter named in args, and writes each one through a pointer the caller supplies.

#include "clad/Differentiator/Differentiator.h"
#include <iostream>

double f(double x, double y, double z) { return x * y * z; }

int main() {
  auto d_f = clad::gradient(f, "x, y");
  // One adjoint per differentiated parameter, and clad accumulates into them,
  // so they start at zero.
  double dx = 0, dy = 0;
  d_f.execute(/*x=*/2, /*y=*/3, /*z=*/4, &dx, &dy);
  std::cout << "dx: " << dx << ", dy: " << dy << "\n"; // prints: dx: 12, dy: 8
}

The example differentiates with respect to x and y. Naming one parameter, as in clad::gradient(f, "x"), differentiates with respect to that one; naming none, as in clad::gradient(f), differentiates with respect to all of them.

The Hessian Mode

Clad can also produce a hessian matrix through the clad::hessian API call. It returns the matrix as a flattened array in row major order, so \(n\) independent variables need \(n^2\) elements.

#include "clad/Differentiator/Differentiator.h"
#include <iostream>

double f(double x, double y, double z) { return x * y * z; }

int main() {
  // Two independent variables, so the hessian is 2 x 2 and needs 4 elements.
  auto f_hess = clad::hessian(f, "x, y");
  double matrix_f[4] = {0};
  f_hess.execute(3, 4, 5, matrix_f);
  std::cout << "[" << matrix_f[0] << ", " << matrix_f[1] << "\n"
            << " " << matrix_f[2] << ", " << matrix_f[3] << "]\n";
  // prints: [0, 5
  // prints:  5, 0]
}

When an array is involved, say which elements to differentiate with respect to: for double f_arr(double x, double y, double z[2]), the call clad::hessian(f_arr, "x, y, z[0:1]") uses four independent variables and so needs sixteen elements.

The Jacobian Mode

Clad can produce the jacobian of a function using its vectorized forward mode. It returns the jacobian as a clad::matrix for every pointer or array parameter.

#include "clad/Differentiator/Differentiator.h"
#include <iostream>

void f(double x, double y, double z, double* output) {
  output[0] = x * y;
  output[1] = y * y * x;
  output[2] = 6 * x * y * z;
}

int main() {
  auto f_jac = clad::jacobian(f);

  // One row per output element; one column per independent scalar, counting
  // the three elements output itself contributes.
  clad::matrix<double> d_output(3, 6);
  double output[3] = {0};
  f_jac.execute(3, 4, 5, output, &d_output);
  for (int row = 0; row < 3; ++row)
    std::cout << d_output[row][0] << " " << d_output[row][1] << " "
              << d_output[row][2] << "\n";
  // prints: 4 3 0
  // prints: 16 24 0
  // prints: 120 90 72
}

The matrix has one row per element of the output and one column per independent scalar. Here the output has three elements, and the independent scalars are x, y, z and the three elements of output itself, which is why it is \(3 \times 6\). The last three columns come out zero here, because output does not depend on its own previous contents.

Error Estimation API

Clad is capable of annotating a given function with floating point error estimation code using reverse mode AD.

#include "clad/Differentiator/Differentiator.h"
#include <iomanip>
#include <iostream>

double func(double x, double y) { return x * y; }

int main() {
  auto dfunc_error = clad::estimate_error(func);

  // The gradient's arguments, plus a double& that receives the error estimate.
  double x = 3, y = 5, d_x = 0, d_y = 0, final_error = 0;
  dfunc_error.execute(x, y, &d_x, &d_y, final_error);

  std::cout << std::setprecision(3) << final_error << "\n"; // prints: 5.36e-06
}

The signature is the one clad::gradient would generate, with one extra argument of type double& at the end, which receives the total floating point error.