# Numerical Implementation: Solution for the Euler Lagrange Equation Of the Rudin Osher Fatemi (ROF) Total Variation Denoising Model

I am watching some wonderful videos on Variational Methods in Image Processing, and the presenter talks about how variational methods are used to de-blur or denoise images, as well as other applications--including "active contours." Variational image processing has been an active research program for some time. One example model the presenter discusses is the Rudin, Osher, Fatemi model for image denoising. The energy function for this minimization problem is given below.

$$E(u) = \int_{\Omega}(u - f)^2 + \lambda |\nabla u|dx$$

In the equation, $$u$$ is a function which generates some data and exhibits some properties as per the first and second terms. The presenter indicated that the Variational problem is minimized using gradient descent methods for optimization.

My question was, does anyone know of a good tutorial on numerically implementing gradient descent on this type of optimization problem? I want to understand exactly how the computer will compute the solution to this problem. Gradient Descent is easy enough to compute, but I was not sure about optimizing over the discretized infinite dimensional space of functions. The functional analysis is not the problem, just understand how the computer algorithm is coded to make this work.

I found an article by Getreuer which talks about the optimization algorithm using the Split Bregman method, which also has some C code. In looking at the code though, it is a bit hard to understand the underlying algorithm versus all of the different function calls and typedefs, and all.

I was hoping someone might know a good online tutorial or video, etc., that discusses how to implement this type of optimization routine. Thanks.

• Working on answer + code for you. – Royi Mar 27 at 6:40
• @Royi Oh that is so great. Man, I imagine so many people would find your code and explanation helpful. It is so important to find some sort of readable code or explanation of how these types of algorithms work for the sake of reproducibility and validation. Anything that you can provide would be appreciated. – krishnab Mar 27 at 7:39
• I added a code for the solution. – Royi Mar 29 at 23:41

## 1 Answer

The basic idea is to right the problem in a discrete manner.
Let's start with a 1D example.

Imagine the discrete realization of $$f$$ is given by the vector $$\boldsymbol{f} \in \mathbb{R}^{n}$$.

The problem becomes:

$$\hat{ \boldsymbol{u} } = \arg \min_{ \boldsymbol{u} } \frac{1}{2} {\left\| \boldsymbol{u} - \boldsymbol{f} \right\|}_{2}^{2} + \lambda \left\| D \boldsymbol{u} \right\|_{1}$$

Where $$D$$ is finite differences operator in Matrix form such that:

$$\left\| D \boldsymbol{u} \right\|_{1} = \sum_{i = 2}^{n} \left| \boldsymbol{u}_{i} - \boldsymbol{u}_{i - 1} \right|$$

This is basically Linear Least Squares with Total Variation Regularization.

There are may methods to solve this, for instance, see The Dual Norm of Total Variation Norm (Form of ⟨⋅,⋅⟩) By Smoothing.

Modern way of solving this is using the Alternating Direction Method of Multipliers (ADMM) with great resource on this subject would be Professor Stephen P. Boyd page on ADMM.

Assume the solution to that is given by a Black Box - 1D TV Denoising.
I will point out more than that, let's assume our solution can solve the problem above for any matrix $$D$$ (Not specifically Forward / Backward Differences operator).
We'll get to a specific solution and code for it later, but let's assume we have it.

We have in image in 2D space but we could vectorize it into a column stack form using the Vectorization Operator.
It is easy to see the Fidelity Term (The $${L}_{2}$$ term) in our optimization works well. But we have to deal with the 2nd term - The TV Operator.
The TV operator could be defined in (At least?) 2 ways in 2D.
The 2 form of the TV Operators in 2D is $$\mathbf{TV} \left( X \right) = \mathbb{R}^{m \times n} \to \mathbb{R}$$

1. 2D Anistoropic Total Variation
$$\mathbf{TV}_{I} \left( X \right) = \sum_{i = 1}^{m - 1} \sum_{j = 1}^{n - 1} \sqrt{ {\left( {X}_{i, j} - {X}_{i + 1, j} \right)}^{2} + {\left( {X}_{i, j} - {X}_{i, j + 1} \right)}^{2} } + \sum_{i = 1}^{m - 1} \left| {X}_{i, n} - {X}_{i + 1, n} \right| + \sum_{j = 1}^{n - 1} \left| {X}_{m, j} - {X}_{m, j + 1} \right|$$
2. 2D Isotropic Total Variation

$$\mathbf{TV}_{A} \left( X \right) = \sum_{i = 1}^{m - 1} \sum_{j = 1}^{m} \left| {X}_{i, j} - {X}_{i + 1, j} \right| + \sum_{i = 1}^{m} \sum_{j = 1}^{n - 1} \left| {X}_{i, j} - {X}_{i, j + 1} \right|$$

Remark
One could come up with many more variations. For instance, include the diagonals as well.

I will ignore the Isotropic version for this case (Though it is not hard to solve it as well).
The Anisotropic Case is very similar to the 1D case with on difference, we need to create 2 subtractions for each pixel.
So with some tweaking of the matrix we can write is in the same model of the 1D case.
We just need to update $$D$$ to be $$D : \mathbb{R}^{m n} \to \mathbb{R}^{2 m n - m - n}$$.

So we need a code to build the matrix $$D$$ in the 2D case and the solver to the 1D case and we're done.

Work in Progress

## Results

So, the implemented code yields as following (The reference being CVX based solution):

As can be seen above, the numerical solver can achieve the optimal solution.

The code is available at my StackExchange Mathematics Q3164164 GitHub Repository.

## Black Box Solver

Relatively simple yet effective method is An Algorithm for Total Variation Minimization and Applications.

The idea is the use of Dual Norm / Support Function which suggests that:

$${\left\| D \boldsymbol{x} \right\|}_{1} = \max_{\boldsymbol{p}} \left\{ \boldsymbol{p}^{T} D \boldsymbol{x} \mid {\left\| \boldsymbol{p} \right\|}_{\infty} \leq 1 \right\}$$

Hence our problem becomes:

$$\arg \min_{ \boldsymbol{u} } \frac{1}{2} {\left\| \boldsymbol{u} - \boldsymbol{f} \right\|}_{2}^{2} + \lambda \left\| D \boldsymbol{u} \right\|_{1} = \arg \min_{ \boldsymbol{u} } \max_{ {\left\| p \right\|}_{\infty} \leq 1 } \frac{1}{2} {\left\| \boldsymbol{u} - \boldsymbol{f} \right\|}_{2}^{2} + \lambda \boldsymbol{p}^{T} D \boldsymbol{u}$$

Using the Min Max Theorem (The objective is Convex with respect to $$\boldsymbol{u}$$ and Concave with respect to $$\boldsymbol{p}$$) one could switch the order of the Maximum and Minimum which yields:

$$\arg \max_{ {\left\| p \right\|}_{\infty} \leq 1 } \min_{ \boldsymbol{u} } \frac{1}{2} {\left\| \boldsymbol{u} - \boldsymbol{f} \right\|}_{2}^{2} + \lambda \boldsymbol{p}^{T} D \boldsymbol{u}$$

Which can be solved pretty easily with Projected Gradient Descent.

### Remarks

• Wow, this is amazing. Thank you so much for your effort. I keep seeing this TV norm minimization everywhere, but did not have a sense of the implementation. This is really really helpful. Hopefully you get a lot of good karma :). – krishnab Mar 29 at 23:49