Differentiation of Noisy Data

This Mathematica notebook implements Chartrand's algorithm for differentiating noisy data (see Chartrand, R. "Differentiation of Noisy, Nonsmooth Data", 2007).  Given a function f, we compute the approximate smoothed derivative of f as the minimizer of the functional

F(u) = α ∫|u'|dx + (1/2)∫(((|A u - f|))^2)dx

where A is the antidifferentiation operator Au(x) = (∫_0^x)u(y)dy  and α is a regularization parameter. The idea here is to simultaneously penalize candidate functions for being ``too rough'' and for having antiderivatives too far from the original function value.  We will compare this numerically to the derivative using Richardson estimators and to the derivatives of smoothed functions:

To find this minimizer numerically, we start with a couple of matrices. The idea is that we will define the function on a vector of equally spaced x locations, and then compute derivatives and antiderivatives at the midpoints of these intervals. The first matrix we need is the differentiation operator:

DMatrix[n_,h_] := SparseArray[{Band[{1,1}] -> -1/h, Band[{1,2}] -> 1/h},{n-1,n}];

DMatrix[5,h] // MatrixForm

We now need to define the antidifferentiation operator K:

KMatrix[n_,h_] := Sum[SparseArray[{Band[{i,1}] -> h/2, Band[{i,2}] -> h/2},{n-1,n}],{i,1,n-1}];

KMatrix[5,h] // MatrixForm

We will assume that our current guess for the minimizer is a list of (n-1) values named u. We construct two utility matrices E and L using these as input. By definition, L = h D^t E D, while E is the diagonal matrix with entries ((u(x_i) - u(x_{i-1}))^2 + epsilon)^(-1/2)

EMatrix[U_,eps_] := SparseArray[Band[{1,1}] -> (((U - RotateRight[U])^2 + eps)^(-1/2))[[2;;Length[U]]] ];

EMatrix[Table[U[[i]],{i,1,10}],eps] // MatrixForm

LMatrix[U_,h_,eps_] := Module[{DMat,EMat},
DMat = DMatrix[Length[U],h];
h*Transpose[DMat].EMatrix[U,eps].DMat
];

LMatrix[Table[U[[i]],{i,1,10}],h,eps] // MatrixForm

We now construct an approximation to the Hession of the functional F. This is supposed to be

K^t K + alpha L

We get

H[U_,h_,eps_,alpha_] := Transpose[KMatrix[Length[U],h]].KMatrix[Length[U],h] + alpha*LMatrix[U,h,eps];

H[Table[U[[i]],{i,1,10}],h,eps,alpha]

We also work out the right hand side of our solution: K^t (K u - f ) + alpha L u. We need to keep in mind that the vector u is defined in between the gridpoints of f. Therefore, we can't compare the integral K u to the entire list of values of f, since it has two less values. We solve this problem by passing the corresponding values of f, which is a list one shorter than the U list. The derivative of U will be defined on the same gridpoints as f.

GVector[U_,f_,h_,eps_,alpha_] := Transpose[KMatrix[Length[U],h]].(KMatrix[Length[U],h].U -f) + alpha*LMatrix[U,h,eps].U;

GVector[Table[U[[i]],{i,1,10}],Table[f[[i]],{i,1,9}],h,eps,alpha] // ColumnForm

UpdateFn[U_,f_,h_,eps_,alpha_] := LinearSolve[H[U,h,eps,alpha],-GVector[U,f,h,eps,alpha]];
Iterator[U_,f_,h_,eps_,alpha_] := U + UpdateFn[U,f,h,eps,alpha];

We now test this out on some easy data.

fTest = Table[Abs[i/100 - 1/2] - 1/2 + RandomReal[NormalDistribution[0,0.05]],{i,1,100}];
dfTest = Table[Cos[i/100 + 0.05],{i,0,100}];
ListLinePlot[fTest,PlotMarkers->Automatic,Filling->Axis]
U =(fTest - RotateRight[fTest])/0.01;
U = AppendTo[U,0];
ListLinePlot[U,PlotMarkers->Automatic,Filling->Axis]

Manipulate[
Its = Nest[Iterator[#, fTest,0.01,0.00001,10^(-alpha)] &, U, 40];
{ListLinePlot[{Its,dfTest},PlotMarkers->Automatic,Filling->Axis,ImageSize->800],
ListLinePlot[{KMatrix[Length[Its],0.01].Its,fTest},PlotMarkers->Automatic,Filling->Axis,ImageSize->800],
ListLinePlot[{DMatrix[Length[Its],0.01].Its},PlotMarkers->Automatic,Filling->Axis,ImageSize->800]},
{alpha,-1,7},ContinuousAction->False]

?Nest

Its = Nest[Iterator[#, fTest,0.01,0.00001,0.01] &, U, 40]

KMatrix[Length[Its],0.01].Its

DMatrix[Length[Its],0.01].Its