Using Clad

This section briefly describes all the key functionalities offered by Clad. If you are just getting started with Clad, then this is the best place to start. You may want to skim some sections on the first read.

In case you haven’t installed Clad already, then please do before proceeding with this guide. Visit Installation and usage to know more about installing clad.

Let’s get started.

Automatic Differentiation

Clad differentiation functions takes a function as an input and returns a clad::CladFunction object that contains information about the generated derived function. Generated derived function can be called by calling the .execute method on the corresponding clad::CladFunction object.

Clad consists of five primary automatic differentiation functions:

  • clad::differentiate – Primary forward mode automatic differentiation

  • clad::gradient – Primary reverse mode automatic differentiation

  • clad::hessian

  • clad::jacobian

  • clad::estimate_error

Each of these functions will be explored in this guide.

Which one you want depends on two things you already know about your function: how many inputs and outputs it has, and whether you need first or second derivatives.

        flowchart TD
  Q1{"Do you need<br/>second derivatives?"}
  Q2{"How many inputs<br/>and outputs?"}
  H["clad::hessian"]
  D["clad::differentiate"]
  G["clad::gradient"]
  J["clad::jacobian"]

  Q1 -- "yes" --> H
  Q1 -- "no" --> Q2
  Q2 -- "one in, one out" --> D
  Q2 -- "many in, one out" --> G
  Q2 -- "many out" --> J

Each generates a function that keeps the original parameters and adds to them – a pointer per differentiated parameter in reverse mode, a result matrix in Hessian and Jacobian mode. Which parameters are differentiated is chosen with the second argument, so clad::gradient(f, "x") generates a derivative with one extra pointer rather than two. Derived function types gives the full signature of each, and Core concepts explains why forward mode suits few inputs and reverse mode suits many.

Three variants modify that choice rather than replacing it. clad::estimate_error generates the gradient and, with it, an estimate of the floating-point error. clad::differentiate<clad::opts::vector_mode> computes the same gradient in one forward pass instead of a reverse one. clad::differentiate<clad::immediate_mode> produces a derivative usable in a constant expression.

Forward Mode Automatic Differentiation

Forward mode AD computes derivatives of all the output parameters of a program with respect to an input parameter. The input parameter with respect to which differentiation takes place is termed as independent parameter. Mathematically, the forward mode AD allows the efficient computation of columns of the jacobian matrix.

clad::differentiate provides the forward mode differentiation functionality. It takes as input, a source function and independent parameter information, parameter with respect to which differentiation should take place, and returns a clad::CladFunction object. A call to clad::differentiate tells Clad to generate a function that computes the derivatives of the source function with respect to the independent parameter. We will call the function generated by Clad that computes derivatives as derived function. clad::CladFunction is a simple wrapper over the derived function and provides convenient access to it. Generated derived function is executed by calling the .execute member function on the associated clad::CladFunction object.

An example that demonstrates the usage of clad::differentiate:

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

double fn(double x, double y) { return x * x + y * y; }

int main() {
  // Differentiate 'fn' with respect to 'x'.
  auto d_fn_1 = clad::differentiate(fn, "x");

  // Computes the derivative of 'fn' with respect to 'x' at (x, y) = (3, 4).
  std::cout << d_fn_1.execute(3, 4) << "\n"; // prints: 6
}

The independent parameter can be specified either using the parameter name or the parameter index. Parameter indexing starts from 0. Therefore for the example shown above, the two differentiation calls shown below are equivalent:

clad::differentiate(fn, "x");

and:

clad::differentiate(fn, 0);

The derived function is executed by calling the .execute method on the associated clad::CladFunction object. Clad forward mode differentiated functions returns the computed derivatives, this behaviour is unique among other types of derived functions generated by Clad.

Clad can also differentiate with respect to an array element. The following example demonstrates this:

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

double fn_arr(double* arr, int n) {
  double res = 0;
  for (int i = 0; i < n - 1; ++i)
    res += arr[i] * arr[i + 1];
  return res;
}

int main() {
  // Differentiate 'fn_arr' with respect to element '1' of the 'arr' parameter.
  auto d_fn_arr = clad::differentiate(fn_arr, "arr[1]");
  double arr[5] = {1, 2, 3, 4, 5};
  std::cout << d_fn_arr.execute(arr, 5) << "\n"; // prints: 4
}

Arbitrary higher-order derivatives can also be computed using forward mode differentiation. The following example demonstrates computation of higher-order derivatives:

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

double fn(double i) { return i * i * i * i; }

int main() {
  // Differentiate the 3rd order derivative of 'fn' with respect to 'i'.
  auto d_fn_3 = clad::differentiate<3>(fn, "i");
  std::cout << d_fn_3.execute(3) << "\n"; // prints: 72
}

Clad also differentiates a call to a function it has a custom derivative for, such as std::sin, more than once. A usage example can be something like:

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

