# Automatic balancing parameter selection
for Tikhonov-TV regularization

Ali Gholami and Silvia Gazzola

###### Abstract

This paper considers large-scale linear ill-posed inverse problems whose
solutions can be represented as sums of smooth and piecewise constant
components. To solve such problems we consider regularizers consisting of two
terms that must be balanced. Namely, a Tikhonov term guarantees the smoothness
of the smooth solution component, while a total-variation (TV) regularizer
promotes blockiness of the non-smooth solution component. A scalar parameter
allows to balance between these two terms and, hence, to appropriately
separate and regularize the smooth and non-smooth components of the solution.
This paper proposes an efficient algorithm to solve this regularization
problem by the alternating direction method of multipliers (ADMM).
Furthermore, a novel algorithm for automatic choice of the balancing parameter
is introduced, using robust statistics. The proposed approach is supported by
some theoretical analysis, and numerical experiments concerned with different
inverse problems are presented to validate the choice of the balancing
parameter.

##  1 Introduction

This paper is concerned with the solution of ill-conditioned linear systems of
equations of the form

$\mathbb{d=Gm+e},$ (1)

arising from the discretization of continuous inverse problems in different
fields of applied sciences and engineering; see, e.g., [4, 25]. Here
$\mathbb{G}$ is an $M\times N$ real matrix representing a forward operator,
$\mathbb{d}$ and $\mathbb{e}$ are real vectors of lengths $M$ representing
known measured data and unknown Gaussian white noise, respectively, and the
unknown $\mathbb{m}$ is a real vector of length $N$ representing a quantity of
interest. The singular values of $\mathbb{G}$ may concentrate close to the
origin and decay to zero, making the system ill-conditioned. Both one-
dimensional (1D) and two-dimensional (2D) problems are considered in this
paper; in the 2D case, the model parameters $\mathbb{m}$ are obtained by
stacking the columns of a rectangular $N_{z}\times N_{x}$ array, so that
$N=N_{z}\cdot N_{x}$. Due to the presence of noise and the ill-conditioning of
the system, a naive (e.g., least squares) solution of (1) possibly gives a
meaningless estimate of the unknown $\mathbb{m}$. Thus we seek to determine an
approximation of $\mathbb{m}$ by regularizing problem (1), i.e., by
incorporating some a priori information about $\mathbb{m}$ into the problem
formulation; see [14].

In this paper, we assume that the solution to the problem (1) is piecewise
smooth, e.g., it consists of a smooth background while containing regions with
discontinuities. This allows us to decompose $\mathbb{m}$ as the sum of two
components $\mathbb{m}_{1}$ and $\mathbb{m}_{2}$, i.e.,

$\mathbb{m}=\mathbb{m}_{1}+\mathbb{m}_{2},$ (2)

where $\mathbb{m}_{1}$ is the piecewise constant component of $\mathbb{m}$ and
$\mathbb{m}_{2}$ is the smooth component of $\mathbb{m}$. Following [12], we
incorporate regularizarization by considering a Tikhonov regularization term
for the smooth component and a(n anisotropic) total variation (TV) term for
the piecewise constant component. A scalar balancing parameter $\beta>0$ is
included to control the amount of regularization associated to each term and,
hence, to properly separate the two complementary components of the solution.
Assuming that the value $\varepsilon=\|\mathbb{e}\|_{2}^{2}$ is available, the
pair ($\mathbb{m}_{1},\mathbb{m}_{2}$) is obtained by solving the following
Tikhonov-TV constrained optimization problem m_1,m_2∥D_1m_1∥_1+β2∥D_2 m_2∥_2^2
∥G(m_1+m_2)-d∥_2^2=ε, where $\|\cdot\|_{p}$, $p=1,2$, denotes the vector
$p$-norm, and $\mathbb{D}_{1}$ and $\mathbb{D}_{2}$ are scaled finite
difference discretizations of the gradient and partial second derivative
operators, respectively.

The Tikhonov-TV regularization method (2) has gained considerable attention in
scientific applications, for example: image denoising, inpainting and
deblurring [12], signal decomposition [15], optical flow computation [17],
signal filtering [22], seismic travel time tomography [3], and seismic full
waveform inversion [2]. This paper proposes an alternative algorithm to the
one in [12], featuring a fully automatic balancing parameter choice strategy,
which can be reliably used to solve large scale problems arising in the
applications mentioned above. Note that, in [12], $\xi=\beta^{-1}$ was used as
balancing parameter; using $\beta$ is more convenient for the approach
presented in this paper.

In this paper we solve (2) using the alternating direction method of
multipliers (ADMM) [11, 8]. This method has proved to be particularly
advantageous for dealing with large scale problems such as those arising in
geophysics and other scientific applications [6]. Specifically, ADMM splits
the original optimization problem into reduced subproblems that are solved in
sequence and whose solution is efficient and simple to implement. Furthermore,
ADMM naturally handles non-smooth terms (such as the $1$-norm term in (2)) by
employing proximal methods. Very recent research has focussed on applying ADMM
to nonlinear, nonconvex equality-constrained optimization problems, similar to
problem formulation (2). For instance, [27] analyzes the convergence of a
nonlinear objective function subject to a 2-norm squared equality constraint;
because of our approach in handling each subproblem (see Section 2), such
analysis can be adapted to our problem formulation (2).

The value of the parameter $\beta$ has a great impact on the quality of the
regularized solution for the problem at hand. The authors of [12] proposed to
determine $\beta$ via L-curve analysis, still assuming that an estimate of the
error bound $\varepsilon$ is known a priori. Specifically, an unconstrained
counterpart of (2) is considered, which involves two parameters (one for
weighting the fit-to-data and the regularization terms, and one for balancing
the two regularization terms); the problem is solved by applying both the
discrepancy principle and the L-curve (taking a range of values assigned to
$\beta$), requiring the solution of many instances of the considered
optimization problem (one for each value of $\beta$) to generate the L-curve.
Recently, the authors of [15] presented a procedure for tuning $\beta$ in
image decomposition problems, i.e., problems like (2) where $\mathbb{G}$ is
the identity matrix. However, their method assumes that a good estimate of the
minimum non-zero gradient norm of the blocky component $\mathbb{m}_{1}$ is
known a priori. Therefore, their method is limited to problems for which this
quantity can be estimated directly from the input data.

In this paper, we present an extremely simple and efficient strategy for
automatic selection of $\beta$ using robust statistics. Our basic idea is
that, since the desired model $\mathbb{m}$ consists of the blocky component
$\mathbb{m}_{1}$ and the smooth component $\mathbb{m}_{2}$, from a statistical
point of view the model gradient $\mathbb{D}_{1}\mathbb{m}$ is from a mixture
of non-Gaussian and Gaussian distributions. Specifically, the entries of
$\mathbb{D}_{1}\mathbb{m}_{1}$ are considered as anomalies/outliers in the
model gradient and therefore statistical anomaly detection tools are used to
identify them. The robust z-score [16] allows us to optimally determine the
lower limit of the anomalous entries (or the upper limit of normal/Gaussian
distributed entries) in the gradient vector. In particular, the optimum
$\beta$ is determined such that $\|\mathbb{D}_{1}\mathbb{m}_{2}\|_{\infty}$ is
equal to the minimum value of the anomalies (or, equivalently, the maximum
value of the normal entries) determined by the z-score. A simple algorithm is
presented to achieve this, which updates the value of $\beta$ at each ADMM
iteration, leading to an optimal separation of the model components. Extensive
numerical examples from imaging, tomography and geophysical inverse problems
are presented to show the excellent performance of the proposed method for
solving ill-posed inverse problems with piecewise smooth solutions.

The rest of the paper is organized as follows. In Section 2, we present more
details about the formulation and ADMM implementation of the Tikhonov-TV
regularization method (2). In Section 3, we introduce the new technique for
the selection of the balancing parameter $\beta$. Experimental results are
displayed in Section 4. Finally, some concluding remarks are proposed in
Section 5.

##  2 Solving the Tikhonov-TV regularized problem using ADMM

In this section we provide some details about the formulation and the solution
of problem (2), where the parameter $\beta$ is fixed. We start by precisely
defining the matrices $\mathbb{D}_{1}$ and $\mathbb{D}_{2}$ introduced in
Section LABEL:intro. For problems (1) formulated in 1D, we take

$\mathbb{D}_{1}:=\mathbb{D}_{1,\text{\scriptsize{1D}}}^{(N)}:=\begin{pmatrix}-1&1&&&\\\
&-1&1&&\\\ &&\ddots&\ddots&\\\ &&&-1&1\\\
&&&&0\end{pmatrix}\in\mathbb{R}^{N\times N}$ (3)

and

$\mathbb{D}_{2}:=\mathbb{D}_{2,\text{\scriptsize{1D}}}^{(N)}:=\underbrace{(\mathbb{D}_{1,\text{\scriptsize{1D}}}^{(N)})}_{=:\bar{\mathbb{D}}_{1}}\underbrace{\mathbb{D}_{1,\text{\scriptsize{1D}}}^{(N)}}_{=\mathbb{D}_{1}}=\left(\begin{array}[]{rrrrrr}1&-2&1&&&\\\
&1&-2&1&&\\\ &&\ddots&\ddots&\ddots&\\\ &&&1&-2&1\\\ &&&&1&-1\\\
&&&&&0\end{array}\right)\in\mathbb{R}^{N\times N},$ (4)

where the alternative notations $\mathbb{D}_{1,\text{\scriptsize{1D}}}^{(N)}$
and $\mathbb{D}_{2,\text{\scriptsize{1D}}}^{(N)}$ highlight the
differentiation order, the dimensionality of the problem, and the size of the
matrices. Although, for the 1D case,
$\bar{\mathbb{D}}_{1}=\mathbb{D}_{1}=\mathbb{D}_{1,\text{\scriptsize{1D}}}^{(N)}$,
we keep these notations to better match the 2D case (see (6)). Note that
$\mathbb{D}_{1}$ and $\mathbb{D}_{2}$ as expressed above are associated to a
rescaled forward first order and second order finite difference scheme,
respectively. For problems (1) formulated in 2D, denoting Kronecker products
by $\otimes$ and identity matrices of size $N$ by $\mathbb{I}_{N}$, we take

