Partial Derivatives

Question

I am trying to write an algorithm that will perform N-Dimensional mixed partial derivatives. I have an idea of what I need to be able to achieve, but I cannot seem to come up with the correct loops/recursion that are required to realize the N-dimensional case.

Here is the pattern for the first 4 dimensions:

| 1D  wzyx  | 2D           | 3D           | 4D           |
----------------------------------------------------------
| dx (0001) | dx    (0001) | dx    (0001) | dx    (0001) |
|           | dy    (0010) | dy    (0010) | dy    (0010) |
|           | dyx   (0011) | dyx   (0011) | dyx   (0011) |
|           |              | dz    (0100) | dz    (0100) |
|           |              | dzx   (0101) | dzx   (0101) |
|           |              | dzy   (0110) | dzy   (0110) |
|           |              | dzyx  (0111) | dzyx  (0111) |
|           |              |              | dw    (1000) |
|           |              |              | dwx   (1001) |
|           |              |              | dwy   (1010) |
|           |              |              | dwyx  (1011) |
|           |              |              | dwz   (1100) |
|           |              |              | dwzx  (1101) |
|           |              |              | dwzy  (1110) |
|           |              |              | dxyzw (1111) |

The number of derivatives for each dimension (because it follows a binary pattern) is (2^dim)-1; e.g., 2^3 = 8 - 1 = 7.

The derivative that is dyx is the dx value of the adjacent points in the y dimension. That holds true for all of the mixed partials. So that dzyx is dyx of the adjacent points in the z dimension. I'm not sure if this paragraph is relevant information for the question, just thought I'd put here for completeness.

Any help pointers suggestions are welcome. The part in bold is the part I need to realize.

::EDIT::

I'm going to to try and be a bit more explicit by providing an example of what I need. This is only a 2D case but it kind of exemplifies the whole process I think.

I need help coming up with the algorithm that will generate the values in columns dx, dy, dyx, et. al.

|  X  |  Y  | f(x, y) |  dx             |  dy       | dyx               |
-------------------------------------------------------------------------
|  0  |  0  |    4    |  (3-4)/2 = -0.5 |  (3-4)/2  | (-0.5 - (-2.0))/2 |
|  1  |  0  |    3    |  (0-4)/2 = -2.0 |  (2-3)/2  | (-2.0 - (-2.0))/2 |
|  2  |  0  |    0    |  (0-3)/2 = -1.5 | (-1-0)/2  | (-1.5 - (-1.5))/2 |
|  0  |  1  |    3    |  (2-3)/2 = -0.5 |  (0-4)/2  | (-0.5 - (-0.5))/2 |
|  1  |  1  |    2    | (-1-3)/2 = -2.0 | (-1-3)/2  | (-1.5 - (-2.0))/2 |
|  2  |  1  |   -1    | (-1-2)/2 = -1.5 | (-4-0)/2  | (-1.5 - (-1.5))/2 |
|  0  |  2  |    0    | (-1-0)/2 = -0.5 |  (0-3)/2  | (-0.5 - (-0.5))/2 |
|  1  |  2  |   -1    | (-4-0)/2 = -2.0 | (-1-2)/2  | (-2.0 - (-2.0))/2 |
|  2  |  2  |   -4    |(-4--1)/2 = -1.5 |(-4--1)/2  | (-1.5 - (-1.5))/2 |

f(x, y) is unknown, only its values are known; so analytic differentiation is of no use, it must be numeric only.

Any help pointers suggestions are welcome. The part in bold is the part I need to realize.

::EDIT - AGAIN::

Started a Gist here: https://gist.github.com/1195522

Answer 1

This problem is cleanly solved by functional programming. Indeed, \partial_{xy}f is the partial derivative along x of \partial_y f.

I suppose you have a black box function (or function object) f, taking its values as a pointer to a memory buffer. Its signature is assumed to be

double f(double* x);

Now, here is a code to get the (second order finite difference) partial derivative to f:

template <typename F>
struct partial_derivative
{
    partial_derivative(F f, size_t i) : f(f), index(i) {}

    double operator()(double* x)
    {
        // Please read a book on numerical analysis to tweak this one
        static const double eps = 1e-4;

        double old_xi = x[index];
        x[index] = old_xi + eps;
        double f_plus = f(x);

        // circumvent the fact that a + b - b != a in floating point arithmetic
        volatile actual_eps = x[index];
        x[index] = old_xi - eps;
        actual_2eps -= x[index]
        double f_minus = f(x);

        return (f_plus - f_minus) / actual_2eps;
    }

private:
    F f;
    size_t index;
};

template <typename F>
partial_derivative<F> partial(F f, index i)
{
    return partial_derivative<F>(f, i);
}

Now, to compute \partial_{123}f, you do:

boost::function<double(double*)> f_123 = partial(partial(partial(f, 0), 1), 2);

If you need to compute them all:

template <typename F>
boost::function<double(double*)> mixed_derivatives(F f, size_t* i, size_t n_indices)
{
    if (n_indices == 0) return f;
    else return partial(mixed_derivatives(f, i + 1, n_indices - 1), i[0]);
}

and so you can do:

size_t indices[] = { 0, 1, 2 };
boost::function<double(double*)> df_123 
    = mixed_derivatives(f, indices, sizeof(indices) / sizeof(size_t));

Answer 2

If I understood you correctly, I think the following can work:

function partial_dev(point, dimension):
    neighbor_selector = top_bit(dimension)
    value_selector = dimension XOR neighbor_selector
    prev_point = point_before(point,neighbor_selector)
    next_point = pointafter(point,neighbor_selector)
    if value_selector == 0:
        return (f[prev_point] - f[next_point])/2
    else:
        return ( partial_dev(prev_point, value_selector) -
                 partial_dev(next_point, value_selector) )/2

The idea is: your top bit of the dimension value is the coordinate in which the before and after points are selected. If the rest of the dimension value is 0, you use the f values for the point for partial derivative calculation. If it is not, you get the partial derivative represented by the rest of the bits to calculate the values.

If you need all the values of all the dimension values calculated, then you don't need recursion at all: just use the dimension selector as an array index, where each array element contains the full value set for that dimension. The array is initialized such that vals[0][coords] = f(coords). Then you calculate vals[1], vals[2], and when calculating vals[3], you use vals[1] as the value table instead of vals[0] (because 3 = 0b11, where neighbor selector is 0b10 and value_selector is 0b01).