double mysin(double x) { return std::sin(x); }

int main() {
  auto d_sin_2 = clad::differentiate<2>(mysin);
  std::cout << d_sin_2.execute(3) << "\n"; // prints: -0.14112
}

Note

This chaining currently stops at the second derivative. Asking for a third derivative through a built-in pushforward compiles and then fails to link, because clad declares the pushforward of the pushforward of the pushforward without defining it. See issue #2081.

Note

For derivative orders up to 3, Clad defines enums that can be used instead of the integer template parameter, so the third-order example above can also be written:

auto d_fn_3 = clad::differentiate<clad::order::third>(fn, "i");

Note

Forward mode AD can only be used to differentiate with respect to a single value. For differentiating with respect to multiple values (parameters), reverse-mode AD must be used..

Visit API reference of clad::differentiate for more details.

Reverse Mode Automatic Differentiation

Reverse mode AD computes the derivative of one output with respect to every input at once. Mathematically it gives a row of the jacobian matrix, where forward mode gives a column, which is why it is the mode to reach for when a function has many inputs and few outputs, which is what a cost or a likelihood usually is.

clad::gradient provides the reverse mode differentiation functionality. It takes a source function and, optionally, the parameters to differentiate with respect to, and returns a clad::CladFunction the same way clad::differentiate does. The derived function differs in its signature: it returns nothing, and takes one extra pointer per differentiated parameter, in the order the parameters are declared. Clad accumulates into those pointers rather than assigning to them, so the caller allocates them and sets them to zero.

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

double fn(double x, double y) { return x * x + y * y; }

int main() {
  // Differentiate 'fn' with respect to every parameter.
  auto fn_grad = clad::gradient(fn);

  // The derivative is accumulated into these, so they start at zero.
  double dx = 0, dy = 0;

  // The arguments of 'fn' first, then one pointer per differentiated
  // parameter, in the same order.
  fn_grad.execute(3, 4, &dx, &dy);

  printf("dfn/dx = %g, dfn/dy = %g\n", dx, dy);
  // prints: dfn/dx = 6, dfn/dy = 8
}

Passing the parameters explicitly narrows the gradient: clad::gradient(fn, "x") generates a function taking one extra pointer rather than two. As in forward mode, a parameter can be named or given by index, so clad::gradient(fn, "x, y") and clad::gradient(fn, "0, 1") are the same request.

Visit the API reference of gradient() for more details, and Core concepts for what Clad generates and why it needs to store values along the way.

Hessian Computation

Clad can directly compute the hessian matrix of a function using the clad::hessian function.

For parameters \(x_1, x_2, \ldots, x_n\) it is the matrix of second derivatives:

\[\begin{split}\mathbf{H}_f = \begin{bmatrix} \pdv[2]{f}{x_1} & \pdv{f}{x_1}{x_2} & \cdots & \pdv{f}{x_1}{x_n} \\[6pt] \pdv{f}{x_2}{x_1} & \pdv[2]{f}{x_2} & \cdots & \pdv{f}{x_2}{x_n} \\[6pt] \vdots & \vdots & \ddots & \vdots \\[6pt] \pdv{f}{x_n}{x_1} & \pdv{f}{x_n}{x_2} & \cdots & \pdv[2]{f}{x_n} \end{bmatrix}\end{split}\]

clad::hessian provides the hessian computation functionality. The clad::hessian function takes a source function as input, and optionally, information about independent parameters – with respect to which hessian should be computed.

Internally, clad::hessian uses both the forward mode AD and the reverse mode AD to efficiently compute hessian matrix.

An example that demonstrates the usage of clad::hessian:

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

double kinetic_energy(double mass, double velocity) {
  return 0.5 * mass * velocity * velocity;
}

int main() {
  // Tells clad to generate a function that computes the hessian matrix of
  // 'kinetic_energy' with respect to all the input parameters.
  auto hessian_one = clad::hessian(kinetic_energy);

  // The independent arguments can also be named explicitly.
  auto hessian_two = clad::hessian(kinetic_energy, "mass, velocity");

  // A matrix per call to store the hessian in. Each must have enough space:
  // 2 independent variables require 4 elements (2 * 2 = 4). clad adds into
  // the matrix, so each one starts at zero.
  double matrix_one[4] = {0}, matrix_two[4] = {0};

  // Substitutes these values into the hessian function and pipes the result
  // into the matrix.
  hessian_one.execute(10, 2, matrix_one);
  printf("%g %g %g %g\n", matrix_one[0], matrix_one[1], matrix_one[2],
         matrix_one[3]);
  // prints: 0 2 2 10

  hessian_two.execute(5, 1, matrix_two);
  printf("%g %g %g %g\n", matrix_two[0], matrix_two[1], matrix_two[2],
         matrix_two[3]);
  // prints: 0 1 1 5
}

