# Solving optimization problem used for high quality denoising

The highest voted answer to this question suggests that to denoise a signal while preserving sharp transitions one should

minimize the objective function:

$$|x-y|^2 + b|f(y)|$$

where $x$ is the noisy signal, $y$ is the denoised signal, $b$ is the regularziation parameter, and $|f(y)|$ is some L1 norm penalty. Denoising is accomplished by finding the solution $y$ to this optimization problem, and $b$ depends on the noise level.

However, there is no indication of how one might accomplish that in practice as that is a problem in a very high dimensional space, especially if the signal is e.g. 10 million samples long. In practice how is this sort of problem solved computationally for large signals?

-
Are you concerned with run time? Otherwise, the iterature on how to minimize a function is quite extensive (Levenberg-Marquardt, Nelder-Mead, etc. come to mind). There are even some modified versions that are made specifically for this. –  thang Feb 2 '13 at 1:09
Actually, I have a question for the people answering below. Besides being slow, what is wrong with just something like Levenberg-Marquardt or Nelder-Mead? These are generalized optimizers, so you can even numerically approximate $f$. –  thang Feb 2 '13 at 1:13
Yes, I am concerned with run time, but thanks for pointing these methods out. –  John Robertson Feb 4 '13 at 22:06

Boyd has A Matlab Solver for Large-Scale ℓ1-Regularized Least Squares Problems. The problem formulation in there is slightly different, but the method can be applied for the problem.

Classical majorization-minimization approach also works well. This corresponds to iteratively perform soft-thresholding (for TV, clipping).

The solutions can be seen from the links. However, there are many methods to minimize these functionals by the extensive use of optimisation literature.

PS: As mentioned in other comments, FISTA will work well. Another 'very fast' algorithm family is primal-dual algorithms. You can see the interesting paper of Chambolle for an example, however there are plethora of research papers on primal-dual methods for linear inverse problem formulations.

-
What does 'primal-dual' refer to exactly? –  Mohammad Oct 4 '12 at 17:28
Mohammad, I did not implement any primal-dual algorithm for inverse problems. However, you can see an example from the link I mentioned in the answer: the paper of Chambolle. From this paper, you can see what a primal-dual algorithm means precisely. These methods provide just another (and fastly convergent) solution to inverse problems. –  Deniz Oct 5 '12 at 13:41
I thought primal dual is combinatorial optimization? How can you transform this problem generically (for a generic $f$) into that framework? –  thang Feb 2 '13 at 1:12
thang, as I mentioned before, I am not an expert on this domain. You can see the paper of Chambolle and see that how primal-dual methods can be used to solve the problems like $\ell_1$ or TV regularization. –  Deniz Feb 3 '13 at 0:04

To solve optimization problems with TV penalty, we use a recently proposed algorithm called FISTA, which has better convergence rate than conventional iterative methods, such as ASD-POCS.

-
Is it possible for you to add some more information on the algorithm, since the only reference you linked requires purchasing the article? –  Jason R Oct 2 '12 at 12:35

In the particular case where $f(y)=\|y\|_1$, the objective function can be written as

$$\|x-y\|^2 + b\|y\|_1 = \sum_i(x_i - y_i)^2 + b\sum_i |y_i|,$$

minimizing it requires to minimize each entry of the sum:

$$\hat{y_i} = argmin \{(x_i-y_i)^2 + b|y_i|\}$$

Using subdifferentials it is possible to show that the minimizer is the soft-thresholding operator with threshold $b$. That's the method proposed by Donoho and Johnstone for signal denoising. See their paper Ideal spatial adaptation by wavelet shrinkage for more details.

So in this case, I think that you don't need a more sophisticated solver to estimate your signal.

-
You have an $L^1$ norm penalty $\vert y_i\vert$ rather than a total variation penalty $\vert y_{i+1}-y_i\vert$. Is that a typo? –  John Robertson Oct 4 '12 at 16:35
In the question it says: "and |f(y)| is some L1 norm penalty", so I just plugged the $\ell_1$ norm, which is the classic case in signal denoising. But may be I'm misunderstanding the question. –  Alejandro Oct 4 '12 at 16:49
Yah, that could have been more clear. In that quote $f$ is a function on the entire signal, not necessarily a function being run on each component of the signal, that is $f$ can combine different signal samples together e.g. $f(x_0,x_1,...)=(x_1-x_0,x_2-x_1,...)$ is perfectly legitimate. –  John Robertson Oct 4 '12 at 17:05
I see. I will add that my answer if for the particular case where $f(y)$ is the $\ell_1$ norm. –  Alejandro Oct 4 '12 at 17:32

Added: if $f(x) = \ell_1(x) = \sum |x_i|$, the terms are all independent — as @Alejandro points out, you can just minimize each term by itself. It's more interesting to minimize
$\qquad\qquad$ $\| Ax - b \|_2^2 + \lambda \|x\|_1$
where $\|x\|_1$ instead of $\|x\|_2$ is intended to push many $x_i$ to 0.
The following notes are for this case. (I call the variables $x$, not $y$.)

(A year later) another name for this for the case $f(x) = \ell_1$ norm is Elastic net regularization.
Hastie et al., Elements of Statistical Learning p. 661 ff. discuss this for classification.

A fast simple way to get an approximate solution with many $x_i = 0$ is to alternate

1. minimize $\|Ax - b\|$ by plain least squares
2. shrink a.k.a. soft-threshold: set small $x_i = 0$.

This is a form of Iteratively reweighted least squares, with weights 0 or 1. I'd expect that methods in papers cited in previous answers will give better results; this is simple.

(When minimizing a sum $f() + \lambda g()$, it's a good idea to plot $f()$ and $\lambda g()$ on a log-log scale for iter 1 2 3 ... Otherwise, one term may swamp the other, and you won't even notice —
especially when they scale differently.)

-