$\mathbb{D}_{1}:=\begin{pmatrix}\mathbb{D}_{1,\text{\scriptsize{1D}}}^{(N_{x})}\;&\\!\\!\\!\otimes\\!\\!\\!\\!&\mathbb{I}_{N_{z}}\\\
\mathbb{I}_{N_{x}}&\\!\\!\otimes\\!\\!\\!&\,\mathbb{D}_{1,\text{\scriptsize{1D}}}^{(N_{z})}\end{pmatrix}\in\mathbb{R}^{2N\times
N},$ (5)

which represents the discrete partial first order derivatives in the $x$ (or
horizontal) direction (matrix on the top) and in the $z$ (or vertical)
direction (matrix on the bottom). Similarly to the 1D case, we take

$\displaystyle\mathbb{D}_{2}$
$\displaystyle:=\begin{pmatrix}(\mathbb{D}_{1,\text{\scriptsize{1D}}}^{(N_{x})})\mathbb{D}_{1,\text{\scriptsize{1D}}}^{(N_{x})}&\\!\\!\otimes\\!\\!&\mathbb{I}_{N_{z}}\\\
\mathbb{I}_{N_{x}}&\\!\\!\otimes\\!\\!&(\mathbb{D}_{1,\text{\scriptsize{1D}}}^{(N_{z})})\mathbb{D}_{1,\text{\scriptsize{1D}}}^{(N_{z})}\end{pmatrix}$
$\displaystyle\;=\underbrace{\begin{pmatrix}(\mathbb{D}_{1,\text{\scriptsize{1D}}}^{(N_{x})})\otimes\mathbb{I}_{N_{z}}&\mathbb{0}\\\
\mathbb{0}&\mathbb{I}_{N_{x}}\otimes(\mathbb{D}_{1,\text{\scriptsize{1D}}}^{(N_{z})})\end{pmatrix}}_{=:\bar{\mathbb{D}}_{1}}\underbrace{\begin{pmatrix}\mathbb{D}_{1,\text{\scriptsize{1D}}}^{(N_{x})}&\\!\\!\otimes\\!\\!&\mathbb{I}_{N_{z}}\\\
\mathbb{I}_{N_{x}}&\\!\\!\otimes\\!\\!&\mathbb{D}_{1,\text{\scriptsize{1D}}}^{(N_{z})}\end{pmatrix}}_{=\mathbb{D}_{1}}\in\mathbb{R}^{2N\times
N}.$ (6)

Note that the one adopted above is not the only way of defining discrete
differential operators of the first and second order. Indeed, one may consider
alternative operators obtained varying boundary conditions and/or varying the
way first-order operators are combined [14].

In order to compute the regularized solution $\mathbb{m}$, one can solve
problem (2) for $\mathbb{m}_{1}$ and $\mathbb{m}_{2}$ separately, or
reformulate the problem to be solved directly for $\mathbb{m}$. Note that, in
general, the two components $\mathbb{m}_{1}$ and $\mathbb{m}_{2}$ are not
unique: indeed, for any constant signal $\mathbb{c}$, we have that
$\mathbb{D}_{1}\mathbb{c}=\mathbb{0}$ and
$\mathbb{D}_{2}\mathbb{c}=\mathbb{0}$, so that, if the pair ($\mathbb{m}_{1}$,
$\mathbb{m}_{2}$) is a solution of (2), then
$(\mathbb{m}_{1}\pm\mathbb{c},\mathbb{m}_{2}\mp\mathbb{c})$ is also a solution
of (2). However, the sum $\mathbb{m}=\mathbb{m}_{1}+\mathbb{m}_{2}$ is unique;
see [26]. Because of this, in this paper, we solve directly for $\mathbb{m}$.

Let us introduce the auxiliary variables

$\mathbb{e}=\mathbb{d-Gm},\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
\mathbb{g}_{1}=\mathbb{D}_{1}\mathbb{m}_{1},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\mathbb{g}_{2}=\mathbb{D}_{1}\mathbb{m}_{2}\,,$ (7)

so that problem (2) can be written in the following equivalent form
m,g_1,g_2,e∥g_1∥_1+β2∥¯D_1 g_2∥_2^2 D_1m=g_1+g_2 Gm+e=d ∥e∥_2^2=ε.

The augmented Lagrangian function associated with problem (7) is

$\displaystyle\mathcal{L}(\mathbb{m},\mathbb{g}_{1},\mathbb{g}_{2},\mathbb{e},\widehat{\boldsymbol{\lambda}}_{1},\widehat{\boldsymbol{\lambda}}_{2},\widehat{{\lambda}}_{3})=$
$\displaystyle\|\mathbb{g}_{1}\|_{1}+\frac{\beta}{2}\|\bar{\mathbb{D}}_{1}\mathbb{g}_{2}\|_{2}^{2}-\langle\widehat{\boldsymbol{\lambda}}_{1},\mathbb{D}_{1}\mathbb{m}-\mathbb{g}_{1}-\mathbb{g}_{2}\rangle-\langle\widehat{\boldsymbol{\lambda}}_{2},\mathbb{G}\mathbb{m}+\mathbb{e}-\mathbb{d}\rangle-\langle\widehat{\lambda}_{3},\|\mathbb{e}\|_{2}^{2}-\varepsilon\rangle$
$\displaystyle+\frac{\mu_{1}}{2}\|\mathbb{D}_{1}\mathbb{m}-\mathbb{g}_{1}-\mathbb{g}_{2}\|_{2}^{2}+\frac{\mu_{2}}{2}\|\mathbb{G}\mathbb{m}+\mathbb{e}-\mathbb{d}\|_{2}^{2}+\frac{\mu_{3}}{2}(\|\mathbb{e}\|_{2}^{2}-\varepsilon)^{2},$
(8)

where $\langle\cdot,\cdot\rangle$ denotes the canonical inner product in
$\mathbb{R}^{d}$ ($d=1,M,2N$),
$\widehat{\boldsymbol{\lambda}}_{1},\,\widehat{\boldsymbol{\lambda}}_{2},\,\widehat{\lambda}_{3}$
are the Lagrange multipliers, and $\mu_{1},\,\mu_{2},\,\mu_{3}>0$ are the
penalty parameters; see [18].

The $k$th ADMM iteration has the form

$\displaystyle\mathbb{m}^{k}$
$\displaystyle=\arg\min_{\mathbb{m}}\>\mathcal{L}(\mathbb{m},\mathbb{g}_{1}^{k-1},\mathbb{g}_{2}^{k-1},\mathbb{e}^{k-1},\widehat{\boldsymbol{\lambda}}_{1}^{k-1},\widehat{\boldsymbol{\lambda}}_{2}^{k-1},\widehat{\lambda}_{3}^{k-1}),$
(9a) $\displaystyle\mathbb{g}_{1}^{k}$
$\displaystyle=\arg\min_{\mathbb{g}_{1}}\>\mathcal{L}(\mathbb{m}^{k},\mathbb{g}_{1},\mathbb{g}_{2}^{k-1},\mathbb{e}^{k-1},\widehat{\boldsymbol{\lambda}}_{1}^{k-1},\widehat{\boldsymbol{\lambda}}_{2}^{k-1},\widehat{\lambda}_{3}^{k-1}),$
(9b) $\displaystyle\mathbb{g}_{2}^{k}$
$\displaystyle=\arg\min_{\mathbb{g}_{2}}\>\mathcal{L}(\mathbb{m}^{k},\mathbb{g}_{1}^{k},\mathbb{g}_{2},\mathbb{e}^{k-1},\widehat{\boldsymbol{\lambda}}_{1}^{k-1},\widehat{\boldsymbol{\lambda}}_{2}^{k-1},\widehat{\lambda}_{3}^{k-1}),$
(9c) $\displaystyle\mathbb{e}^{k}$
$\displaystyle=\arg\min_{\mathbb{e}}\>\mathcal{L}(\mathbb{m}^{k},\mathbb{g}_{1}^{k},\mathbb{g}_{2}^{k},\mathbb{e},\widehat{\boldsymbol{\lambda}}_{1}^{k-1},\widehat{\boldsymbol{\lambda}}_{2}^{k-1},\widehat{\lambda}_{3}^{k-1}),$
(9d) $\displaystyle\widehat{\boldsymbol{\lambda}}_{1}^{k}$
$\displaystyle=\widehat{\boldsymbol{\lambda}}_{1}^{k-1}-\mu_{1}(\boldsymbol{\nabla}\mathbb{m}^{k}-\mathbb{g}_{1}^{k}-\mathbb{g}_{2}^{k}),$
(9e) $\displaystyle\widehat{\boldsymbol{\lambda}}_{2}^{k}$
$\displaystyle=\widehat{\boldsymbol{\lambda}}_{2}^{k-1}-\mu_{2}(\mathbb{G}\mathbb{m}^{k}+\mathbb{e}^{k}-\mathbb{d}).$
(9f) $\displaystyle\widehat{\lambda}_{3}^{k}$
$\displaystyle=\widehat{\lambda}_{3}^{k-1}-\mu_{3}(\|\mathbb{e}^{k}\|_{2}^{2}-\varepsilon)\,.$
(9g)

By combining the linear and quadratic terms in the augmented Lagrangian
function (2) as

$-\langle\widehat{\boldsymbol{\lambda}}_{i},\mathbb{x}\rangle+\frac{\mu_{i}}{2}\|\mathbb{x}\|_{2}^{2}=\frac{\mu_{i}}{2}\|\mathbb{x}-(\frac{1}{\mu_{i}})\widehat{\boldsymbol{\lambda}}_{i}\|_{2}^{2}-\frac{\mu_{i}}{2}\|\widehat{\boldsymbol{\lambda}}_{i}\|_{2}^{2},$
(10)