Few important things to note about clad::hessian:

  • If no independent variable information is provided, then hessian matrix is computed by taking all differentiable function parameters as independent variables.

  • Independent argument information is provided as a string literal with comma separated names of function parameters. For example:

    double product(double i, double j, double k) {
  return i*j*k;
}
// Tells Clad to generate a function that computes hessian matrix of the
// 'product' function with respect to the input parameters 'i' and 'j'.
auto d_fn = clad::hessian(product, "i, j");

  • clad::hessian also supports differentiating w.r.t multiple parameters.

auto d_fn = clad::hessian(fn);
double hessian_matrix[4] = {};
d_fn.execute(3, 5, hessian_matrix);

  • hessian matrix array should be passed as the last argument to the call to the CladFunction::execute as shown in the above code sample. The derived function adds into that array rather than overwriting it, so it should start zeroed. The array size should at least be as big as the size required to store the hessian matrix. Passing an array less than the required size will result in undefined behaviour.

Consider the case of the input being an array we need to specify the array index that needs to be differentiated even when we want to differentiate w.r.t entire array.

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

double fn(double x, double arr[2]) { return x * arr[0] * arr[1]; }

int main() {
  auto fn_hessian = clad::hessian(fn, "x, arr[0:1]");

  // We have 3 independent variables, thus we require space for 9 elements.
  double mat_fn[9] = {0};
  double num[2] = {1, 2};
  fn_hessian.execute(3, num, mat_fn);

  printf("%g %g %g\n%g %g %g\n%g %g %g\n", mat_fn[0], mat_fn[1], mat_fn[2],
         mat_fn[3], mat_fn[4], mat_fn[5], mat_fn[6], mat_fn[7], mat_fn[8]);
  // prints: 0 2 1
  // prints: 2 0 3
  // prints: 1 3 0
}

Jacobian Computation

The Jacobian matrix is the generalization of the gradient for vector-valued functions of several variables.

Clad can compute the jacobian matrix of a function through the clad::jacobian interface.

For a function with \(n\) parameters \(x_1, \ldots, x_n\) and \(m\) outputs \(f_1, \ldots, f_m\), it holds one row per output and one column per parameter:

\[\begin{split}\mathbf{J}_f = \begin{bmatrix} \pdv{f_1}{x_1} & \pdv{f_1}{x_2} & \cdots & \pdv{f_1}{x_n} \\[6pt] \pdv{f_2}{x_1} & \pdv{f_2}{x_2} & \cdots & \pdv{f_2}{x_n} \\[6pt] \vdots & \vdots & \ddots & \vdots \\[6pt] \pdv{f_m}{x_1} & \pdv{f_m}{x_2} & \cdots & \pdv{f_m}{x_n} \end{bmatrix}\end{split}\]

A self-explanatory example that demonstrates the usage of clad::jacobian:

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

void fn_jacobian(double i, double j, double* res) {
  res[0] = i * i;
  res[1] = j * j;
  res[2] = i * j;
}

int main() {
  // Generates all first-order partial derivatives columns of a jacobian matrix
  // and stores CallExprs to them inside a single function.
  auto jacobian = clad::jacobian(fn_jacobian);

  // An empty matrix to store the jacobian in. It must have enough space:
  // 5 columns (the sum of the independent variable sizes) and 3 rows (the
  // size of res).
  clad::matrix<double> d_res(3, 5);

  // Substitutes these values into the jacobian function and pipes the result
  // into the d_res variable.
  double res[3] = {0, 0, 0};
  jacobian.execute(3, 5, res, &d_res);

  // Now the derivatives are available as d_res[i][j].
  printf("%g %g\n%g %g\n%g %g\n", d_res[0][0], d_res[0][1], d_res[1][0],
         d_res[1][1], d_res[2][0], d_res[2][1]);
  // prints: 6 0
  // prints: 0 10
  // prints: 5 3
}

Few important things to note through this example:

  • clad::jacobian supports differentiating w.r.t multiple parameters.

  • The clad::matrix args are generated for all array/pointer parameters. They need to be passed after the original parameters to CladFunction::execute call. The size of every matrix needs to be exactly the size required to store the derivative of the corresponding parameter w.r.t. all input parameters. Passing a matrix of a different size will result in undefined behaviour.

Array Support

Clad currently supports differentiating arrays for forward, reverse, hessian and error estimation modes. The interface for these vary a bit.

Forward mode: The interface requires the user to provide the exact index of the array for which the function is to be differentiated. The interface of the diff function remains the same as before. An example is given below:

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

double f(double arr[4]) { return arr[0] * arr[1] * arr[2] * arr[3]; }

