Greetings.

This week, the fourth post in a series on the new finite-difference framework. If you’re one of those people that prefer binge reading, the first post is here. Unless you’re watching the whole first season of Luke Cage on Netflix, that is.

In other news, Hacktoberfest 2016 has begun, in which you can contribute to QuantLib (or any other project on GitHub) and earn a t-shirt. I didn’t get one last year, since I would have felt a bit like cheating if I had opened pull requests on QuantLib. What can I say? I’m old-fashioned. However, you shouldn’t have such qualms, so code away. (And I have a couple of projects to which I’ll try to contribute, too, so I might get a t-shirt after all. And you might hear shortly about one of them…)

Follow me on Twitter if you want to be notified of new posts, or add me to your Google+ circles, or subscribe via RSS: the buttons for that are in the footer. Also, make sure to check my Training page.

Basic operators

And now, some operators. As for the old framework, the most basic would be those that represent the first and second derivative along one of the axes. In this case, “basic” doesn’t necessarily mean “simple”; the formula for a three-point first derivative when the grid \( {x_1, \ldots, x_i, \ldots, x_N} \) is not equally spaced is

\[\frac{\partial f}{\partial x}(x_i) \approx -\frac{h_+}{h_-(h_- + h_+)} f(x_{i-1}) +\frac{h_+ - h_-}{h_+ h_-} f(x_i) +\frac{h_-}{h_+(h_- + h_+)} f(x_{i+1})\]

where \( h_- = x_i - x_{i-1} \) and \( h_+ = x_{i+1} - x_i \), and the formula for the second derivative is just as long (if you’re interested, you can find the derivation of the formulas in Bowen and Smith [1]; preprints of the paper are available on the web). Furthermore, even though I’ve used the notation \( x_i \) to denote the \( i \)-th point along the relevant dimension, the points \( x_{i-1} \), \( x_i \) and \( x_{i+1} \) in the formula might actually be \( x_{i,j,k-1,l} \), \( x_{i,j,k,l} \) and \( x_{i,j,k+1,l} \) if we’re using a multi-dimensional array.

In any case, we’ll be using a three-point stencil such as the ones shown for two different directions in cases (a) and (b) of the figure below. Thus, both first- and second-derivative operators can be represented as instantiations of the same underlying class, called TripleBandLinearOp.

Stencils used for basic finite-difference derivatives.

The basic definition of triple-band operators is similar to that of tridiagonal ones, which you might remember from a previous post. In that case, when applying the operator \( T \) to an array \( x \), the \( i \)-th element of the result \( y \) could be written as:

\[y_i = l_i x_{i-1} + d_i x_i + u_i x_{i+1}\]

that is, a linear combination of the \( i \)-th element of \( x \) and of its two neighbors (\( l \), \( d \) and \( u \) are the lower diagonal, diagonal and upper diagonal). The triple-band operator generalizes the idea to multiple dimensions; here we have, for instance, that the element with indices \( i,j,k,l \) in the result is

\[y_{i,j,k,l} = l_k x_{i,j,k-1,l} + d_k x_{i,j,k,l} + u_k x_{i,j,k+1,l}\]

which is still the combination of the element at the same place in the input and of its two neighbors in a given direction, which is the same for all elements. (I’m glossing over the elements at the boundaries, but those are easily managed by ignoring their non-existent neighbors.)

When the input and result array are flattened, and thus the operator is written as a 2-dimensional matrix, we find the three bands that give the operator its name. Let’s refer again to the figure above, and assume that the horizontal direction is the first and the vertical direction is the second, and therefore that the elements are enumerated first left to right and then top to bottom. In case (a), the two neighbors are close to the center point in the flattened array, and the operator turns out to be tridiagonal; in case (b), the neighbors are separated by a whole row of elements in the flattened array and the operator has two non-null bands away from the diagonal (which is also not null).

The interface of the TripleBandLinearOp class is shown in the following listing.

class TripleBandLinearOp : public FdmLinearOp {
  public:
    TripleBandLinearOp(Size direction,
                       const shared_ptr<FdmMesher>& mesher);

    Array apply(const Array& r) const;
    Array solve_splitting(const Array& r, Real a,
                          Real b = 1.0) const;

    TripleBandLinearOp mult(const Array& u) const;
    TripleBandLinearOp multR(const Array& u) const;
    TripleBandLinearOp add(const TripleBandLinearOp& m) const;
    TripleBandLinearOp add(const Array& u) const;

    void axpyb(const Array& a, const TripleBandLinearOp& x,
               const TripleBandLinearOp& y, const Array& b);
  protected:
    TripleBandLinearOp() {}
    Size direction_;
    shared_array<Size> i0_, i2_;
    shared_array<Size> reverseIndex_;
    shared_array<Real> lower_, diag_, upper_;
    shared_ptr<FdmMesher> mesher_;
};

It contains a few methods besides the required apply. The solve_splitting method, like the solveFor method in the interface of old-style operators, is used to solve linear systems in which the operator multiplies an array of unknowns, and thus enables to write code using implicit schemes (which are usually stabler or faster than explicit schemes using apply); and other methods such as mult or add define a few operations that make it possible to perform basic algebra and to build operators based on existing pieces (I’ll leave it to you to check the code and see what’s what).

The constructor performs quite a bit of preliminary computations in order to make the operator ready to use.