and by making a change of variables
$\boldsymbol{\lambda}_{i}=(1/\mu_{i})\widehat{\boldsymbol{\lambda}}_{i},i=1,2,3$,
the ADMM iteration in (9a)-(9f) can be written in the simpler scaled form

$\displaystyle\mathbb{m}^{k}$
$\displaystyle=\arg\min_{\mathbb{m}}\\{\frac{\mu_{1}}{2}\|\mathbb{D}_{1}\mathbb{m}-\mathbb{g}_{1}^{k-1}-\mathbb{g}_{2}^{k-1}-\boldsymbol{\lambda}_{1}^{k-1}\|_{2}^{2}+\frac{\mu_{2}}{2}\|\mathbb{Gm}+\mathbb{e}^{k-1}-\mathbb{d}-\boldsymbol{\lambda}_{2}^{k-1}\|_{2}^{2}\\},$
(11a) $\displaystyle\mathbb{g}_{1}^{k}$
$\displaystyle=\arg\min_{\mathbb{g}_{1}}\\{\frac{\mu_{1}}{2}\|\mathbb{D}_{1}\mathbb{m}^{k}-\mathbb{g}_{1}-\mathbb{g}_{2}^{k-1}-\boldsymbol{\lambda}_{1}^{k-1}\|_{2}^{2}+\|\mathbb{g}_{1}\|_{1}\\},$
(11b) $\displaystyle\mathbb{g}_{2}^{k}$
$\displaystyle=\arg\min_{\mathbb{g}_{2}}\\{\frac{\mu_{1}}{2}\|\mathbb{D}_{1}\mathbb{m}^{k}-\mathbb{g}_{1}^{k}-\mathbb{g}_{2}-\boldsymbol{\lambda}_{1}^{k-1}\|_{2}^{2}+\frac{\beta}{2}\|\bar{\mathbb{D}}_{1}\mathbb{g}_{2}\|_{2}^{2}\\},$
(11c) $\displaystyle\mathbb{e}^{k}$
$\displaystyle=\arg\min_{\mathbb{e}}\\{\frac{\mu_{2}}{2}\|\mathbb{Gm}^{k}+\mathbb{e}-\mathbb{d}-\boldsymbol{\lambda}_{2}^{k-1}\|_{2}^{2}+\frac{\mu_{3}}{2}\left(\|\mathbb{e}\|_{2}^{2}-\varepsilon-\lambda_{3}^{k-1}\right)^{2}\\},$
(11d) $\displaystyle\boldsymbol{\lambda}_{1}^{k}$
$\displaystyle=\boldsymbol{\lambda}_{1}^{k-1}+\mathbb{g}_{1}^{k}+\mathbb{g}_{2}^{k}-\mathbb{D}_{1}\mathbb{m}^{k},$
(11e) $\displaystyle\boldsymbol{\lambda}_{2}^{k}$
$\displaystyle=\boldsymbol{\lambda}_{2}^{k-1}+\mathbb{d}-\mathbb{e}^{k}-\mathbb{G}\mathbb{m}^{k},$
(11f) $\displaystyle{\lambda}_{3}^{k}$
$\displaystyle={\lambda}_{3}^{k-1}+\varepsilon-\|\mathbb{e}^{k}\|_{2}^{2}\,.$
(11g)

In the following, we explain how the subproblems (11a)-(11d) can be solved
efficiently.

Regarding $\mathbb{m}$, since the associated subproblem (11a) is a
differential Tikhonov regularization problem in standard form, the optimality
conditions for (11a) lead to

$\mathbb{m}^{k}\\!=\\!\left(\mu_{1}\mathbb{D}_{1}^{T}\mathbb{D}_{1}+\mu_{2}\mathbb{G}^{T}\mathbb{G}\right)^{-1}\\!\\!\\!\left(\mu_{1}\mathbb{D}_{1}^{T}(\mathbb{g}_{1}^{k-1}\\!\\!+\mathbb{g}_{2}^{k-1}+\boldsymbol{\lambda}_{1}^{k-1})\\!+\mu_{2}\mathbb{G}^{T}(\mathbb{d}-\mathbb{e}^{k-1}\\!\\!+\boldsymbol{\lambda}_{2}^{k-1})\right)\\!.$
(12)

Although $\mathbb{m}^{k}$ has the closed form expression given above, in
practice and in a large-scale setting, depending on the properties of
$\mathbb{G}$, it may be necessary to apply an iterative solver (such as CG or
CGLS) to compute $\mathbb{m}^{k}$; additional details about this will be
provided in Section 4.

Regarding $\mathbb{g}_{1}$, since the associated subploblem (11b) is, by
definition, the proximal operator associated with the $\ell_{1}$\- norm
penalty, this can be computed as follows

$\displaystyle\mathbb{g}_{1}^{k}=\mathcal{T}_{\frac{1}{\mu_{1}}}(\mathbb{D}_{1}\mathbb{m}^{k}-\mathbb{g}_{2}^{k-1}-\boldsymbol{\lambda}_{1}^{k-1}),$
(13)

where the operator $\mathcal{T}_{\frac{1}{\mu_{1}}}$ is the soft-thresholding
operator defined component-wise as:

$[\mathcal{T}_{\frac{1}{\mu_{1}}}(\mathbb{x})]_{i}=x_{i}\leavevmode\nobreak\
\max(1-\frac{1}{\mu_{1}|x_{i}|},0).$ (14)

Regarding $\mathbb{g}_{2}$, again the associated subproblem (11c) is
differentiable, and the optimality condition for (11c) leads to

$\mathbb{g}_{2}^{k}=\left(\mathbb{I}+(\beta/\mu_{1})\bar{\mathbb{D}}_{1}^{T}\bar{\mathbb{D}}_{1}\right)^{-1}\left({\mathbb{D}_{1}}\mathbb{m}^{k}-\mathbb{g}_{1}^{k}-\boldsymbol{\lambda}_{1}^{k-1}\right).$
(15)

This is a first-order Tikhonov filter and can be calculated efficiently by
different methods including those based on the fast Fourier transform or the
discrete cosine transform; see [24] and the references therein.

Regarding $\mathbb{e}$, the associated subproblem (11d) is differentiable but
nonconvex, and the optimality condition of (11d) guarantees that
$\mathbb{e}^{k}$ satisfies

$\left(1+\frac{2\mu_{3}}{\mu_{2}}(\|\mathbb{e}^{k}\|_{2}^{2}-\varepsilon-\lambda_{3}^{k-1})\right)\mathbb{e}^{k}=\mathbb{d}-\mathbb{Gm}^{k}+\boldsymbol{\lambda}_{2}^{k-1}\,.$
(16)

The solution to this system of equations may not be unique. However, we claim
that

$\mathbb{e}^{k}=\gamma_{k}(\mathbb{e}^{k})(\mathbb{d}-\mathbb{Gm}^{k}+\boldsymbol{\lambda}_{2}^{k-1})\,,$
(17)

where the scale parameter $\gamma_{k}=\gamma_{k}(\mathbb{e}^{k})$ is the
maximum real root of the depressed cubic equation (see [7])

$\gamma^{3}\,+\,\frac{\mu_{2}-2\mu_{3}(\varepsilon+\lambda_{3}^{k-1})}{2\mu_{3}\|\mathbb{d}-\mathbb{Gm}^{k}+\boldsymbol{\lambda}_{2}^{k-1}\|_{2}^{2}}\,\gamma\,+\,\frac{-\mu_{2}}{2\mu_{3}\|\mathbb{d}-\mathbb{Gm}^{k}+\boldsymbol{\lambda}_{2}^{k-1}\|_{2}^{2}}=0\,,$
(18)

satisfies (16) and is also the global minimizer of (11d). Indeed, expression
(17) comes from the fact that $\mathbb{e}^{k}$ on the left-hand side of (16)
is premultiplied by a real scalar dependent on $\|\mathbb{e}^{k}\|_{2}^{2}$.
By substituting $\mathbb{e}^{k}$ expressed as in (17) into (16) leads to the
determination of $\gamma_{k}$ as a root of the depressed cubic equation (18).
Furthermore, plugging $\mathbb{e}^{k}$ from (17) into the objective function
(11d) (here denoted by $f(\gamma_{k})$), and letting

$E_{k}=\|\mathbb{d}-\mathbb{Gm}^{k}+\boldsymbol{\lambda}_{2}^{k-1}\|_{2}^{2},\quad
p_{k}=\frac{\mu_{2}-2\mu_{3}(\varepsilon+\lambda_{3}^{k-1})}{E_{k}},\quad
q_{k}=-\frac{\mu_{2}}{E_{k}}$

to simplify the notations, we get that