int main() {
  // Differentiating the function f w.r.t arr[1]:
  auto f_diff = clad::differentiate(f, "arr[1]");

  double arr[4] = {1, 2, 3, 4};
  // Pass the input to f to the execute function.
  // The output is stored in a variable with the same type as the return type
  // of f.
  double f_dx = f_diff.execute(arr);

  printf("df/darr[1] = %g\n", f_dx); // prints: df/darr[1] = 12
}

Reverse mode: The interface doesn’t require any specific index to be mentioned. The interface of the diff function requires you to pass T* for the independent variables after you pass the inputs to the original function. The T here is type of the original variable. The example below will explain it better:

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

double g(double x, double arr[2]) { return x * arr[0] + x * arr[1]; }

int main() {
  // Differentiating g w.r.t all the input variables (x, arr).
  auto g_grad = clad::gradient(g);

  double x = 2, arr[2] = {1, 2};
  // Create memory for the output of differentiation. clad adds into it, so
  // it starts at zero.
  double dx = 0, darr[2] = {0};

  // The inputs to the original function g (i.e x and arr) are passed
  // followed by the variables to store the output (i.e dx and darr).
  g_grad.execute(x, arr, &dx, darr);

  printf("dg/dx = %g \ndg/darr = { %g, %g } \n", dx, darr[0], darr[1]);
  // prints: dg/dx = 3
  // prints: dg/darr = { 2, 2 }
}

Hessian Mode: The interface requires the indexes of the array being differentiated to be mentioned explicitly even if you are trying to differentiate w.r.t the whole array. The interface of the diff function requires you to pass an an array T* after passing the inputs to the original function. The T is the return type of the original function and the size of the array should be at least the square of the number of independent variables (each index of an array is counted as one independent variable). Example:

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

double h(double x, double arr[3]) { return x * arr[0] * arr[1] * arr[2]; }

int main() {
  // Differentiating h w.r.t all the input variables (x, arr).
  // Note that the array and the indexes are explicitly mentioned even though
  // all the indexes (0, 1 and 2) are being differentiated.
  auto h_hess = clad::hessian(h, "x, arr[0:2]");

  double x = 2, arr[3] = {1, 2, 3};

  // Create memory for the hessian matrix. The minimum required size of the
  // matrix is the square of the number of independent variables. Since there
  // are 3 indexes of the array and a scalar variable, the total number of
  // independent variables is 4. clad adds into the matrix, so it starts at
  // zero.
  double mat[16] = {0};

  // The inputs to the original function h (i.e x and arr) are passed
  // followed by the output matrix.
  h_hess.execute(x, arr, mat);

  printf("hessian matrix: \n"
         "{ %g, %g, %g, %g\n"
         "  %g, %g, %g, %g\n"
         "  %g, %g, %g, %g\n"
         "  %g, %g, %g, %g }\n",
         mat[0], mat[1], mat[2], mat[3], mat[4], mat[5], mat[6], mat[7],
         mat[8], mat[9], mat[10], mat[11], mat[12], mat[13], mat[14], mat[15]);
  // prints: hessian matrix:
  // prints: { 0, 6, 3, 2
  // prints:   6, 0, 6, 4
  // prints:   3, 6, 0, 2
  // prints:   2, 4, 2, 0 }
}

Error estimation: This interface is the same as with reverse mode.

Differentiating Functors and Lambdas

Despite significant differences, differentiating functors and lambda expressions is remarkably similar to differentiating ordinary functions. Similarly, computing the hessian matrix and jacobian matrix of functors and lambda expressions is also similar to computing hessian matrix and jacobian matrix of ordinary functions.

Differentiating functors and lambdas means differentiating the call operator (operator()) member function defined by the functor and lambda type and executing the differentiated function using a reference to the functor object.

An example that demonstrates the differentiation of functors:

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

// A class type with user-defined call operator.
class Equation {
  double m_x, m_y;

public:
  Equation(double x = 0, double y = 0) : m_x(x), m_y(y) {}
  double operator()(double i, double j) { return m_x * i * j + m_y * i * j; }
  void setX(double x) { m_x = x; }
};

int main() {
  Equation E(3, 5);

  // A functor is an object of any type which has a user defined call
  // operator.
  //
  // Clad differentiation functions can directly differentiate functors.
  // Functors can be passed to clad differentiation functions in two distinct
  // ways:

  // 1) Pass by reference.
  // Differentiates 'E' with respect to parameter 'i'. The object 'E' is saved
  // in the 'CladFunction' object 'd_E'.
  auto d_E = clad::differentiate(E, "i");

  // 2) Pass as pointers.
  // Differentiates 'E' with respect to parameter 'i'. The object 'E' is saved
  // in the 'CladFunction' object 'd_E_pointer'.
  auto d_E_pointer = clad::differentiate(&E, "i");

  // Calculate the derivative of 'E' when (i, j) = (7, 9).
  std::cout << d_E.execute(7, 9) << "\n";         // prints: 72
  std::cout << d_E_pointer.execute(7, 9) << "\n"; // prints: 72
}