TripleBandLinearOp::TripleBandLinearOp(
            Size direction,
            const shared_ptr<FdmMesher>& mesher) : /* ... */ {

   shared_ptr<FdmLinearOpLayout> layout = mesher->layout();
   vector<Size> newDim(layout->dim());
   iter_swap(newDim.begin(), newDim.begin()+direction_);
   vector<Size> newSpacing = FdmLinearOpLayout(newDim).spacing();
   iter_swap(newSpacing.begin(), newSpacing.begin()+direction_);

   for (FdmLinearOpIterator iter = layout->begin();
                            iter != layout->end(); ++iter) {
       const Size i = iter.index();
       i0_[i] = layout->neighbourhood(iter, direction, -1);
       i2_[i] = layout->neighbourhood(iter, direction,  1);
       const vector<Size>& coordinates = iter.coordinates();
       const Size newIndex =
           inner_product(coordinates.begin(), coordinates.end(),
                         newSpacing.begin(), Size(0));
       reverseIndex_[newIndex] = i;
   }
}

On the one hand, for each point in the layout, it stores in the two arrays i0_ and i2_ the flattened indices of the two neighbors in the given direction; this is done by means of the two calls to layout->neighborhood, and allows the core of the apply method to be written as a simple loop similar to:

for (Size i=0; i < layout->size(); ++i) {
    result[i] =
        lower_[i]*a[i0_[i]] + diag_[i]*a[i] + upper_[i]*a[i2_[i]];
}

where a is the array to which the operator is applied, and which translates into code the equation I wrote above for \( y_{i,j,k,l} \).

(As a side note, the loop in the apply method would seem like an obvious candidate for parallelization, and in fact we tried to speed it up by adding a #pragma omp parallel for. However, this didn’t seem to make the code any faster, and in some cases it actually hurt performance. It might be that the compiler is using parallel instructions already and making a better job when left alone. I’m an absolute beginner at this, though, so I’ll leave it to you to investigate the matter if you’re interested.)

On the other hand, the constructor fills the reverseIndex_ array, which provides support for the implementation of solve_splitting. In short, and oversimplifying a bit, the algorithm used for solving the system is the same used for tridiagonal operators and involves a double sweep of the array that yields the result in linear time. However, the algorithm is required to sweep the array in the correct order, that is, along the direction of the derivative. Referring again to the figure, looping over the usual flattened index would sweep the array horizontally, and therefore would only work in case (a) but not in case (b).

The code solves this by introducing a second flattened index. Let’s say we have a four-dimensional array with indices \( i,j,k,l \) and corresponding sizes \( M,N,P,Q \), and that we want to take the derivative along the \( k \) index. The spacing of the layout is \( (1, M, MN, MNP) \), so the triple \( (i,j,k,l) \) yields the flattened index \( I = i+M\cdot j + MN\cdot k + MNP\cdot l \). Now, what the constructor does is to create a new spacing with its first dimension swapped with the one corresponding to the derivative; in this example, it would be \( (MN, M, 1, MNP) \). This defines a new flattened index \( \tilde{I} = MN\cdot i + M\cdot j + k + MNP\cdot l \) that sweeps the array so that consecutive values of \( \tilde{I} \) map to adjacent elements in the \( k \) direction. Finally, the values of \( I \) corresponding to each \( \tilde{I} \) are saved in the reverseIndex_ array.

This allows the solve_splitting method to perform it work correctly. The two loop will enumerate all the values of \( \tilde{I} \) in order, and thus will sweep the array in the correct direction; and the reverseIndex_ array will give us the corresponding ordinary index that can be used to access the elements into the passed array.

Array TripleBandLinearOp::solve_splitting(const Array& r,
                                          Real a, Real b) const {
    // ... initializations ...
    Size rim1 = reverseIndex_[0];
    Real bet=1.0/(a*diag_[rim1]+b);
    retVal[reverseIndex_[0]] = r[rim1]*bet;
    for (Size j=1; j<=layout->size()-1; j++){
        const Size ri = reverseIndex_[j];
        tmp[j] = a*upper_[rim1]*bet;
        bet=1.0/(b+a*(diag_[ri]-tmp[j]*lower_[ri]));
        retVal[ri] = (r[ri]-a*lower_[ri]*retVal[rim1])*bet;
        rim1 = ri;
    }
    for (Size j=layout->size()-2; j>0; --j)
        retVal[reverseIndex_[j]] -=
                tmp[j+1]*retVal[reverseIndex_[j+1]];
    retVal[reverseIndex_[0]] -= tmp[1]*retVal[reverseIndex_[1]];

    return retVal;
}

Based on TripleBandLinearOp, the framework defines generic operators for the first and second derivative along a given direction; they are sketched in the listing that follows. As you see, they’re simple enough: they only declare the constructor, which fills the lower_, diag_ and upper_ arrays with the correct values (for instance, the coefficients of \( f(x_{i-1}) \), \( f(x_i) \) and \( f(x_{i+1}) \) in the equation for \( \frac{\partial f}{\partial x}(x_i) \)) based on the passed layout.

class FirstDerivativeOp : public TripleBandLinearOp {
  public:
    FirstDerivativeOp(Size direction,
                      const shared_ptr<FdmMesher>& mesher);
};

class SecondDerivativeOp : public TripleBandLinearOp {
  public:
    SecondDerivativeOp(Size direction,
                       const shared_ptr<FdmMesher>& mesher);
};

To round up the basic operators, I’ll just mention that the framework also defines a NinePointLinearOp class, managing stencils like the one in case (c) of the figure above, and its derived class SecondOrderMixedDerivativeOp, modeling the cross-derivative \( \frac{\partial^2 f}{\partial x_i \partial x_j} \) along two directions.

In next post: an example or two of real-world operators.

Bibliography

[1] M. K. Bowen and R. Smith, Derivative formulae and errors for non-uniformly spaced points. In Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, volume 461, pages 1975–1997. The Royal Society, 2005.