$\displaystyle f(\gamma_{k})$
$\displaystyle=\frac{\mu_{2}}{2}\|\mathbb{Gm}^{k}+\gamma_{k}(\mathbb{d}-\mathbb{Gm}^{k}+\boldsymbol{\lambda}_{2}^{k-1})-\mathbb{d}-\boldsymbol{\lambda}_{2}^{k-1}\|_{2}^{2}$
$\displaystyle\quad+\frac{\mu_{3}}{2}\left(\|\gamma_{k}(\mathbb{d}-\mathbb{Gm}^{k}+\boldsymbol{\lambda}_{2}^{k-1})\|_{2}^{2}-\varepsilon-\lambda_{3}^{k-1}\right)^{2}$
$\displaystyle=\frac{\mu_{2}E_{k}}{2}(\gamma_{k}-1)^{2}+\frac{\mu_{3}}{2}(E_{k}\gamma_{k}^{2}-\varepsilon-\lambda_{3}^{k-1})^{2}$
$\displaystyle=\mu_{3}E_{k}^{2}\gamma_{k}\left(\frac{1}{2}\gamma_{k}^{3}+\frac{\mu_{2}-2\mu_{3}(\varepsilon+\lambda_{3}^{k-1})}{2\mu_{3}E_{k}}\gamma_{k}-\frac{\mu_{2}}{\mu_{3}E_{k}}\right)+\frac{1}{2}(\mu_{2}E_{k}+\mu_{3}(\varepsilon+\lambda_{3}^{k-1})^{2})$
$\displaystyle=\mu_{3}E_{k}^{2}\gamma_{k}\left(-\frac{1}{2}\gamma_{k}^{3}+\gamma_{k}^{3}+p_{k}\gamma_{k}+q_{k}-\frac{\mu_{2}}{2\mu_{3}E_{k}}\right)+\frac{1}{2}(\mu_{2}E_{k}+\mu_{3}(\varepsilon+\lambda_{3}^{k-1})^{2})$
$\displaystyle=-\mu_{3}E_{k}^{2}\gamma_{k}\left(\frac{1}{2}\gamma_{k}^{3}+\frac{\mu_{2}}{2\mu_{3}E_{k}}\right)+\frac{1}{2}(\mu_{2}E_{k}+\mu_{3}(\varepsilon+\lambda_{3}^{k-1})^{2})\,,$
(19)

where we have used the fact that $\gamma_{k}$ solves (18). It follows that
$f(\gamma_{k})$ is minimized when evaluated at the largest root of (18).

##  3 Balancing Parameter Selection by Robust Statistics

The Tikhonov-TV regularized problem (2) introduced in Section LABEL:intro
stems from the statistical assumption that the model gradient
$\mathbb{g}=\mathbb{D}_{1}\mathbb{m}$ is a mixture of two components resulting
from different distributions: a non-Gaussian distributed (sparse) component
$\mathbb{g}_{1}=\mathbb{D}_{1}\mathbb{m}_{1}$ and a Gaussian distributed (non-
sparse) component $\mathbb{g}_{2}=\mathbb{D}_{1}\mathbb{m}_{2}$. We determine
the balancing parameter $\beta$ by using statistical tools that allow to
optimally separate these two components of the gradient. Specifically, we
assume that the non-zero entries of $\mathbb{g}_{1}$, which are associated
with jumps in the regularized solution, can be considered as anomalies (or
outliers) in the gradient vector $\mathbb{g}$. The lower value of these
anomalous entries can be optimally determined by robust statistics [20].

According to the classical z-score statistic, an element of a data set is
considered an anomaly if it falls outside a distance (e.g., 2.5 standard
deviations) from the mean. A smaller distance (e.g., 2.0 standard deviations)
can be used if the model size is small and a larger distance (e.g., 3.0 or 3.5
standard deviations) can be used if the model size is large. Therefore, its
value can be determined according to the user’s perspective; through this
paper we denote this threshold distance by ${\tau_{\tiny{\text{nrm}}}}$ and,
unless otherwise stated, we set its value to 2.5 for all numerical examples. A
main difficulty in anomaly detection is that the z-score is sensitive to
extreme values because the mean and standard deviation are sensitive to
extreme values.

The robust z-score [16] is based on robust estimators rather than the mean and
the standard deviation. A robust measure of mean is the median. For a given
vector $\mathbb{g}$, the median of $\mathbb{g}$, denoted by
$\text{median}(\mathbb{g})$, is the value such that at least half of the
entries in $\mathbb{g}$ are smaller than or equal to
$\text{median}(\mathbb{g})$, and that at least half of the entries in
$\mathbb{g}$ are larger than or equal to $\text{median}(\mathbb{g})$. Thus, if
$\mathbb{g}$ has an odd number of sorted elements, then
$\text{median}(\mathbb{g})$ is the middle value; otherwise it can be the
average of the two middle values. A robust measure of standard deviation is
known as the median absolute deviation (MAD). For a given vector $\mathbb{g}$,
the MAD of $\mathbb{g}$, denoted by $\text{MAD}(\mathbb{g})$, is given by

$\text{MAD}(\mathbb{g})=1.4826\leavevmode\nobreak\
\text{median}(\mathbb{g}-\text{median}(\mathbb{g}))\,,$ (20)

where the constant value 1.4826 is used to make the estimator consistent for
Gaussian distributions; see [19]. Other alternatives to the MAD can also be
used for estimation of the standard deviation; see, again, [19]. Using these
robust measures of the mean and standard deviation, a robust z-score is
assigned to each gradient sample, calculated as

$z_{i}=\frac{[\mathbb{g}]_{i}-\text{median}(\mathbb{g})}{\text{MAD}(\mathbb{g})}\,,$
(21)

and every element for which $|z_{i}|$ is larger than the threshold value
${\tau_{\tiny{\text{nrm}}}}$ is considered as anomaly. Consequently, we define
$\text{nrm}(\mathbb{g})$ as the vector comprising the ‘normal’ (i.e., the
‘not-anomalous’) entries of $\mathbb{g}$, i.e.,

$\text{nrm}(\mathbb{g})=\\{[\mathbb{g}]_{i}:\leavevmode\nobreak\
|z_{i}|\leq{\tau_{\tiny{\text{nrm}}}}\\}.$ (22)

According to the insight given at the beginning of this section, the optimal
balancing parameter to be used in (2) should be such that the smooth
components $\mathbb{g}_{2}$ of the computed solution gradient $\mathbb{g}$ are
classified as ‘normal’ (22) according to the robust z-score (21). We impose
this condition by defining $\beta_{\text{\scriptsize{opt}}}$ as a root of the
equation

$0=\|\mathbb{g}_{2}(\beta)\|_{\infty}-\|\text{nrm}(\mathbb{g}(\beta))\|_{\infty}=:\phi(\beta)\,,\quad\mbox{i.e.,}\quad\phi(\beta_{\text{\scriptsize{opt}}})=0\,.$
(23)

In the above equation, the notations $\mathbb{g}_{2}(\beta)$ and
$\mathbb{g}(\beta)$ have been employed to highlight the dependence of the
vectors $\mathbb{g}_{2}$ and $\mathbb{g}$ on $\beta$; in the following, both
these notations will be freely used depending on the context. In (23), the
vectors $\mathbb{g}_{2}$ and $\text{nrm}(\mathbb{g})$ are evaluated in the max
norm to enforce that the value of the maximum entry in the smooth and normal
components of $\mathbb{g}$ coincide.

######  Proposition 3.1.

The function $\phi(\beta)$ defined in (23) has at least one positive root.

###### Proof.

Once the ADMM iterations (11a)-(11g) have converged, we can write the computed
$\mathbb{g}_{1}$ and $\mathbb{g}_{2}$ as

$\displaystyle\mathbb{g}_{1}$ $\displaystyle=$
$\displaystyle\mathcal{T}_{\frac{1}{\mu_{1}}}(\mathbb{D}_{1}\mathbb{m}-\mathbb{g}_{2}-\boldsymbol{\lambda}_{1})$
(24) $\displaystyle\mathbb{g}_{2}$ $\displaystyle=$
$\displaystyle\left(\mathbb{I}+(\beta/\mu_{1})\bar{\mathbb{D}}_{1}^{T}\bar{\mathbb{D}}_{1}\right)^{-1}\left({\mathbb{D}_{1}}\mathbb{m}-\mathbb{g}_{1}-\boldsymbol{\lambda}_{1}\right)$
(25)

where the updates rules (13) and (15) have been exploited. From the above
equalities it can be clearly seen that, when $\beta=0$,
$\mathbb{g}_{2}(0)=\mathbb{D}_{1}\mathbb{m}-\mathbb{g}_{1}-\boldsymbol{\lambda}_{1}$,
so that

$\mathbb{g}_{1}=\mathcal{T}_{\frac{1}{\mu_{1}}}(\mathbb{D}_{1}\mathbb{m}-\mathbb{g}_{2}-\mathbf{\lambda}_{1})=\mathcal{T}_{\frac{1}{\mu_{1}}}(\mathbb{g}_{1})\,,$

which holds if and only if $\mathbb{g}_{1}=\mathbf{0}$. As a consequence,
$\mathbb{g}(0)=\mathbb{g}_{2}(0)$. After defining the diagonal matrix
$\mathbb{D}(\mathbb{g}_{2})$ whose diagonal entries are

$[\mathbb{D}(\mathbb{g}_{2})]_{i,i}=\begin{cases}1&\mbox{if
$[\mathbb{g}_{2}]_{i}\in\text{norm}(\mathbb{g})$}\\\
0&\mbox{otherwise}\end{cases}\,,$

and using standard norm inequalities, we get

$\|\text{nrm}(\mathbb{g}(0))\|_{\infty}=\|\mathbb{D}(\mathbb{g}_{2})\mathbb{g}_{2}(0)\|_{\infty}\leq\|\mathbb{D}(\mathbb{g}_{2})\|_{\infty}\|\mathbb{g}_{2}(0)\|_{\infty}=\|\mathbb{g}_{2}(0)\|_{\infty}\,,$

so that

$\phi(0)=\|\mathbb{g}_{2}(0)\|_{\infty}-\|\text{nrm}(\mathbb{g}_{2}(0))\|_{\infty}\geq
0\,.$

Let us denote by $\mathbb{v}_{i}$ and $\sigma_{i}$, $i=1,\dots,N$ the
eigenvectors and nonnegative eigenvalues of
$\bar{\mathbb{D}}_{1}^{T}\bar{\mathbb{D}}_{1}$, respectively. It follows that
$\mathbb{g}_{2}$ can be expressed as

$\mathbb{g}_{2}=\sum_{i=1}^{N}\left(1+\frac{\beta}{\mu_{1}}\sigma_{i}\right)^{-1}(\mathbb{v}_{i}^{T}(\mathbb{D}_{1}\mathbb{m}-\mathbb{g}_{1}-\boldsymbol{\lambda}_{1}))\mathbb{v}_{i},$