Functors and lambda expressions can be passed both by reference and by pointers. Therefore, the two differentiation calls shown below are equivalent:

Experiment E;  // a functor
// passing function by reference
auto d_E = clad::differentiate(E, "i");

and:

Experiment E;  // a functor
// passing function by pointer
auto d_E = clad::differentiate(&E, "i");

An example that demonstrates differentiation of lambda expressions:

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

int main() {
  auto lambda = [](double i, double j) { return i * j; };
  // Pass by reference.
  auto lambda_grad = clad::gradient(lambda);
  // Can be passed by pointer as well!
  auto lambda_grad_pointer = clad::gradient(&lambda);

  double d_i_1, d_j_1, d_i_2, d_j_2;
  d_i_1 = d_j_1 = d_i_2 = d_j_2 = 0;

  lambda_grad.execute(3, 5, &d_i_1, &d_j_1);
  lambda_grad_pointer.execute(3, 5, &d_i_2, &d_j_2);

  std::cout << d_i_1 << " " << d_j_1 << "\n"; // prints: 5 3
  std::cout << d_i_2 << " " << d_j_2 << "\n"; // prints: 5 3
}

Note

Functor class should not contain multiple overloaded call operators. This restriction will be removed in the future.

Differentiating Templates and Overloaded Functions

Clad differentiates a function, not a name. When a name stands for more than one function – a function template, or a set of overloads – the call has to say which one is meant, exactly as taking its address would.

For a function template, name the instantiation:

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

template <typename T> T poly(T x, T y) { return x * x * y; }

int main() {
  auto d_poly = clad::gradient(poly<double>);
  double d_x = 0, d_y = 0;
  d_poly.execute(3, 5, &d_x, &d_y);
  std::cout << d_x << " " << d_y << "\n"; // prints: 30 9
}

Clad differentiates that instantiation, with the template arguments already substituted, so it sees ordinary code in which T is double. Each instantiation is a different function and gets its own derivative: clad::gradient(poly<float>) derives a second one. The same holds for member functions of class templates and for explicit specializations.

For an overload set, cast to the signature you want:

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

double area(double r) { return 3.0 * r * r; }
double area(double w, double h) { return w * h; }

int main() {
  auto d_disk = clad::gradient(static_cast<double (*)(double)>(area));
  auto d_rect = clad::gradient(static_cast<double (*)(double, double)>(area));

  double d_r = 0, d_w = 0, d_h = 0;
  d_disk.execute(2, &d_r);
  d_rect.execute(3, 5, &d_w, &d_h);
  std::cout << d_r << " " << d_w << " " << d_h << "\n"; // prints: 12 5 3
}

Overloaded member functions take a pointer-to-member cast:

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

class A {
public:
  float f1(float x) { return x + x + x; }
  double f1(double x) { return x + x + x + x; }
};

int main() {
  auto d_f1 = clad::differentiate(static_cast<float (A::*)(float)>(&A::f1), 0);

  A a;
  std::cout << d_f1.execute(a, 2.0f) << "\n"; // prints: 3
}

The cast picks the overload by ordinary C++ rules, before Clad is involved. In particular, a name declared in a derived class hides the base class overloads of that name. A cast to a signature only the base declares therefore has to use the base in the type, as in static_cast<double (A::*)(double)>(&B::f1), and that is true whether or not a using declaration makes the overload visible in the derived class.

If the call does not narrow the name down to a single function, the error comes from C++ overload resolution rather than from Clad, and does not mention differentiation at all:

error: no matching function for call to 'gradient'
note: candidate template ignored: couldn't infer template argument 'F'

That message means the argument named a set of functions instead of one function: add the template argument or the cast.

Differentiable Class Types

Note

This feature is currently experimental. Expect some adventures while using it. If you do decide to go on this adventure, please consider giving your review as well on how we can improve this functionality.

One of the main goals of Clad is to be able to differentiate existing codebases with minimal boilerplate code. Most existing codebases invariably use several data structures for representing different information at various stages of computations. Thus, we are adding support for differentiating class type objects for convenient and effective differentiation of existing codebases.

Before going any further, let’s first understand how data structures and functions relate to their mathematical counterparts, and what does it mean for a data structure or a function to be differentiable. As per the principles of calculus and from purely mathematical perspective, only mathematical functions, are differentiable. Intuitively, data structures represent a mathematical space, and are thus, not differentiable from mathematical point of view. On the other hand, functions represent mathematical functions and can thus, be differentiable. Please note that not all functions represent a mathematical space. We will call a function, a differentiable function, if it represents a mathematical function and can be differentiated. We will call a type, a differentiable type, if Clad can perform calculus on the values of this type. A differentiable function should only use differentiable types in its definition.

