# Implementation of Total Variation Regularization Algorithm (Lagged Diffusivity Algorithm)

I am trying to compute the derivative of an experimentally-measured quantity as a function of time. The data are fairly noisy, which causes problems. For instance, using finite differences (central differences) to compute the derivative results in a derivative data set which is even noisier than the original data set. The idea ultimately is to plot the derivative of the quantity versus the original quantity to identify whether or not a power law relationship exists (we expect it to show up). Therefore, I would like to have a fairly smooth derivative if I am to have any hope of extracting something close to the actual power law relationship from the data.

I have researched this subject fairly extensively over the past couple of days. I have found several useful references (see, e.g., here, here (see Chapter 8), and here) on using Total Variation Regularization (which I have decided to try using), but I'm running into trouble when I try to actually implement an algorithm. (I'm using Mathematica 9.0, and as far as I can tell, I will need to write my own function to do this.)

From the first reference, I'm trying to solve for the minimizer $u^*$, which is the derivative of a function $f$ on some closed interval domain, of the following functional: $$F(u)=\alpha R(u)+DF(Au-f),$$ where $R(u)$ is a term that penalizes irregularity (or roughness, or noisiness) in $u$, $A$ is an antidifferentiation operator, $DF(Au-f)$ is a data fidelity term that penalizes differences between $Au$ and $f$, and $\alpha$ is an adjustable parameter to control the relative importance of the two terms.

I'm attempting to follow the discrete implementation I in the first reference, page 3. The author says that "we assume $u$ is defined on a uniform grid $\{x_i\}_0^L=\{0,\Delta x,2\Delta x,...,L\}$. Derivatives of $u$ are computed halfway between grid points as $Du(x_i+\Delta x/2)=u(x_{i+1})-u(x_i)$. This defines our $L\times(L+1)$ differentiation matrix $D$...Integrals of $u$ are likewise computed halfway between grid points, using the trapezoidal rule to obtain an $L\times(L+1)$ matrix $A$. Let $E_n$ be the diagonal matrix whose $i$th entry is $((u_n(x_i)-u_n(x_{i-1}))^2+\epsilon)^{-1/2}$ [where $\epsilon$ is some small positive number used to avoid division by zero]. Let $L_n=\Delta x D^TE_nD$, $H_n=K^tK+\alpha L_n$. The matrix $H_n$ is an approximation of the Hessian of $F$ at $u_n$. The update $s_n=u_{n+1}-u_n$ is the solution to $H_ns_n=-g_n$, where $g_n=K^T(Ku_n-f)+\alpha L_nu_n$..."

I think there are some typos in there, but this reference seems to be about the most detailed I've found on how to actually perform the algorithm. My biggest problem with this is that the dimensions don't work out, which makes me wonder if I'm not understanding what the variables are standing for. First, shouldn't $K$ be replaced by $A$ (or $A$ with $K$)? Second, the derivative definition is missing a factor of $1/\Delta x$, right? Now, on to the dimensions. For fun, let's say I have just 5 data points. (I actually have ~1500 pts.) Then, $f$ will be a vector with dimension 5x1. Similarly, I can get a $u$ vector with dimension 5x1 by using 2-pt central differencing for the interior points, and 3-pt forward and backward differencing at the two endpoints (to maintain the same order of error ($\Delta x^2$)). The $D$ matrix will have dimension 4x5, since it gives derivatives of the derivative of $f$ (or derivatives of $u$) at the midpoints of the 4 segments by considering the 5 computed values of $u$. Similarly, the $K$ (or $A$) matrix will have dimension 4x5. Matrix $E_n$ has dimension 4x4. Thus, $L_n$ has dimension 5x5, based on its definition above and consideration of the dimensions of the factors used to computed it. Likewise, $H_n$ will have dimension 5x5. The problem comes up during computation of $g_n$. $u_n$ has dimension 5x1. $K$ has dimension 4x5, so $Ku_n$ has dimension 4x1. But we are then supposed to subtract from this $f$, which has dimension 5x1, so the subtraction is undefined.

Obviously I have done something wrong somewhere. I would really appreciate it if someone could help with this problem. Thanks.