so that, when $\beta\rightarrow+\infty$,
$\mathbb{g}_{2}(\beta)\rightarrow\mathbf{0}$ and
$\|\mathbb{g}_{2}(\beta)\|_{\infty}\rightarrow 0$. As a consequence,

$\phi(\beta)=\|\mathbb{g}_{2}(\beta)\|_{\infty}-\|\text{nrm}(\mathbb{g}(\beta))\|_{\infty}\rightarrow-\|\text{nrm}(\mathbb{g}(\beta))\|_{\infty}\leq
0\quad\mbox{as}\quad\beta\rightarrow+\infty\,.$

The result follows from the continuity of $\phi(\beta)$, applying the
intermediate value theorem. ∎

Under the reasonable assumption that
$\|\text{nrm}(\mathbb{g}(\beta))\|_{\infty}$ is an increasing function of
$\beta$, it is possible to prove that $\phi(\beta)$ defined in (23) is a
strictly decreasing function of $\beta>0$ and, therefore, the zero of
$\phi(\beta)$ is unique; see the arguments in Appendix A for more details.

In order to solve (23), we first reformulate it as an equivalent fixed point
problem of the form

$\beta\left(\|\text{nrm}(\mathbb{g}(\beta))\|_{\infty}+\|\mathbb{g}_{2}(\beta)\|_{\infty}\right)=\beta\left(2\|\mathbb{g}_{2}(\beta)\|_{\infty}\right)\,.$
(26)

Then, starting from an initial value $\beta^{0}>0$, we solve (26) via an
iteration of the form

$\beta^{j+1}=\frac{1}{2}\beta^{j}+\frac{1}{2}\underbrace{\left(\frac{4\|\mathbb{g}_{2}(\beta^{j})\|_{\infty}}{\|\mathbb{g}_{2}(\beta^{j})\|_{\infty}+\|\text{nrm}(\mathbb{g}(\beta^{j}))\|_{\infty}}-1\right)\beta^{j}}_{=:\psi(\beta^{j})},\quad
j=0,1,\dots\,.$ (27)

Since the above update rule can be regarded as an averaged iteration algorithm
and since, under mild assumptions, it can be proved that the function
$\psi(\beta)$ is not expansive (see Appendix A), thanks to the theory
presented in [21] the iteration (27) is guaranteed to globally converge to a
fixed point of (26) and, therefore, to a solution to (23).

Within the ADMM scheme (11a)-(11e) we dynamically tune the value of $\beta$
according to rule (27), in such a way that condition (23) is (approximately)
satisfied for the final regularized solution (i.e., at the end of the ADMM
iterations). Namely, suppose that, at iteration $k$, we have access to the
balancing parameter $\beta^{k-1}$ and to the quantities
$\mathbb{g}_{1}^{k-1}$, $\mathbb{g}_{2}^{k-1}$, $\mathbb{e}^{k-1}$,
$\boldsymbol{\lambda}_{1}^{k-1}$, $\boldsymbol{\lambda}_{2}^{k-1}$. Then:
$\mathbb{m}^{k}$ is computed from $\mathbb{g}_{1}^{k-1}$,
$\mathbb{g}_{2}^{k-1}$, $\mathbb{e}^{k-1}$, $\boldsymbol{\lambda}_{1}^{k-1}$,
$\boldsymbol{\lambda}_{2}^{k-1}$; $\mathbb{g}_{1}^{k}$ is computed from
$\mathbb{m}^{k}$, $\mathbb{g}_{2}^{k-1}$, $\boldsymbol{\lambda}_{1}^{k-1}$;
$\mathbb{g}_{2}^{k}$ is computed from $\mathbb{m}^{k}$, $\mathbb{g}_{1}^{k}$,
$\boldsymbol{\lambda}_{1}^{k-1}$ and $\beta^{k-1}$. We set the value of
$\beta^{k}$ by applying one iteration of the update rule (27) and using the
updated quantities available at iteration $k$ (i.e.,
$\mathbb{g}_{2}^{k}(\beta^{k-1})$ and $\mathbb{g}^{k}(\beta^{k-1})$). Although
a convergence proof for the ADMM scheme (9a)-(9f), with $\beta$ updated
according to (27) is outside the scope of this paper, numerical experiments
consistently show that both ADMM applied with the fixed value
$\beta=\beta_{\text{\scriptsize{opt}}}$ obtained when (27) has converged and
ADMM applied with adaptive $\beta$ selected by (27) converge to the same
solution; see also Section 4. The resulting method is summarized in Algorithm
1.

Inputs: | penalty parameters $\mu_{1},\,\mu_{2},\,\mu_{3}>0$ for the Lagrangian, $\varepsilon=\|\mathbb{e}\|_{2}^{2}$, an initial value $\beta^{0}$,
---|---
| threshold $\tau_{\tiny\text{nrm}}>0$ for (22)

Set $\boldsymbol{\lambda}_{1}^{0}=\mathbb{0}$,
$\boldsymbol{\lambda}_{2}^{0}=\mathbb{0}$, ${\lambda}_{{3}}^{0}=0$,
$\mathbb{g}_{1}^{0}=\mathbb{0}$, $\mathbb{g}_{2}^{0}=\mathbb{0}$.

for  _$k=1,2,\cdots$ until a stopping criterion is satisfied _ do

Compute $\mathbb{m}^{k}$ using (12)

Compute $\mathbb{g}_{1}^{k}$ using (13)

Compute $\mathbb{g}_{2}^{k}$ using (15) with $\beta=\beta^{k-1}$

Determine $\gamma_{k}$ by solving (18) and compute $\mathbb{e}^{k}$ using (17)

Update
$\boldsymbol{\lambda}_{1}^{k}=\boldsymbol{\lambda}_{1}^{k-1}+\mathbb{g}_{1}^{k}+\mathbb{g}_{2}^{k}-\mathbb{D}_{1}\mathbb{m}^{k}$

Update
$\boldsymbol{\lambda}_{2}^{k}=\boldsymbol{\lambda}_{2}^{k-1}+\mathbb{d}-\mathbb{e}^{k}-\mathbb{G}\mathbb{m}^{k}$

Update
${\lambda}_{3}^{k}={\lambda}_{3}^{k-1}+\varepsilon-\|\mathbb{e}^{k}\|_{2}^{2}$

Update
$\beta^{k}=\dfrac{1}{2}\beta^{k-1}+\dfrac{1}{2}\left(\dfrac{4\|\mathbb{g}_{2}^{k}(\beta^{k-1})\|_{\infty}}{\|\mathbb{g}_{2}^{k}(\beta^{k-1})\|_{\infty}+\|\text{nrm}(\mathbb{g}^{k}(\beta^{k-1}))\|_{\infty}}-1\right)\beta^{k-1}$

end for

Algorithm 1 ADMM for solving the Tikhonov-TV regularized problem (7), with
authomatic update of $\beta$.

##  4 Numerical Experiments

In this section, we assess the performance of the proposed automated ADMM
algorithm in solving inverse problems arising in different applications.

Concerning the inputs of Algorithm 1, for all the experiments we set
$\mu_{1}>\mu_{2}$ and $\mu_{2}=\mu_{3}=1$. Moreover, we assume that the value
$\varepsilon=\|\mathbb{e}\|_{2}^{2}$ is known in all but one experiments,
where we test the robustness of the proposed solver with respect to
inaccuracies in the value of $\varepsilon$; similarly, we take
$\tau_{\text{\tiny{nrm}}}=2.5$ in all but one experiments, where we test the
robustness of the proposed solver with respect to the choice of
$\tau_{\text{\tiny{nrm}}}$. The iterations are stopped either when a maximum
number of iterations is reached or when the computed $\mathbb{m}$ stabilizes,
i.e.,
$\nicefrac{{\|\mathbb{m}^{k+1}-\mathbb{m}^{k}\|_{2}}}{{\|\mathbb{m}^{k}\|_{2}}}<10^{{-4}}$.

We compare the results obtained by applying the new method to the results
obtained applying TV regularization and Tikhonov regularization (with a
regularization term of the form
$\|\bar{\mathbb{D}}_{1}\mathbb{D}_{1}\mathbb{m}\|_{2}^{2}$). Both TV and
Tikhonov regularization are implemented through ADMM and can be regarded as
special cases of Algorithm 1. Specifically, to recover the TV formulation, we
set $\mathbb{g}_{2}=\mathbf{0}$, so that minimization in (7) happens on
$\mathbb{m}$, $\mathbb{g}_{1}$ and $\mathbb{e}$ only, and the first constraint
reduces to $\mathbb{D}_{1}\mathbb{m}=\mathbb{g}_{1}$; as a consequence, the
updates (12), (13), (16) and (11e)-(11g) are performed taking
$\mathbb{g}_{2}^{k-1}=\mathbb{g}_{2}^{k}=\mathbf{0}$, and update (15) is
discarded. Similarly, to recover the Tikhonov formulation, we set
$\mathbb{g}_{1}=\mathbf{0}$, so that minimization in (7) happens on
$\mathbb{m}$, $\mathbb{g}_{2}$ and $\mathbb{e}$ only, and the first constraint
reduces to $\mathbb{D}_{1}\mathbb{m}=\mathbb{g}_{2}$; as a consequence, the
updates (12), (15), (16) and (11e)-(11g) are performed taking
$\mathbb{g}_{1}^{k-1}=\mathbb{g}_{1}^{k}=\mathbf{0}$, and update (13) is
discarded. To quantitatively evaluate the accuracy of each algorithm we
compute the 2-norm relative error at the $k$th iteration, $k=1,2,\dots$,
defined as
$\nicefrac{{\|\mathbb{m}^{k}-\mathbb{m}\|_{2}}}{{\|\mathbb{m}\|_{2}}}$.

###  4.1 Compressed sensing and denoising