\[F : X \longrightarrow Y\]

where

\[\begin{split}X = (x_0, x_1, x_2, \ldots) \\ Y = (y_0, y_1, y_2, \ldots)\end{split}\]

For a class type to be differentiable, it should satisfy the following rules:

  • Represent a real vector space.

  • Have a default constructor that zero initialises the object of the class. Class objects initialised by the default constructor should represent 0 tangent vector – that is, all real data members should be equal to 0.

  • Copy initialisation should perform deep copy initialisation. For example, after performing the initialisation a(b):, 1) there should be no shared resource between a and b, and 2) values of all the associated data members of a and b should be equal.

  • Assignment operator should performs deep copy. For example, after performing the assignment a = b;: there should be no shared resource between a and b; and values of all the associated data members of a and b should be equal.

In general, the type of the derivative of a variable of type YType with respect to a variable of type XType is a function of both. Intuitively, the derivative type has to be able to represent every derivative obtained by differentiating a variable y with respect to a variable x, and there is more than one of those as soon as either is an aggregate type.

In case when both y and x are built-in scalar numerical type, as your intuition probably suggests, the derivative type is also a built-in scalar numerical type. Things get more complicated when one or both the types are aggregate types. When used in a computation, aggregate types – or more generally, class types – represent a mathematical space. In the specific case when either of y or x is a built-in numerical scalar type, the derivative type can be considered to be the other type. For example, consider the following aggregate type:

struct Vector {
  double data[5];
};

If we are differentiating a variable v of type Vector with respect to a variable x of type double. Then, the derivative is of the Vector type. Mathematically, let Vector represents a vector space \(V\), then the following relations holds true:

\[\begin{split}x \in \mathbb{R} \\ v \in V \\ \pdv{v}{x} \in V\end{split}\]

If \(\pdv{v}{x}\) is stored in a variable d_v, the individual derivatives are reached as follows:

d_v.data[0];  // derivative of v.data[0] w.r.t x
d_v.data[1];  // derivative of v.data[1] w.r.t x
.. and so on ..

Similarly, in the case of differentiating a variable y of type double with respect to a variable v of type Vector, the derivative, d_v, is again of the Vector type. But the derivatives represented by the elements of d_v.data have changed. In this case, the elements of d_v.data represent derivative of y with respect to each of the elements of v.data.

If \(\pdv{y}{v}\) is stored in a variable d_v, the individual derivatives are reached as follows:

d_v.data[0];  // derivative of y w.r.t v.data[0]
d_v.data[1];  // derivative of y w.r.t v.data[1]
.. and so on ..

Currently, class type support have the following limitations:

  • Calls to member functions are not supported

  • Calls to overloaded operators are not supported

  • Initialising class type variables using non-default constructors that take non-literal arguments (non-zero derivatives) are not supported.

  • Class cannot have pointer or reference data members.

Class type support is under active development and thus, most of these limitations will be removed soon.

The non_differentiable Attribute

Occasionally, you may want to skip differentiating a specific variable or function call. For example, some variables might be used purely for logging, as constants, or as standalone metrics. Clad provides the non_differentiable annotation attribute to safely omit generating derivatives for these components.

Clad ships the CLAD_NONDIFFERENTIABLE macro (defined in clad/Differentiator/BuiltinDerivatives.h, pulled in by Differentiator.h) for this; it expands to __attribute__((annotate("non_differentiable"))).

If CLAD_NONDIFFERENTIABLE is applied to a variable, Clad skips generating a derivative counterpart for it:

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

class PointData {
public:
  double x;
  double y;
  CLAD_NONDIFFERENTIABLE double weight; // not differentiated
};

double energy(PointData p) { return p.weight * (p.x * p.x + p.y * p.y); }

int main() {
  auto d_energy = clad::gradient(energy, "p");

  PointData p{3, 4, 2}, d_p{};
  d_energy.execute(p, &d_p);

  // No derivative is generated for 'weight', so its adjoint stays zero.
  std::cout << d_p.x << " " << d_p.y << " " << d_p.weight << "\n";
  // prints: 12 16 0
}

If the attribute is applied to a function declaration, Clad refrains from producing any derivative expressions for that specific function. Instead, calls to the primal function are injected directly, behaving as if the result has a zero derivative:

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

CLAD_NONDIFFERENTIABLE double get_scaling_factor(double i, double j) {
  return i * j;
}

double compute(double i, double j) {
  // get_scaling_factor will skip differentiation completely.
  return get_scaling_factor(i, j) + i * j;
}

int main() {
  auto d_compute = clad::gradient(compute);

  double d_i = 0, d_j = 0;
  d_compute.execute(3, 5, &d_i, &d_j);

  // Only the second term contributes: the call behaves as a constant.
  std::cout << d_i << " " << d_j << "\n";
  // prints: 5 3
}

Marking a type you do not own

