Goal: Vector mode with the forward and reverse mode AD CoDiPack types.

Prerequisite: Tutorial 1 - Forward mode AD, Tutorial 2 - Reverse mode AD

Function: Simple vector valued function

template<typename Real>
void func(const Real* x, size_t l, Real* y) {
  y[0] = 0.0;
  y[1] = 1.0;
  for(size_t i = 0; i < l; ++i) {
    y[0] += x[i];
    y[1] *= x[i];
  }
}
 

Full code:

 
#include <codi.hpp>
#include <iostream>
 
template<typename Real>
void func(const Real* x, size_t l, Real* y) {
  y[0] = 0.0;
  y[1] = 1.0;
  for(size_t i = 0; i < l; ++i) {
    y[0] += x[i];
    y[1] *= x[i];
  }
}
 
 
void forwardVectorMode() {
 
  using Real = codi::RealForwardVec<5>;   // Step 1: Use the vector mode type
 
  Real x[5];
  Real y[2];
  x[0] = 1.0;
  x[1] = 2.0;
  x[2] = 3.0;
  x[3] = 4.0;
  x[4] = 5.0;
 
  // Step 2: Set the seeding for each vector direction
  for(size_t i = 0; i < 5; ++i) {
    x[i].gradient()[i] = 1.0;
  }
 
  func(x, 5, y);
 
  // Step 3: Get the gradients from the outputs.
  codi::Jacobian<double> jacobian(2,5);
  for(size_t i = 0; i < 5; ++i) {
    jacobian(0,i) = y[0].getGradient()[i];
    jacobian(1,i) = y[1].getGradient()[i];
  }
 
  std::cout << "Forward vector mode:" << std::endl;
  std::cout << "f(1 .. 5) = (" << y[0] << ", " << y[1] << ")" << std::endl;
  std::cout << "df/dx (1 .. 5) = \n" << jacobian << std::endl;
}
 
void reverseVectorMode() {
 
  using Real = codi::RealReverseVec<2>;   // Step 1: Use the vector mode type
  using Tape = typename Real::Tape;
 
  Real x[5];
  Real y[2];
  x[0] = 1.0;
  x[1] = 2.0;
  x[2] = 3.0;
  x[3] = 4.0;
  x[4] = 5.0;
 
  Tape& tape = Real::getTape();
  tape.setActive();
 
  for(size_t i = 0; i < 5; ++i) {
    tape.registerInput(x[i]);
  }
 
  func(x, 5, y);
 
  tape.registerOutput(y[0]);
  tape.registerOutput(y[1]);
 
  tape.setPassive();
 
  // Step 2: Set the seeding for each vector direction
  y[0].gradient()[0] = 1.0;
  y[1].gradient()[1] = 1.0;
 
  tape.evaluate();
 
  // Step 3: Get the gradients from the inputs.
  codi::Jacobian<double> jacobian(2,5);
  for(size_t i = 0; i < 5; ++i) {
    jacobian(0,i) = x[i].getGradient()[0];
    jacobian(1,i) = x[i].getGradient()[1];
  }
 
  std::cout << "Reverse vector mode:" << std::endl;
  std::cout << "f(1 .. 5) = (" << y[0] << ", " << y[1] << ")" << std::endl;
  std::cout << "df/dx (1 .. 5) = \n" << jacobian << std::endl;
 
  tape.reset();
}
 
int main(int nargs, char** args) {
  forwardVectorMode();
 
  std::cout << std::endl;
 
  reverseVectorMode();
 
  return 0;
}
 

The vector mode of AD is introduced by using matrices for the adjoint values and tangent values $\bar y$ , $\bar x$ , $\dot x$ and $\dot y$ . If the vector mode dimension is denoted by , the matrices have the dimension $\bar Y \in \R^{m \times d}$ , $\bar X \in \R^{n \times d}$ , $\dot X \in \R^{n \times d}$ and $\dot Y \in \R^{m \times d}$ . The basic idea is that derivatives for multiple directions are computed with one tape evaluation or forward mode run.

With CoDiPack, the vector mode is enabled by exchanging the type used for the gradient computation. A fixed size vector is implemented in the Direction class. The codi::Real*Vec types use this class to provide the vector mode. Only the vector dimension needs to be provided by the user, e.g. codi::RealReverseVec<5>.

This is nearly the only change one has to make when switching to the vector mode in CoDiPack. In order to access the elements of the vector mode, the Direction class offers array access operators. Gradients can also be initialized with initializer lists that are also supported by the Direction class.

The code example above highlights the necessary changes for the vector mode with the comments Step 1 to Step 3. In step 1, the CoDiPack type is defined as a vector type. Step 2 shows the seeding for the vector directions. In step 3, the Jacobian matrix is extracted.

Performance notes

CoDiPack offers only a fixed size vector mode that has to be defined at compile time. This decision against providing a vector mode that can be defined at run time was made deliberately. In the forward mode, all operations have to check if the vectors have the correct size and need to allocate new memory accordingly. This reduces the performance of a forward vector mode. In addition, if the vector size is known a priori, the compiler can optimize the forward and reverse vector mode even further by using e.g. SIMD instructions.

In the reverse mode, it is possible to introduce a run time decision on the vector size by using the custom evaluation methods described in Example 8 - Generalized adjoint vector interface and Example 2 - Custom adjoint vector evaluation