The compressed sensing theory allows to recover a sparse signal from a small
number of random projections [10]. In this subsection, we first show the
performance of the new automated Tikhonov-TV regularization method in
recovering 1D signals with different features, in the framework of the
compressed sensing theory. We consider four test signals, ranging from very
smooth to very rough: these are shown at the top row of Fig. 1 (dashed red
lines). The observation matrix $\mathbb{G}$ consists of 1024 random column
vectors of length 250, all drawn from a standard Gaussian distribution and
then normalized. Each observation vector is generated as $\mathbb{d=Gm+e}$
where $\mathbb{m}$ is the test signal and $\mathbb{e}$ is some Gaussian white
noise (such that the ‘noise level’
$\|\mathbb{e}\|_{2}/\|\mathbb{G}\mathbb{m}\|_{2}$ is 0.1%). To evaluate the
long-term behavior of Algorithm 1, we stop the method when it has computed 500
iterations. The first row of Fig. 1 shows the reconstructed signals (in blue)
overlaid by the ground truth (dashed, in red), while the nonsmooth and smooth
components of the regularized solution are shown in the second and third rows,
respectively. We can observe that nearly optimal reconstructions and
decompositions are obtained for all of the signals. The model gradients of the
recovered signals are shown in the last row of Fig. 1 (in blue), with the
gradient of the recovered smooth component (dashed red curve) overlaid. In
each frame in the bottom row of Fig. 1, the lower bounds for the anomalous
gradient entries as determined by the z-score (21), (22) are shown by
horizontal dashed lines. We can observe that, in all of the cases, the z-score
determined meaningful bounds, which resulted in an optimal determination of
$\beta$. Fig. 2 shows the evolution of $\beta$ and $\phi(\beta)$ versus
iteration count for all the four signals, with different initial guesses
$\beta_{0}$. We can observe that formula (27) stably converges to the root of
$\phi$, irrespective of the initial value $\beta_{0}$.

Figure 1: Compressed sensing problem. First row: reconstructed signal (in
blue) overlaid by the ground truth (in red). Second and third row: the blocky
component $\mathbb{m}_{1}$ and the smooth component $\mathbb{m}_{2}$,
respectively. Last row: the gradient of the reconstruction (in blue) overlaid
by the gradient of the smooth component $\mathbb{m}_{2}$ (in red); the
horizontal dashed lines show the lower value determined for anomalous
components. Figure 2: Compressed sensing problem. Evolution of the balancing
parameter $\beta$ (top row) and $\phi(\beta)$ (bottom row) versus iteration
count, for different initial values $\beta_{0}$ of $\beta$, corresponding to
the signals shown in Fig. 1.

We then turn to a 2D image denoising example. The image which is to be
denoised, displayed in Fig. 3(a), is obtained by adding some Gaussian white
noise of noise level $\|\mathbb{e}\|_{2}/\|\mathbb{m}\|_{2}=30\%\ $ to a clean
synthetic image which consists of piecewise-smooth regions, displayed in Fig.
3(b). We use the combined Tikhonov-TV method and, to evaluate its
regularization performance, we compare with the TV-only and the Tikhonov-only
methods to denoise the image (as explained above). We run each method for 500
iterations (to evaluate the long-term behavior of Algorithm 1 even after the
stopping criterion based on the stabilization of the solution is satisfied)
and the obtained results are depicted in Fig. 3(c)-(e). As expected, the
Tikhonov-TV method estimates a more accurate image by balancing between
restoration of the sharp edges and of the smooth regions; instead, the TV-only
and Tikhonov-only methods give more emphasis on the restoration of the sharp
edges and smooth regions, respectively, and hence provide a suboptimal
estimate of the noiseless image. The evolution of the discrepancy for all
three methods is depicted in Fig. 4(a). Also, Figs. 4(b)-(c) show the
evolution of $\beta$ and $\phi(\beta)$ versus iteration count for the
Tikhonov-TV method.

Figure 3: Image denoising example. (a) degraded image by a Gaussian noise of
(noise level 30%); (b) noise-free image; (c) estimated image by the proposed
method (error = $1.98\times 10^{-4}$); (d) estimated image by the TV method
(error = $3.18\times 10^{-4}$); (e) estimated image by the Tikhonov method
(error = $4.08\times 10^{-4}$). Figure 4: Image denoising example. Evolution
of (a) discrepancy, (b) $\beta$ and (c) $\phi(\beta)$ versus iteration for the
proposed method. The asterisk highlights the quantities computed at the 135th
iteration, i.e., when the stopping criterion on the stabilization of the
solution is satisfied.

For this test problem, we experimentally study the sensitivity of the new
Algorithm 1 with respect to a couple of algorithmic parameters that should be
given in input, namely: the value of the noise magnitude
$\varepsilon=\|\mathbb{e}\|_{2}^{2}$ and the threshold
$\tau_{\text{\tiny{nrm}}}$ for the normal entries of $\mathbb{g}$ (22).
Concerning the former, we consider cases where the true value of $\varepsilon$
is over- and under-estimated by 5% and 10%. The results of these tests are
displayed in Fig. 5, (a)-(d). In particular, looking at frame (a), it is
evident that, as the iterations proceed, the value of the squared discrepancy
$\|\mathbb{G}\mathbb{m}-\mathbb{d}\|_{2}^{2}$ stabilizes around the inputted
values of $\varepsilon$ (leading to over- and under-fitted data when
$\varepsilon$ is under- and over-estimated, respectively). Looking at frames
(b)-(c) we can see that the update rule (27) for $\beta$ converges to a zero
of $\phi(\beta)$ for all the considered values of $\varepsilon$, although the
values of $\beta$ determined at the end of the iterations differ. The impact
of an inaccurate value of $\varepsilon$ on the quality of the solution is
visible in frame (d): it is clear that, for this test problem, an under-
estimation of $\varepsilon$ (corresponding to an over-fitting) of the data
still leads to results comparable to the case where an accurate value of
$\varepsilon$ is considered; naturally, the higher the error in the estimate
of $\varepsilon$, the lower the quality of the computed solution. In Fig. 5
(e) we illustrate some of the entries of the recovered noise vector
$\mathbb{e}$ at the end of the ADMM iterations. The new ADMM formulation
allows recovery of the random noise vector $\mathbb{e}$ (through subproblem
(9d)), alongside the inverse problem solution $\mathbb{m}$ and its gradient
components $\mathbb{g}_{1}$ and $\mathbb{g}_{2}$. We can clearly see that, for
this test problem, the computed approximation of $\mathbb{e}$ quite carefully
reproduces the behavior of the unknown noise corrupting the original image.

Figure 5: Image denoising example. Sensitivity of Algorithm 1 with respect to
over and under estimations of $\varepsilon$. (a)-(d) evolution of (a)
discrepancy, (b) $\beta$, (c) $\phi(\beta)$, and (d) 2-norm of the estimated
model error. Frame (e) shows approximation of the entries 166 to 346 (at
column 256) of noise vector $\mathbb{e}$ at the 500th (final) ADMM iteration.
In this frame, the black line shows the true noise.

Concerning the sensitivity of Algorithm 1 with respect to the value of
$\tau_{\text{\tiny{nrm}}}$ we consider the default value
$\tau_{\text{\tiny{nrm}}}=2.5$ and $\tau_{\text{\tiny{nrm}}}=2$,
$\tau_{\text{\tiny{nrm}}}=3$; the behavior of the solver for these tests are
displayed in Figure 6. Looking at frame (a) we can clearly see that a higher
value of $\tau_{\text{\tiny{nrm}}}$ results in a lower value of the computed
$\beta$ at the end of the iterations of Algorithm 1: this is expected, as a
higher $\tau_{\text{\tiny{nrm}}}$ implies more entries of $\mathbb{g}$ to be
regarded as ‘normal’ (see equation (22)), which is achieved by imposing less
penalisation on the smoothness-enforcing term in (7). Similarly, a lower
$\tau_{\text{\tiny{nrm}}}$ results in a higher value of $\beta$, i.e., less
entries of $\mathbb{g}$ are regarded as ‘normal’ and the smoothness-enforcing
term in (7) is penalised more. Looking at frame (b) we can see that such
variations in $\tau_{\text{\tiny{nrm}}}$ do not impact the quality of the
solution computed by Algorithm 1. Even if not reported, a similar behavior is
observed testing a wider range of $\tau_{\text{\tiny{nrm}}}$ values.

Figure 6: Image denoising example. Evolution of (a) $\beta$ and (b) 2-norm
relative error versus iteration for the proposed method for different values
of the threshold $\tau_{\tiny{\text{nrm}}}$ in (22).

Finally, we experimentally assess the convergence properties of the new ADMM
scheme (9a)-(9f) applied with iteration-dependent choice of the balancing
parameter $\beta$ according to (27), as well as with a fixed value
$\beta=\beta_{500}$, i.e., picking the value the fixed point iteration (27)
converged to after 500 iterations. Fig. 7 shows the residual and error values
versus iteration number: we can clearly see that, although the two instances
of ADMM display some discrepancies for approximately the first 100 iterations,
they eventually converge to the same value. Even if not reported, this
behavior is observed in all the considered test problems.

Figure 7: Image denoising example. Evolution of (a) discrepancy and (b) 2-norm
of the estimated model error for Algorithm 1 with iteration-dependent choice
of $\beta$ and for the ADMM method (9a)-(9f) with fixed $\beta=\beta_{500}$.

###  4.2 Inverse problems in Geophysics

We apply the new automated Tikhonov-TV regularization method for subsurface
interval-velocity estimation from root-mean-square (RMS) velocities. For a
horizontally layered earth model, the RMS velocity $V(t)$ is a continuous
function of time defined as

$V(t)=\sqrt{\frac{1}{t}\int_{0}^{t}[v(\tau)]^{2}d\tau},$ (28)

