Introduction

A while back, I wrote a post about using Devito for finite difference calculations in finance. In that post, I priced a European Call using the Black Scholes model in Numpy and Devito. I wanted to revisit that calculation, and reproduce the results using the Eigen Matrix Library.

Eigen is a well known C++ library for matrix manipulation and linear algebra. It's the basis for several open source numerical libraries including Google's TensorFlow and Ceres, Stan and many others. One of the reasons for it's widespread adoption is that allows users to achieve great numerical performance from a high level interface that resembles mathematics. I won't go into great detail about the library's design in this post, but I wanted to highlight two key design features. You can find more details elsewhere.

Eigen Library Features

Expressions Templates

As I mentioned earlier, Eigen's main focus is linear algebra. In C++, this class of computations can be efficiently handled by a meta-programming technique known as Expression Templates. Operations in linear algebra map naturally to C++'s concept of operator overloading. For example, a simple vector sum expression

Vector a, b, c;
c = a + b;

can handled by custom use of the operator+() method, which implements a for-loop over vector components. However, an ineffeciency arises when would considers compound expressions, such as

Vector a, b, c;
d = a + b + c;

A naive implemenation of operator+() would introduce a temporary variable tmp, such that the expressions compiles to

Vector a, b, c;
tmp = a + b;
d = tmp + c;

This requires additional, and unnecessary, memory alocations and for-loop traversals. To avoid this problem, Eigen uses expression templates to delay the evaluation of compound expressions like the one above. Instead of evaluating the expression directly, Eigen returns a custom expression object that is used to construct parse trees at compile time. These trees are evaluated when the expression is assigned to a variable. This is an example of lazy evaluation.

Compiler Vectorization

Eigen also allows users to take advantage of modern CPU vectorization funcionality. Eigen supports many types of vectorization instructions sets such as SSE and AVX Simply passing an additional flag to the compiler activaties this funcionality automatically. See more information here.

Demo

Now on to the demo. As a reminder, we are using the finite difference technique to compute the value of a European Call in a Black Scholes model. As with the Numpy demo, we use two 1-d arrays as buffers to propogate the differential equation backwards in time from the payoff defined at option maturity. We use central differening to approximate the differentials.

In Eigen, we use ArrayXd objects. These allow component-wise operations instead of linear algebra, which we don't need for this demo.

#include <Eigen/Dense>

using namespace Eigen;

// C , Cn are the two price buffers used to propogate the Call value backwards in time
// delta, gamma and theta are intermediate differential terms
ArrayXd S(21), C(21), Cn(21), delta(19), gamma(19), theta(19);

S is the spatial discretization. After initializing C at the payoff, we propogate the difference equation backwards in time from maturity

// Propagate system backwards in time,
// apply descretized diff eqs and spatial boundary conditions
for (size_t i = 0; i < 10; i++)
{
  delta = (0.5/dx)*(C.tail<19>() - C.head<19>());
  gamma = (1/std::pow(dx,2))*(C.tail<19>() - 2*C.segment<19>(1) + C.head<19>());
  theta = -(0.5*std::pow(s,2))*square(S.segment<19>(1))*gamma - r*S.segment<19>(1)*delta + r*C.segment<19>(1);
  
  // Move C into Cn
  Cn = std::move(C);

  C.segment<19>(1) = Cn.segment<19>(1) - dt*theta;
  //spatial bc's
  C(0) = Cn(0)*(1 - r*dt);
  C(nx-1) = 2*C(nx-2) - C(nx-3);

}

We take advantage of the Eigen array move sematics to move data from one array to another. Running the code generates the value at $t_0$

C: [6.77083e-15, 2.47797e-11, 1.38987e-08, 2.31379e-06, 0.00015092,
    0.00448551, 0.0667669, 0.533585, 2.44971, 7.05424, 14.3368,
    23.2468, 32.8816, 42.7766, 52.7503, 62.7446, 72.7435, 82.7433, 
    92.7431, 100]

which matches our results from the earlier blog post.

Conclusion

In this post, I reproduced the results for a simple finite difference calculation using the Eigen library. As always, you can find my demo code at by github repo.

In future posts, I'd like the build a simple JIT compiler in manner similar to the one that it used in Devito. My long term goal is develop a pricing library that using computer algebra and JIT techniques to price derivatives from high level symbolic expression in python. Hopefully Eigen can provide a foundation for the numerical calculations in this future library.