The attribute must sit on the declaration you want to mark, so it cannot be attached to a library type declared in a header you do not control (a stream, an allocator, a third-party handle). For those, mark the type from the outside with CLAD_NONDIFFERENTIABLE_TYPE at global scope; Clad then treats every construction of and call on that type as opaque:

CLAD_NONDIFFERENTIABLE_TYPE(ThirdParty::Handle);
// template arguments (and commas) are fine:
CLAD_NONDIFFERENTIABLE_TYPE(std::map<int, double>);

It expands to a clad::Tag specialization carrying CLAD_NONDIFFERENTIABLE. A concrete type is marked as written; a template family (e.g. every Boxed<T>) needs a partial specialization the macro cannot express – write it directly:

namespace clad {
template <class T> class CLAD_NONDIFFERENTIABLE Tag<Boxed<T>> {};
}

The standard-library primitives Clad already knows are non-differentiable (std::string, std::allocator, the stream family) are marked this way in STLBuiltins.h.

Specifying Custom Derivatives

At times Clad may be unable to differentiate your function (e.g. if its definition is in a library and the source code is not available) or an efficient/more numerically stable expression for derivatives may be known that couldn’t be computed just by applying the rules of automatic differentiation. In such cases, it is useful to be able to specify custom derivatives for the function.

Clad supports this functionality by allowing the user to specify their own custom derivatives pushforward and pullback functions in the namespace clad::custom_derivatives::. For a function FNAME one can specify:

  • a custom derivative pushforward function by defining a function FNAME_pushforward inside the namespace clad::custom_derivatives::.

  • a custom derivative pullback function by defining a function FNAME_pullback inside the namespace clad::custom_derivatives::.

When Clad will encounter a function FNAME, it will first search for a suitable custom derivative function definition within the custom_derivatives namespace. Provided no definition was found, Clad will proceed to automatically derive the function.

Pushforward and Pullback functions in Core Concepts describes what each one is and when Clad asks for it.

Note

Currently, there is no way of specifying custom derivative function for member functions. This limitation will be removed soon.