where $v(t)$ is the instantaneous/interval velocity function and $t$ is the
wave propagation time; see [9]. In practice, the velocity analysis of common-
depth-point (CDP) gathers gives an estimate of the RMS velocity; when both the
RMS velocity $V(t)$ and the interval velocity $v(t)$ are needed, the
estimation of the latter from the RMS velocity is a severely ill-conditioned
problem, thus proper regularization is required to stabilize the solution.
Upon discretization, equation (28) reads

$V_{i}=\sqrt{\frac{1}{i}\sum_{j=1}^{i}v^{2}_{j}},$ (29)

where $i=1,...,N$ is the sample number and $N$ is the number of earth layers.
Squaring both sides of equation (29) results in a system of linear equations
of the form $\mathbb{d}=\mathbb{G}\mathbb{m}$, where $d_{i}=iV_{i}^{2}$,
$m_{j}=v_{j}^{2}$ and the discrete forward operator $\mathbb{G}$ is a causal
integration matrix, i.e., a lower triangular matrix of ones

$\mathbb{G}=\begin{pmatrix}1&&&&\\\ 1&1&&&\\\ 1&1&1&&\\\
\vdots&\vdots&\ddots&\ddots&\\\
1&1&\cdots&1&1\end{pmatrix}\in\mathbb{R}^{N\times N}.$ (30)

We simulate a 1D RMS velocity vector from a velocity log of the 2004 BP model
[5]. The interval and RMS velocities are shown in Figs. 8(a) and (b),
respectively. In order to make the simulation more realistic, only 25% of the
RMS velocity elements where used as input to the inversion, mimicking the real
situations where one only picks RMS velocities at the positions of strong
reflections. This leads to a compressed sensing problem for the Dix inversion
[13]. For this example, the forward operator is $\Phi\mathbb{G}$, where $\Phi$
is a 382 by 1911 sampling matrix (made up of 382 rows of an identity matrix of
size 1911, picked at the locations of the velocity model with strong seismic
response/reflectivity) and $\mathbb{G}$ is defined as in (30) and has size
$1911\times 1911$. The estimated interval velocity obtained applying 100
iterations of Algorithm 1 is shown in Fig. 8(c); we can clearly observe that
the interval velocities are estimated accurately.

We then consider a 2D Dix inversion example. The full RMS velocity field is
shown in Fig. 9(a). We used 30% of the traces (vertical lines, picked at
random) (shown in Fig. 9(b)) as the input to the inversion. In this case, the
size of the RMS velocity field is 478 $\times$ 1349 pixels and the forward
operator is of the form $\Phi\otimes\mathbb{G}$, where $\Phi$ is a sampling
matrix of size $404\times 1349$ (made up of 404 rows of an identity matrix of
size 1349, picked at random) and $\mathbb{G}$ is defined as in (30) and has
size $478\times 478$. The estimated full interval velocity field obtained
applying 100 iterations of Algorithm 1 is shown in Fig. 9(c); also for this
test problem we can clearly see that the interval velocity model is
reconstructed accurately. Fig. 10 shows the evolution of $\beta$ and
$\phi(\beta)$ versus iteration for the considered 1D and 2D Dix inversion
examples.

Figure 8: One-dimensional Dix velocity inversion. (a) True interval velocity.
(b) Decimated noisy RMS velocity. (c) Estimated interval velocity. Figure 9:
Two-dimensional Dix velocity inversion. (a) RMS velocity field corresponding
to the 2004 BP model. (b) Input RMS field from velocity scan using 30% of CDP
gathers selected at random. (c) Interval velocity field obtained by inversion
using the proposed method. Figure 10: Dix velocity inversion. Left frame:
values of the balancing parameter $\beta$ versus iteration number. Right
frame: values of $\phi(\beta)$ versus iteration number. Solid curves
correspond to the 1D case, while dashed curves correspond to the 2D case. The
asterisk highlights iteration for which the stopping criterion on the
stabilization of the solution is satisfied.

Finally, still using two popular datasets in the seismic community, we assess
the performance of our algorithm when used for image decomposition. Fig.
11(a)-(b), top row, show the 2004 BP velocity model [5] and the 2007 BP
velocity model [23], respectively; the velocity variation is a piecewise
smooth function of space for both these models. We apply 500 iterations of the
Tikhonov-TV regularization method to decompose each velocity model into smooth
and nonsmooth components: the results are still displayed in Fig. 11, on the
middle and bottom rows, respectively. We can clearly observe that the two
components $\mathbb{m}_{1}$ and $\mathbb{m}_{2}$ are optimally separated and
allow to recover complementary parts of the original velocity model. Fig. 12
shows the evolution of the values of $\beta$ and $\phi(\beta)$ versus
iteration for these two test problems. We can observe that, in both cases, the
simple update scheme in (27) stably converges to a root of $\phi(\beta)$,
leading to an optimal separation between the smooth and nonsmooth components
of the model $\mathbb{m}$.

Figure 11: Model decomposition into smooth and nonsmooth components. Column
(a): the 2004 BP velocity model. Column (b): the 2007 BP velocity model. Top
row: the original models $\mathbb{m}$. Middle row: the nonsmooth component
$\mathbb{m}_{1}$. Bottom row: the smooth component $\mathbb{m}_{2}$. Figure
12: Model decomposition into smooth and nonsmooth components. Left frame:
values of the balancing parameter $\beta$ versus iteration number. Right
frame: values of $\phi(\beta)$ versus iteration number. The asterisk
highlights iteration for which the stopping criterion on the stabilization of
the solution is satisfied.

###  4.3 X-ray tomography

We first consider a parallel X-ray computed tomography test problem. The
unknown $\mathbb{m}$ to be reconstructed is the Shepp-Logan phantom of size
$320\times 320$ pixels. The original phantom is shown on the leftmost frame of
Fig. 13. We consider a full projection setting with 453 parallel rays for each
of the 90 equispaced angles in range $[-90^{\circ},+90^{\circ}]$. This results
in a discrete forward operator $\bf G$ of size $40770\times 102400$. Gaussian
white noise of level $1\%$ is added to the data. During the inversion stage,
at the $k$th ADMM iteration, the subproblem (12) for $\mathbb{m}$ is solved by
using the CG method. The CG iterations are stopped when the relative residual
tolerance $10^{-7}$ is reached, or when the maximum number of 100 iterations
is reached; the estimated solution at the $(k-1)$th ADMM iteration is used as
an initial guess for the CG at current iteration. This resulted in a dynamic
number of inner CG iterations, as shown in Fig. 15. The behavior of other
relevant quantities is displayed in Fig. 14. In particular, we can observe
that for this test problem, whose solution solely contains piecewise constant
features, there are almost no differences in the performance of Tikhonov-TV
and TV regularization methods.

Figure 13: X-ray tomography test problem. From left to right: true image (320
$\times$ 320 pixels), estimation by the proposed balanced Tikhonov-TV method,
estimation by the TV method, and estimation by the Tikhonov method. The
numbers in parenthesis report the relative error of each estimate.  Figure 14:
X-ray tomography test problem. Evolution of the: (a) relative errors, (b)
squared discrepancy versus iteration for the methods considered in Fig. 13.
Evolution of: (c) $\beta$, (d) $\phi(\beta)$ for the Tikhonov-TV method. The
asterisk highlights iteration for which the stopping criterion on the
stabilization of the solution is satisfied.  Figure 15: X-ray tomography test
problem. The number of CG iterations performed at each ADMM iteration.

We then take a phantom containing both piecewise constant and smooth features,
and we resize it to be $128\times 128$ pixels; the resulting phantom is shown
on the leftmost frame of Fig. 16. We consider a setting where projections
along parallel X-rays can be performed only along limited angles:
specifically, we consider 181 parallel rays for each of the 85 equispaced
angles in $[-42^{\circ},+42^{\circ}]$. This results in a discrete forward
operator $\mathbb{G}$ of size $15385\times 16384$. Gaussian white noise of
level $0.1\%$ is added to the data. For this particular test problem, ADMM
applied to solve the Tikhonov-TV and TV-only test problem display a very slow
convergence (this is visible in Fig. 17 (b), where the only method whose
squared residual approximately equals $\varepsilon$ within 600 iterations is
Tikhonov regularization). Despite this, the reconstruction computed by the
Tikhonov-TV method achieves the lowest relative error among the considered
methods, with a good resolution of the background, the piecewise constant and
the smooth features: this is clearly visible in Fig. 16. Fig. 17 reports the
progress of other relevant quantities versus the iteration count.

exact | Tikhonov-TV (0.2470) | TV (0.2656) | Tikhonov (0.2620)
---|---|---|---
|  |  |
Figure 16: X-ray tomography test problem, with limited angles. Exact phantom
along with the reconstructions achieved after 600 iterations of each methods.
The numbers between brackets are the 2-norm relative errors associated to each
reconstruction. To better highlight differences, the square root of the
modulus of each pixel is displayed. Figure 17: X-ray tomography test problem,
with limited angles. Frame (a): 2-norm relative errors versus iteration number
for different methods. Frame (b): squared 2-norm discrepancies versus
iteration number for different methods. Frame (c): values of $\beta$ versus
iteration number for Tikhonov-TV. Frame (d): values of $\phi(\beta)$ versus
iteration number for Tikhonov-TV.

##  5 Conclusions

In this paper, we proposed a method for selecting the optimal balancing
parameter in Tikhonov-TV regularization for the solution of discrete ill-posed
inverse problems with piecewise smooth solutions. Namely, the solution of the
inverse problem is split into the sum of a smooth and a piecewise constant
component, which are separately regularized by a Tikhonov and a TV term that
must be balanced. We have used robust statistical methods for determining the
optimal balance of these terms, motivated by the fact that the gradient
entries of the piecewise constant component at jump locations can be
considered as anomalies/outliers in comparison to the other gradient entries
constituting the smooth background. This led to the characterization of the
best balancing parameter as a root of a scalar function. Finally, an extremely
simple update scheme was proposed for determining the balancing parameter,
which can be naturally coupled with ADMM to solve the Tikhonov-TV regularized
problem. Extensive numerical experiments on different inverse problems
demonstrate that high quality reconstructions can be obtained by applying the
proposed algorithm, and validate the robustness of the proposed selection
strategy of balancing parameter.

