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?

  • $\begingroup$ 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. $\endgroup$ – thang Feb 2 '13 at 1:09
  • $\begingroup$ 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$. $\endgroup$ – thang Feb 2 '13 at 1:13
  • $\begingroup$ Yes, I am concerned with run time, but thanks for pointing these methods out. $\endgroup$ – 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.

| improve this answer | |
  • $\begingroup$ What does 'primal-dual' refer to exactly? $\endgroup$ – Spacey Oct 4 '12 at 17:28
  • $\begingroup$ 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. $\endgroup$ – Deniz Oct 5 '12 at 13:41
  • $\begingroup$ I thought primal dual is combinatorial optimization? How can you transform this problem generically (for a generic $f$) into that framework? $\endgroup$ – thang Feb 2 '13 at 1:12
  • $\begingroup$ 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. $\endgroup$ – Deniz Feb 3 '13 at 0:04

To solve optimization problems with TV penalty, we use a recently proposed algorithm called Fast Gradient Based Algorithms for Constrained Total Variation Image Denoising and Deblurring Problems (FISTA), which has better convergence rate than conventional iterative methods, such as ASD-POCS.

| improve this answer | |
  • 1
    $\begingroup$ Is it possible for you to add some more information on the algorithm, since the only reference you linked requires purchasing the article? $\endgroup$ – Jason R Oct 2 '12 at 12:35
  • $\begingroup$ iew3.technion.ac.il/~becka/papers/71654.pdf $\endgroup$ – chaohuang Oct 2 '12 at 13:47
  • $\begingroup$ @JasonR, It is basically Nesterov Acceleration of the Prox operator. Really nice work. $\endgroup$ – Royi Oct 18 '19 at 13:45

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.

| improve this answer | |
  • $\begingroup$ 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? $\endgroup$ – John Robertson Oct 4 '12 at 16:35
  • $\begingroup$ 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. $\endgroup$ – Alejandro Oct 4 '12 at 16:49
  • $\begingroup$ 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. $\endgroup$ – John Robertson Oct 4 '12 at 17:05
  • $\begingroup$ I see. I will add that my answer if for the particular case where $f(y)$ is the $\ell_1$ norm. $\endgroup$ – 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.)

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.