Regularization operators with non-uniform increments in r

[luisfabib]

(From #83): With non-uniform increments in r, the regularization operator matrices (regoperator) also need to be adjusted. For example, for nonuniform spacing, the current second-order difference operator is not proportional to the second derivative.

[luisfabib]

Indeed, using the finite difference matrices as approximations of the derivate operators is not correct. In regoperator we have been computing them using finite difference formulas with weights of second order accuracy (see tables in this Wikipedia article).

While the accuracy order is more than enough for our purposes, these are not correct anymore for non-uniformly incremented vectors. Using these for calculating non-uniform derivatives leads to wrong results (see Old row in the figure below).

I have found Fornberg's algorithm to be a functional alternative which correctly estimates the finite difference weights for non-uniform and uniform vectors leading to the correct derivatives (see New row in the figure below.

Here's the reference for that paper:

B. Fornberg, "Calculation of weights in finite difference formulas", SIAM Review 40 (1998), pp. 685-691. DOI

There's even a public MATLAB function fdcoeffF for Fornberg's algorithm from:

Finite Difference Methods for Ordinary and Partial Differential Equations, Steady State and Time Dependent Problems, Randall J. LeVeque, 2007, Society for Industrial and Applied Mathematics (SIAM), Philadelphia Link

One (probably small) issue remaining is the fact that the abosule values obtained from our old implementation are different from the ones obtained using Fornberg's algorithm. This difference in absolute values will lead to different absolute values of the regularization parameters which might be confusing.

Since in the Fornberg's algorithm's case, the non-uniform and uniform case they are equal in absolute values, it would make sense to use this method even for the uniform case. This way the regularization parameter values would be more comparable.

[luisfabib]

When I implemented Fornberg's method in DeerLab, a quick run of the test suite revealed that most tests involving regularization where either broken or failing, both for uniform and non-uniform incremented r vectors.

This can be attributed to the different absolute values of the finite difference weights obtained from Fornberg's method [1] compared to tabulated weights from Keller and Pereyra [2] (currently used in DeerLab, and tabulated in Wikipedia under the wrong citation).

After some reading and consideration I realized that the current implementation of DeerLab using the Keller and Pereyra weights are not correctly implemented.

These tables were derived for uniformly spaced grids with dx = 1. For any grid spacing different from 1 and some derivative order d, these weights should be normalized via w = w/dx^d [2]. This normalization is not done correctly in DeerLab (and neither was done in DeerAnalysis).

If this is taken into account one can see that Fornberg's method yields the exact Keller and Pereyra finite difference matrices (if normalization is taken into account).

Therefore, Fornberg's method should be the way to go since not only works exactly for uniformly spaced grids but also works for non-uniformly spaced ones.

The next step is to solve the numerical problems which arise from the different absolute values in the regularization operator.

[1] B. Fornberg, Generation of Finite Difference Formulas on Arbitrarily Spaced Grids, MoC 51 (1988), 184, 699-706 Link [2] H. B. Keller and V. Pereyra, Symbolic generation of finite difference formulas, MoC 32 (1978),144, 955-971Link

[luisfabib]

Regularization operators in DeerLab are now using Fornberg's method, making regularization compatible with non-uniformly incremented r vectors.

Closed with ad46afc2a7f7fe647fedcf3a2be7e249d2474192