##  Appendix A Proof that the function $\psi(\beta)$ in (27) is not expansive.

Let $0<\beta_{1}<\beta_{2}$. Let us assume that
$\|\text{nrm}(\mathbb{g}(\beta_{1}))\|_{\infty}\leq\|\text{nrm}(\mathbb{g}(\beta_{2}))\|_{\infty}$.
This is a meaningful assumption since
$\text{nrm}(\mathbb{g}(\beta_{1}))\subset\text{nrm}(\mathbb{g}(\beta_{2}))$
because $\mathbb{g}(\beta_{1})$ is expected to be less smooth than
$\mathbb{g}(\beta_{2})$ (because the smooth component $\mathbb{g}_{2}$ of
$\mathbb{g}$ is weighted less in the objective function (7)). We know that
$\|\mathbb{g}_{2}(\beta_{1})\|_{2}>\|\mathbb{g}_{2}(\beta_{2})\|_{2}$ (still
as a consequence of the weighting in the objective function (7)). Since the
2-norm and $\infty$-norm are equivalent on $\mathbb{R}^{N}$ (namely, for all
$\mathbb{x}\in\mathbb{R}^{N}$,
$\|\mathbb{x}\|_{\infty}\leq\|\mathbb{x}\|_{2}\leq
N^{1/2}\|\mathbb{x}\|_{\infty}$), the 2-norm inequality above implies the
$\infty$-norm inequality $\|\mathbb{g}_{2}(\beta_{1})\|_{\infty}\geq
N^{-1/2}\|\mathbb{g}_{2}(\beta_{2})\|_{\infty}$. To keep the notations light,
in the following we let

$n_{1}:=\|\text{nrm}(\mathbb{g}(\beta_{1}))\|_{\infty},\;n_{2}:=\|\text{nrm}(\mathbb{g}(\beta_{2}))\|_{\infty},\;d_{1}:=\|\mathbb{g}_{2}(\beta_{1})\|_{\infty},\;d_{2}:=\|\mathbb{g}_{2}(\beta_{2})\|_{\infty}.$

As a consequence of the inequalities above, which can be rewritten as

$n_{1}\leq n_{2}\quad\mbox{and}\quad d_{1}\geq
N^{-1/2}d_{2}\quad\left(\mbox{so
that}\quad\frac{n_{1}}{d_{1}}<\frac{n_{2}}{N^{-1/2}d_{2}}\right)\,,$ (31)

the term appearing in the right-hand side parenthesis in (27) can be rewritten
and bounded as follows

$\displaystyle\left(4+\frac{n_{1}}{d_{1}}\right)^{-1}-1$ $\displaystyle=$
$\displaystyle
N^{-1/2}\left(\left(4N^{-1/2}+\frac{n_{1}}{N^{1/2}d_{1}}\right)^{-1}-N^{1/2}\right)$
$\displaystyle>$ $\displaystyle
N^{-1/2}\left(\left(4N^{-1/2}+\frac{n_{2}}{d_{2}}\right)^{-1}-N^{1/2}\right)$
$\displaystyle\geq$ $\displaystyle
N^{-1/2}\left(\left(4+\frac{n_{2}}{d_{2}}\right)^{-1}-N^{1/2}\right)\,,$

where, in the first inequality, we have exploited (31). Therefore, the modulus
of the difference $\psi(\beta_{1})-\psi(\beta_{2})$ can be rewritten and
bounded as

$\displaystyle\left|\psi(\beta_{1})-\psi(\beta_{2})\right|$ $\displaystyle=$
$\displaystyle\left|\left(\left(4+\frac{n_{1}}{d_{1}}\right)^{-1}-1\right)\beta_{1}-\left(\left(4+\frac{n_{2}}{d_{2}}\right)^{-1}-1\right)\beta_{2}\right|$
$\displaystyle<$
$\displaystyle\left|\left(\left(4+\frac{n_{1}}{d_{1}}\right)^{-1}-1\right)\beta_{1}-\left(N^{-1/2}\left(4+\frac{n_{2}}{d_{2}}\right)^{-1}-1\right)\beta_{2}\right|$
$\displaystyle\leq$
$\displaystyle\left|\left(\left(4+\frac{n_{1}}{d_{1}}\right)^{-1}-1\right)\beta_{2}-\left(\left(4+\frac{n_{1}}{d_{1}}\right)^{-1}-1\right)\beta_{1}\right|$
$\displaystyle=$
$\displaystyle\left|\left(\left(4+\frac{n_{1}}{d_{1}}\right)^{-1}-1\right)\right||\beta_{1}-\beta_{2}|\,.$

Since the first factor on the left-hand side is positive and strictly less
than 1, we can conclude that $\psi(\beta)$ is not expansive.

## References

  * [1]
  * Aghamiry et al. [2019] H. Aghamiry, A. Gholami, S. Operto,  Compound regularization of full-waveform inversion for imaging piecewise media,  IEEE Transactions on Geoscience and Remote Sensing (2020) 58 1192-1204.
  * Alrajawi et al. [2017] M. Alrajawi, H. Siahkoohi, A. Gholami,  Inversion of seismic arrival times with erratic noise using robust Tikhonov-TV regularization,  Geophysical Journal International 211 (2017) 831–842.
  * Aster et al. [2004] R. C. Aster, B. Borchers, C. H. Thurber, Parameter Estimation and Inverse Problems, Academic Press, 2004.
  * Billette and Brandsberg-Dahl [2004] F. J. Billette, S. Brandsberg-Dahl,  The 2004 BP velocity benchmark,  in: Extended Abstracts, 67th Annual EAGE Conference & Exhibition, Madrid, Spain, 2004, p. B035.
  * Boyd et al. [2010] S. Boyd, N. Parikh, E. Chu, B. Peleato, J. Eckstein,  Distributed optimization and statistical learning via the alternating direction of multipliers,  Foundations and trends in machine learning 3 (2010) 1–122.
  * [7] P. Borwein and T. Erdélyi, Polynomials and Polynomial Inequalities, Springer, New York, 1995.
  * Chen and Teboulle [1994] G. Chen, M. Teboulle,  A proximal-based decomposition method for convex minimization problems,  Mathematical Programming 64 (1994) 81–101.
  * Dix [1955] C. H. Dix,  Seismic velocities from surface measurements,  Geophysics 20 (1955) 68–86.
  * Donoho et al. [2006] D. L. Donoho, et al.,  Compressed sensing,  IEEE Transactions on information theory 52 (2006) 1289–1306.
  * Gabay and Mercier [1976] D. Gabay, B. Mercier,  A dual algorithm for the solution of nonlinear variational problems via finite element approximation,  Computers & mathematics with applications 2 (1976) 17–40.
  * Gholami and Hosseini [2013] A. Gholami, S. M. Hosseini,  A balanced combination of Tikhonov and total variation regularizations for reconstruction of piecewise-smooth signals,  Signal Processing 93 (2013) 1945–1960.
  * Gholami and Naeini [2019] A. Gholami, E. Z. Naeini,  3D Dix inversion using bound-constrained TV regularization,  Geophysics 84 (2019) 1–43.
  * Hansen [2010] P. C. Hansen, Discrete Inverse Problems: Insight and Algorithms, SIAM, 2010.
  * Huska et al. [2019] M. Huska, A. Lanza, S. Morigi, I. Selesnick,  A convex-nonconvex variational method for the additive decomposition of functions on surfaces,  Inverse Problems 35 (2019) 124008.
  * Iglewicz and Hoaglin [1993] B. Iglewicz, D. C. Hoaglin, How to detect and handle outliers, Quality Press, 1993.
  * Kalmoun [2020] E. M. Kalmoun,  A new Tikhonov-TV regularization for optical flow computation,  in: Eleventh International Conference on Graphics and Image Processing (ICGIP 2019), volume 11373, International Society for Optics and Photonics, 2020, p. 113732B.
  * Nocedal and Wright [2006] J. Nocedal, S. Wright, Numerical optimization, Springer Science & Business Media, 2006.
  * Rousseeuw and Croux [1993] P. J. Rousseeuw, C. Croux,  Alternatives to the median absolute deviation,  Journal of the American Statistical association 88 (1993) 1273–1283.
  * Rousseeuw and Hubert [2018] P. J. Rousseeuw, M. Hubert,  Anomaly detection by robust statistics,  Wiley Interdisciplinary Reviews: Data Mining and Knowledge Discovery 8 (2018) e1236.
  * [21] E. Ryu, and S. Boyd,  A primer on monotone operator methods,  Appl. Comput. Math 15 1 (2016) 3–43.
  * Selesnick et al. [2014] I. W. Selesnick, H. L. Graber, D. S. Pfeil, R. L. Barbour,  Simultaneous low-pass filtering and total variation denoising,  IEEE Transactions on Signal Processing 62 (2014) 1109–1124.
  * Shah [2007] H. Shah,  The 2007 BP anisotropic velocity-analysis benchmark,  in: Expanded Abstracts, EAGE workshop, 2007.
  * Strang [1999] G. Strang,  The Discrete Cosine Transform,  SIAM Review 41 1 (1999) 135–147.
  * Tarantola [2005] A. Tarantola, Inverse Problem Theory and Methods for Model Parameter Estimation, SIAM, Philadelphia, 2005.
  * Vasin [2018] V. V. Vasin, V. V. Belyaev,  Approximation of solution components for ill-posed problems by the Tikhonov method with total variation,  Doklady Mathematics 97 (2018) 266-270.
  * [27] J. Wang, L. Zhao,  Nonconvex generalization of Alternating Direction Method of Multipliers for nonlinear equality constrained problems,  Results in Control and Optimization 2 (2021) 266-270.