Example:

  • Suppose that you have a function my_pow(x, y) which computes x to the power of y. In this example, Clad is not able to differentiate my_pow’s body (e.g. it calls an external library or uses some non-differentiable approximation):

    double my_pow(double x, double exponent) { // something non-differentiable here... }

However, the analytical formulas of its derivatives are known, thus one can easily specify custom derivatives. Forward mode uses the pushforward, which returns the value of the function alongside the derivative:

namespace clad {
namespace custom_derivatives {
ValueAndPushforward<double, double> my_pow_pushforward(double x,
                                                       double exponent,
                                                       double d_x,
                                                       double d_exponent) {
  return {my_pow(x, exponent),
          exponent * my_pow(x, exponent - 1) * d_x +
              (my_pow(x, exponent) * ::std::log(x)) * d_exponent};
}
} // namespace custom_derivatives
} // namespace clad

Moreover, a custom gradient can be specified:

namespace clad {
namespace custom_derivatives {
void my_pow_pullback(double x, double exponent, double d_y, double* d_x,
                     double* d_exponent) {
  double t = my_pow(x, exponent - 1);
  *d_x += exponent * t * d_y;
  *d_exponent += t * x * ::std::log(x) * d_y;
}
} // namespace custom_derivatives
} // namespace clad

Whenever Clad will encounter my_pow inside a differentiated function, it will first try to find and use the provided custom derivative function before attempting to automatically differentiate it:

double f(double x, double exponent) { return my_pow(x, exponent); }

int main() {
  // Forward mode calls my_pow_pushforward.
  auto d_f = clad::differentiate(f, "x");
  std::cout << d_f.execute(3, 4) << "\n"; // prints: 108

  // Reverse mode calls my_pow_pullback.
  auto f_grad = clad::gradient(f);
  double d_x = 0, d_exponent = 0;
  f_grad.execute(3, 4, &d_x, &d_exponent);
  std::cout << d_x << " " << d_exponent << "\n"; // prints: 108 88.9876
}

Note

Clad provides custom derivatives for some mathematical functions from <cmath> by default.

Numerical Differentiation Fallback

In the cases that Clad is unable to differentiate a function by itself or cannot see the function’s definition, it will numerically differentiate the function. Clad uses the Five-Point Stencil Method with support for differentiating most scalar or array (pointer) types. For a comprehensive demo on numerically differentiating custom/user-defined types, you can checkout the following demo.

This default behavior can be disabled by passing the -DCLAD_NO_NUM_DIFF flag during the compilation of your programs. This will cause Clad to fail and error out if it encounters something non-differentiable. Another interesting bit of information that can be solicited from numerical differentiation calls is error information. Since numerical differentiation is only a way to estimate the derivative, it is essential to keep track of any associated errors. Error estimates from numerical differentiation calls can be printed to stdout using the -fprint-num-diff-errors compilation flag. This flag is overridden by the -DCLAD_NO_NUM_DIFF flag.

Error Estimation

Clad is capable of annotating a given function with floating point error estimation code using reverse mode AD. An interface similar to clad::gradient(f) is provided as follows:

clad::estimate_error(f), where f is a pointer to the function or method to be annotated with floating point error estimation code.

The function signature of the generated code is the same as from clad::gradient(f) with an extra argument at the end of type double&, which returns the total floating point error in the function by reference. For a user function double f(double, double) example usage is described below:

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

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

int main() {
  // Generate the floating point error estimation code for 'f'.
  auto df = clad::estimate_error(f);
  // Print the generated code to standard output.
  df.dump();
  // Declare the necessary variables.
  double x = 3, y = 5, d_x = 0, d_y = 0, final_error = 0;
  // Finally call execute on the generated code.
  df.execute(x, y, &d_x, &d_y, final_error);
  // After this, 'final_error' holds the floating point error in 'f'.
  std::cout << final_error << "\n"; // prints: 1.90735e-06
}

The above example generates the floating point error estimation code using an in-built taylor approximation model. However, Clad is capable of using any user defined custom model, for information on how to use your own custom model, please visit this demo. This tutorial provides a comprehensive guide on building your own custom models and understanding the working behind the error estimation framework.

Debug functionalities

The quickest look at a derivative is through the clad::CladFunction you already hold. dump() prints it and getCode() returns it as a string:

auto df = clad::differentiate(f, "x");
df.dump();

That gives the one derivative you asked for. The switches below cover the cases it does not: everything Clad generated for a translation unit, a file you can step through in a debugger, and what the analyses decided.

Derivatives are built as an AST and belong to no file of yours, so a diagnostic about one names <clad generated code> rather than a path. Clad renders them into a buffer with real lines and columns to make that possible; one buffer holds every derivative in the translation unit, so which derivative a diagnostic means comes from a note under it. All three switches below are plugin arguments, so each needs -Xclang -plugin-arg-clad -Xclang in front.

-Rclad-analysis=<name>

Reports what the named analysis could not remove from the derivative, and where. This is the analogue of -Rpass-missed, which cannot serve clad: clang’s remark machinery runs in the backend over LLVM IR and never sees a plugin working on the AST.

<clad generated code>:4:5: remark: clad keeps this value for the reverse sweep
    4 |     double _t0 = t;
      |     ^~~~~~~~~~~~~~
<clad generated code>:4:5: note: to-be-recorded analysis could not show it unused
doc.cpp:8:12: note: in the derivative of 'f' requested here
    8 |   auto g = clad::gradient(f);
      |            ^
doc.cpp:4:3: note: the value kept is the one this expression had
    4 |   t = t * t;
      |   ^

Expect three positions: the generated statement that costs, the differentiation that asked for the derivative, and the expression in your own code whose value is being kept – the half you can act on. The note giving the reason distinguishes an analysis that ran and could not prove the value dead from one that was switched off, so a remark never claims to have proved something it never ran.

Run -plugin-arg-clad -help for the analysis names this accepts.

-fgenerated-source-dir=<dir>

Writes the code Clad generates into <dir>, one file per translation unit named after it, and names that file in the debug line table. A debugger then has something to open when it stops inside a derivative.

clang++ -g -fplugin=clad.so -Xclang -plugin-arg-clad \
    -Xclang -fgenerated-source-dir=build/ doc.cpp

Without it, the only way to reach the code is -gembed-source, which puts it in the object; lldb reads that and gdb does not support it at all. When debug information is asked for and neither is in place, Clad says so once and names the flag that fits the debugger being tuned for. Giving the flag with no directory writes nothing and turns that advice off, which is the only way to: a plugin’s diagnostic belongs to no -W group.

-fdump-generated-source

Prints a derivative as text together with the position every statement occupies in it, which is how the mapping is checked without a diagnostic to hang it on.

generated-source: f_grad
  1:44-44: {
  2:5-20: double _d_t = 0.;
  3:5-20: double t = x * x;
  4:5-18: double _t0 = t;
  5:5-13: t = t * t;
  5:9-13: t * t;

Each line is line:begin-end followed by the text that starts there. A node that renders on one line reports the columns it spans, so a subexpression such as t * t is reported inside the statement containing it; one that cannot be measured that way, such as a compound statement, reports only where it begins.

This is not -fgenerate-source-file, which appends each derivative to Derivatives.cpp for reading. This one reports where each statement sits, which is what a diagnostic needs in order to point at it.

Other ways to differentiate

Beyond the four entry points above, Clad offers a vectorised forward mode, an immediate mode for constant evaluation, support for CUDA kernels, and Enzyme as an alternative reverse-mode backend.