Solve Efficiently the 1D Total Variation Regularized Least Squares Problem (Denoising / Deblurring)

Question

How to solve a 1D Least Squares with Total Variation Regularization?
I know gradient based methods, I wonder how much faster / efficient I can get.

Answer 1

I will answer Total Variation Regularization:

$$ \arg \min_{\boldsymbol{x}} f \left( \boldsymbol{x} \right) = \arg \min_{\boldsymbol{x}} \frac{1}{2} {\left\| A \boldsymbol{x} - \boldsymbol{y} \right\|}_{2}^{2} + \lambda {\left\| D \boldsymbol{x} \right\|}_{1} $$

Now, we can use the method of Majorization Minimization to solve this:

$$ \forall t \in \mathbb{R}, \; \frac{ {t}^{2} }{ 2 \left| {t}_{k} \right| } + \frac{\left| {t}_{k} \right|}{2} \geq \left| t \right| \Leftrightarrow \forall \boldsymbol{x} \in \mathbb{R}^{n}, \; \frac{1}{2} \boldsymbol{x}^{T} {D}^{T} \Lambda_{k}^{-1} D \boldsymbol{x} + \frac{1}{2} {\left\| D \boldsymbol{x}_{k} \right\|}_{1} \geq {\left\| D \boldsymbol{x} \right\|}_{1} $$

Where:

$$ \Lambda_{k} = \operatorname{Diag} \left( \left| D \boldsymbol{x}_{k} \right| \right) $$

Hence our iteration can be built using:

$$ \boldsymbol{x}_{k + 1} = \arg \min_{\boldsymbol{x}} g \left( \boldsymbol{x}_{k} \right) $$

Where $ g \left( \boldsymbol{x}_{k} \right) = \frac{1}{2} {\left\| A \boldsymbol{x} - \boldsymbol{y} \right\|}_{2}^{2} + \frac{\lambda}{2} \boldsymbol{x}^{T} {D}^{T} \Lambda_{k}^{-1} D \boldsymbol{x} + \frac{\lambda}{2} {\left\| D \boldsymbol{x}_{k} \right\|}_{1} $

One could check that indeed:

1. $ g \left( \boldsymbol{x}_{k} \right) = f \left( \boldsymbol{x}_{k} \right) $.
2. $ g \left( \boldsymbol{y} \right) \geq f \left( \boldsymbol{y} \right) $.

Hence:

$$\begin{aligned} \boldsymbol{x}_{k + 1} & = \arg \min_{\boldsymbol{x}} g \left( \boldsymbol{x}_{k} \right) && \text{} \\ & = \arg \min_{\boldsymbol{x}} \frac{1}{2} {\left\| A \boldsymbol{x} - \boldsymbol{y} \right\|}_{2}^{2} + \frac{\lambda}{2} \boldsymbol{x}^{T} {D}^{T} \Lambda_{k}^{-1} D \boldsymbol{x} + \frac{\lambda}{2} {\left\| D \boldsymbol{x}_{k} \right\|}_{1} && \text{} \\ & = {\left( {A}^{T} A + \lambda {D}^{T} \Lambda_{k}^{-1} D \right)}^{-1} {A}^{T} \boldsymbol{y} && \text{} \\ & = {\left( {A}^{T} A \right)}^{-1} - {\left( {A}^{T} A \right)}^{-1} {D}^{T} \left( \frac{1}{\lambda} \Lambda_{k} + D {\left( {A}^{T} A \right)}^{-1} {D}^{T} \right) D {\left( {A}^{T} A \right)}^{-1} {A}^{T} \boldsymbol{y} && \text{} \end{aligned}$$

Where the last step is from the Matrix Inversion Lemma (Woodbury Identity) and prevents the case where some elements of the diagonal matrix are close to zero.

This is actually another derivation of the Iteratively Reweighted Least Squares method.

This method will be very effective for small to medium sized problems.
For large problems it is better to reformulate it as a Regularized Least Squares:

$$\begin{aligned} \boldsymbol{x}_{k + 1} & = \arg \min_{\boldsymbol{x}} g \left( \boldsymbol{x}_{k} \right) && \text{} \\ & = \arg \min_{\boldsymbol{x}} \frac{1}{2} {\left\| A \boldsymbol{x} - \boldsymbol{y} \right\|}_{2}^{2} + \frac{\lambda}{2} {\left\| {\Lambda}_{k}^{-\frac{1}{2}} D \boldsymbol{x} \right\|}_{2}^{2} \end{aligned}$$

An use a specific solver (Assuming $ A $ is sparse).

Comparing it to another efficient solver (See How to Solve Image Denoising with Total Variation Prior Using ADMM) for the case $ A = I $ (Denoising):

While the run time is very similar and comparable (ADMM is ~20% faster per iteration) one could see that the MM based method converges faster. It is no surprise as it is basically a Fixed Point method.

The code is available at my StackExchange Signal Processing Q76446 GitHub Repository (Look at the SignalProcessing\Q76446 folder).