text:
\begin{document}

\title{A Gradient-thresholding Algorithm for Sparse Regularization}

\author{Abinash Nayak}

\address{Visiting Assistant Professor, Department of Mathematics, University of Alabama at Birmingham, University Hall, Room 4005, 1402 10th Avenue South, Birmingham AL 35294-1241, (p) 205.934.2154, (f) 205.934.9025}

\email{nash101@uab.edu; avinashnike01@gmail.com}

\subjclass{Primary 65R30, 65R32; Secondary 65R20, 65K10}

\date{\today}

\keywords{Inverse problems, Iterative regularization, Ill-posed problems, Tikhonov regularization, L1 regularization, LASSO, sparsity, ISTA, FISTA.}

\begin{abstract}
Inverse problems arise in a wide spectrum of applications in fields ranging from engineering to scientific computation. Connected with the rise of interest in inverse problems is the development and analysis of regularization methods, such as Tikhonov-type regularization methods or iterative regularization methods, which are a  necessity in most of the inverse problems. In the last few decades, regularization methods motivating sparsity has been the focus of research, due to the high dimensionalty of the real-life data, and $\mathcal{L}^1$-regularization methods (such as LASSO or FISTA) has been in its center (due to their computational simplicity). In this paper we propose a new (semi-) iterative regularization method which is not only simpler than the mentioned algorithms but also yields better results, in terms of accuracy and sparsity of the recovered solution. Furthermore, we also present a very effective and practical stopping criterion to choose an appropriate regularization parameter (here, it's iteration index) so as to recover a regularized (sparse) solution. To illustrate the computational efficiency of this algorithm we apply it to numerically solve the image deblurring problem and compare our results with certain standard regularization methods, like total variation, FISTA, LSQR etc.
\end{abstract}
\maketitle

\section{\textbf{Introduction}}
\subsection{Inverse Problems and Regularization:}
An inverse problem in general is a problem where the effect (output) is known but the source (input) is not, in contrast to a direct problem where we deduce the effect from the source. Mathematically, an inverse problem is often expressed as the problem of finding/estimating a (\textit{source function}) $\varphi$, given an (\textit{effect data}) $g$, which satisfies the following operator equation:
\begin{equation}\label{T Gen.}
T\varphi = g,
\end{equation}
$T$ being some operator describing the underlying physical process, where the domain and range of the operator $T$ varies depending on the problem (usually, it's a subset of $\mathcal{L}^2(\Omega)$, for $\Omega \subset \mathbb{R}^n$). Typically, the solution of \eqref{T Gen.} is approximated by the solution of the least-square problem, i.e., minimizer of the following functional
\begin{equation}\label{LS solution}
    F(\psi) = ||T\psi - g||_2^2.
\end{equation}

However, inverse problems are usually ill-posed, in the sense of violating Hadamard's third condition: ``Continuous dependence of the data", i.e., even for a slightly perturbed data $g_\delta$, such that $||{g - g_\delta}|| \leq \delta$ (usually small), the inverse recovery becomes unstable, $||{\varphi - \varphi_\delta}|| >> \delta$ (very large); typically due to the unboundedness of the (pseudo-) inverse operator $T^\dagger$ and the noise present in the data. To counter such instabilities or the ill-posedness of inverse problems, regularization methods have to be employed. In the last few decades, several regularization methods have been established for linear as well as nonlinear inverse problems. Broadly, there exist two kinds of regularization approaches:

\subsection{Tikhonov-type regularization:}
Tikhonov-type regularization methods are probably the most well known regularization techniques for solving linear as well as nonlinear inverse problems (see \cite{Engl+Hanke+Neubauer,Bakushinsky+Goncharsky, Groetsch, Baumeister,Morozov_a}), where (instead of minimizing the simple least-square problem \eqref{LS solution}) one recovers a \textit{``regularized solution"} by minimizing a (constrained or) penalized functional
\begin{equation}\label{Tikhonov}
F(\psi;\lambda,L,\psi_0,p,q) =  ||{T\psi - g_\delta}||_p^p + \lambda ||L(\psi - \psi_0)||_q^q,
\end{equation}

(for some $p$ and $q$) where $\lambda > 0$ is the called the regularization parameter, $||g - g_\delta|| \leq \delta$ is the error norm, $\psi_0$ is an initial guess and $L$ is a regularization operator\footnote{the space for $\psi$ and the norms involved are defined appropriately.}, with the null spaces of $T$ and $L$ intersecting trivially. In this approach, the ill-posed problem \eqref{T Gen.} is first converted to a family of well-posed problems (depending on a regularization parameter $\lambda$) and once an appropriate choice of the parameter value (say $\lambda_0$) is estimated, a regularized solution $\varphi^\delta$ of \eqref{T Gen.} is recovered by \textbf{completely} minimizing the associated functional, i.e.,
\begin{equation}
    \varphi^\delta_{\lambda_0} = \argmin_{\psi} F(\psi;\lambda_0,L,\psi_0,p,q).
\end{equation}
The choice of an appropriate parameter value $\lambda_0$ is crucial in any regularization method, as it balances between the data fitting term $||{T\psi - g_\delta}||_p^p$ and the regularization term $||L(\psi - \psi_0)||_q^q$. If $\lambda$ is too small then minimizing \eqref{Tikhonov}
\textit{over-fits} the noisy data (thus, leading to a \textit{noisy recovery} or, statistically speaking, a \textit{high-variance recovery}), where as if $\lambda$ is too big then it \textit{under-fits} the data (thus, leading to an \textit{over-smooth recovery}, or statistically, a \textit{high-bias recovery}). Typically, one recovers a bunch of $\lambda$-dependent regularized solutions $\{ \varphi^\delta_{\lambda}: \lambda \in [\lambda_{min}, \lambda_{max}]\}$ and then choose an appropriate regularized solution $(\varphi^\delta_{\lambda_0})$ depending on some parameter choice criterion, say DP (discrepancy principle), GCV (generalized cross-validation), L-curve, monotone error rule etc., see \cite{Morozov1966,Vainikko1982, Gfrerer1987, Hansen1992, Lawson_Hanson,Tautenhahn_Hamarik,Bauer_Hohage, Mathe}.

\subsection{(Semi-) Iterative regularization:}\label{semi-iterative methods}
Where as, in an iterative regularization method one recovers a regularized solution $\varphi^\delta$ of \eqref{T Gen.} by stopping the minimization process of the simple least-square problem \eqref{LS solution} (i.e., no additional penalty term) at an appropriate (early) instance, i.e., starting from an initial guess $\psi_0^\delta$, one updates the $m^{th}$ recovered solution ($\psi_m^\delta$) in a direction ($\xi_m^\delta$) such that the $(m+1)^{th}$ recovered solution ($\psi_{m+1}^\delta$) is better than the previous, meaning $||\psi_{m+1}^\delta-\varphi|| < ||\psi_m^\delta - \varphi||$, and stops the iteration after certain number of steps (i.e., $m  \leq m_0(\delta) < \infty$, for some appropriate $m_0(\delta)$, usually depending on the noise level $\delta$). Mathematically, it can be expressed as
\begin{equation}\label{steepest descent}
    \psi_{m+1}^\delta = \psi_m^\delta + \tau_m \xi_m^\delta, \mbox{ for } m \leq m_0(\delta) < \infty,
\end{equation}{}
where $\tau_m$ is the m-th step-length and the descent direction $\xi_m^\delta$ is usually the (negative) gradient direction ($\xi_m^\delta = -T^*(T\psi_m^\delta - g_\delta)$). When the step-size is fixed ($0<\tau_m = \tau < \frac{2}{||T^*T||}$) for all $m$, it's known as Landweber iterations (also known as Richardson iterations) and has been intensively investigated in the literature (see, \cite{Landweber, Hanke_Neubauer_Scherzer, Engl+Hanke+Neubauer, Bakushinsky+Goncharsky, Hanke1991, Schock1986} and references therein). The main drawback of Landweber iterations is its slow performance, i.e., it takes a large number of iterations to obtain the optimal convergence rates. To circumvent this drawback many extensions (known as polynomial accelerated Landweber methods) have been proposed and studied in the framework of regularization (for an overview, see \cite{Hanke1991, Schock1986}). Such iterative methods are typically known as \textit{semi-iterative methods}, since when dealing with a noisy data ($g_\delta$) one encounters a semi-convergent nature of the recovery process, see \cite{Engl+Hanke+Neubauer}. That is, in the initial stage $||\psi_m^\delta - \varphi|| \rightarrow 0$ for $m \leq m_0(\delta) < \infty$, but upon further iterations $||\psi_m^\delta - \varphi|| \rightarrow \infty$ (or, $||\psi_m^\delta|| \rightarrow \infty$), as $m \rightarrow \infty$.
Hence, to achieve regularization one needs to terminate the descent process at an appropriate early stage, i.e., here the iteration index serves as a regularization parameter. If the iterations are stopped too early ($m$ too small) then that leads to an \textit{over-smooth or high-bias recovery} (i.e, \textit{under-fitting} the data), and if $m$ is large then it leads to a \textit{noisy or high-variance recovery} (i.e., \textit{over-fitting} the noisy data). Typically, the iterations are terminated with the aid of certain stopping rule, such as the \textit{discrepancy principle} (see \cite{Morozov_a}). The advantage of these extended semi-iterative methods over the simple Landweber iteration is that, while Landweber iteration \eqref{steepest descent} uses only the last iterate $\psi_m^\delta$ to construct the new approximate $\psi_{m+1}^\delta$, in a semi-iterative method one make use of the last few iterates (if not all),
\begin{gather}\label{Krylov-space method}
    \psi_{m+1}^\delta = \sum_{i=0}^m \mu_{m,i}\psi_i^\delta + \tau_m \xi_m^\delta,  \\
    \sum_{i=0}^m \mu_{m,i} = 1, \quad \tau_m \neq 0, \notag
\end{gather}{}
where $\xi_m^\delta$ is the descent direction (usually $\xi_m^\delta = -T^*(T\psi_m^\delta - g_\delta)$). Note that, $\psi_{m+1}^\delta - \psi_0^\delta$ belongs to the \textit{Krylov-subspace} $\mathcal{K}_{m+1}(T^*T,\xi_0^\delta)$, which is defined as
\begin{equation}\label{Krylov-subspace}
    \mathcal{K}_{m+1}(T^*T,\xi_0^\delta) := \mbox{span}\{ \xi_0^\delta, (T^*T)\xi_0^\delta, \cdots , (T^*T)^m \xi_0^\delta \}.
\end{equation}{}
Hence, these methods are also called as \textit{Krylov-subspace} methods. One can easily see that choosing $\mu_{m,i} = \delta_{mi}$, where $\delta_{mm}=1$ and $\delta_{mi}=0$, for $m\neq i$, and $\tau_m=\tau$ in \eqref{Krylov-space method} would give back Landweber iterations. It can be further shown (see \cite{Engl+Hanke+Neubauer, Hanke1991}) that there exist sequences $\mu_m$ and $\tau_m$ such that $\psi_{m+1}^\delta$ in \eqref{Krylov-space method} can be reduced to
\begin{equation}
    \psi_{m+1}^\delta = \psi_m^\delta + \mu_m(\psi_m^\delta - \psi_{m-1}^\delta) + \tau_m \xi_m^\delta, \quad m \geq 1.
\end{equation}{}
A specific example of such methods are known as the $\nu$-methods (see \cite{Brakhage1987}), which are defined by $\mu_1=0$, $\tau_1 = \frac{4\nu+2}{4\nu +1}$, and for $m>1$
\begin{align}
    \mu_m &= \frac{(m-1)(2m-3)(2m+2\nu-1)}{(m+2\nu-1)(2m+4\nu-1)(2m+2\nu-3)},\\
    \tau_m &= 4\frac{(2m+2\nu-1)(m+\nu-1)}{(m+2\nu-1)(2m+2\nu-1)}.
\end{align}{}
Even further acceleration is possible by adapting a conjugate-gradient type method (see \cite{Hanke_2017}), where $\mu_{m,i}$ in \eqref{Krylov-space method} depend on data $g_\delta$ (making it a non-linear method).

\subsection{Sparsity Regularization:}
In the last few decades a particular subclass of the regularization methods has drawn a lot of attention, it's called \textit{sparsity regularization} or \textit{regularization with sparsity}, where the support\footnote{the support of a function $\psi$ is defined as $supp_\psi = \{ x \in Domain(\psi): \psi(x) \neq 0\}$} of the solution of \eqref{T Gen.} is significantly smaller than its domain, i.e., $\mu(\{ x: \varphi(x) \neq 0 \} \subset \Omega_\varphi) << \mu(\Omega_\varphi)$, where $\Omega_\varphi$ is the domain of $\varphi$ and $\mu$ is the corresponding measure of the space. In practical situations one only has finitely many observations or measured data and hence, the operator equation \eqref{T Gen.} can be represented in its discretized version (when $T$ is linear, or being approximated by a linear operator) as
\begin{equation}\label{Ax=b}
    A\hat{x} = b,
\end{equation}
where $A \in \mathbb{R}^{m\times n}$ (is usually ill-conditioned), $\hat{x} \in \mathbb{R}^{n}$ and $b \in \mathbb{R}^{m}$. Again, we don't have the exact measurement of the output ($b$), but rather a noisy one ($b_\delta$) such that $||b - b_\delta|| \leq \delta$. From now on (until otherwise stated) we assume $T$ be a linear bounded (injective) operator and, use the symbolic system $(T, \varphi, g)$ to denote the original problem (i.e., $T$ is an operator acting on infinite dimensional spaces, where the functions $\varphi$, $\psi$ and $g$ belongs to) and the system $(A,\hat{x},b)$ to denote its discretized version (i.e., $A$ is $m \times n$ matrix and, $\hat{x}$, $x$ and $b$ are $n$ and $m$ dimensional vectors, respectively). Here, a vector $x \in \mathbb{R}^{n}$ is said to be \textit{sparse} if the number of non-zero elements in $x$ is significantly less than $n$, i.e., $||x||_0 = \sum_i \{ 1 \leq i \leq n: x_i \neq 0, \mbox{ where } x = (x_i)_{i=1}^n \} << n$, this is analogous to the above definition of the support for a function when the function is replaced by a vector and the continuous measure $\mu$ is replaced by a counting measure. Therefore, ideally one would like to minimize the following functional (for an appropriate $\lambda_0$)
\begin{equation}\label{l0-constraint}
    F(x, \lambda_0) = ||Ax - b||_2^2 + \lambda_0 ||x||_0,
\end{equation}
which is known as the \textit{best subset selection} \cite{Beale_Kendall_Mann, Hocking_Leslie}, but due to the non-convexity of $||.||_0$ norm, it's a computational nightmare (NP-hard problem). One can circumvent this problem by having a convex-relaxation of the $||.||_0$-norm, i.e., considering the $||.||_1$ norm (which is convex) instead of the $||.||_0$ norm, where $||x||_1 = \sum_{i=1}^n |x_i|$ for $x = (x_i)_{i=1}^n$. That is, instead of minimizing \eqref{l0-constraint} one minimizes the following functional (for an appropriate $\lambda_0$)
\begin{equation}\label{l1-constraint}
    F(x,\lambda_0) = ||Ax - b||_2^2 + \lambda_0 ||x||_1,
\end{equation}
it's known as \textit{$\mathcal{L}^1$-regularization}, or more popularly \textit{LASSO (least absolute shrinkage and selection operator)} which was proposed by Tibshirani in 1996 \cite{Tibshirani1996}, or \textit{Basis Pursuit Denoising} by Chen et al. in 1998 \cite{Chen_Donoho_Saunders_1998}. $\mathcal{L}^1$-regularization gained its popularity because of its ability to perform variables/features/model selection and its significance in signal and image processing, references on earlier works promoting the use and application of $\mathcal{L}^1$-regularization (together with its relevance to other research areas, such as compressed sensing etc.) can be found in \cite{Figueiredo_Nowak2007,Figueiredo_Nowak2003, Daubechies_Defrise_DeMol2004} and references therein. Note that, the functional defined in \eqref{l1-constraint} (though convex) is not differentiable, due to the non-differentiability of the $||.||_1$-norm, and hence, simple gradient descent algorithm wouldn't work for minimizing \eqref{l1-constraint}, one has to use sub-gradient methods. Also, due to the convexity of the functional \eqref{l1-constraint}, it can be minimized by many well-established numerical algorithms such as Quadratic programming with linear constraints \cite{BenTal_Aharon_Nemirovski}, but the high-dimensionality of the practical problems (such as image deblurring) often precludes the use of such methods. However, because the functional in \eqref{l1-constraint} is a combination of convex-differentiable and convex term, one can make use of the proximal operators to find the derivative of the functional defined in \eqref{l1-constraint}. This leads to one of the most popular methods for minimizing \eqref{l1-constraint}, which falls in the class of \textit{iterative (shrinkage or) soft-thresholding algorithms (ISTA)}, where in each iteration the gradient of the differentiable term is projected and shrinked by a factor of $\lambda$ (known as \textit{soft-thresholding}), see \cite{Bruck_Ronald1977,Passty1979,Daubechies_Defrise_DeMol2004,Combettes_Wajs2005, Figueiredo_Nowak2003, Tikhonov_Arsenin1997, Vonesh_Unser2007}. Specifically, at the $m^{th}$-step the iterate $x_m^\delta$ is improved in the following manner
\begin{equation}
    x_{m+1}^\delta = S_{\lambda \tau}(x_m^\delta + 2 \tau A^*(b_\delta - Ax_m^\delta)),
\end{equation}
where $\tau$ is an appropriate stepsize and the operator $S_{\lambda \tau}: \mathbb{R}^n \rightarrow \mathbb{R}^n$ is the shrinkage operator defined as, for any $x \in \mathbb{R}^n$, $S_{\lambda \tau}(x) = (S_{\lambda \tau}(x)_i)_{i=1}^n \in \mathbb{R}^n$, where
\begin{equation}\label{proximalOP}
    S_{\lambda \tau}(x)_i = sgn(x_i)*(|x_i| - \lambda \tau)_+,
\end{equation}
where $sgn(x_i) = \frac{x_i}{|x_i|}$ and $(|x_i| - \lambda \tau)_+ = max\{|x_i| - \lambda \tau, 0\}$. So one can see that it's a forward-backward splitting method, where at each iteration the iterate ($x_m^\delta$) is first moved in the forward direction (i.e, $x_{m+1}^\delta = x_m^\delta + 2 \tau A^*(b_\delta - Ax_m^\delta)$) and then moved backward by the shrinking operator $S_{\lambda \tau}$, as defined in \eqref{proximalOP}, i.e., truncating out the \textit{``smaller values"} of $x_{m+1}^\delta$. One can also observe from \eqref{proximalOP} that (not only the regularization, but also) the sparsity of the recovered regularized solution ($x^\delta$) will be influenced by the regularization parameter ($\lambda$). Hence, for larger $\lambda$ the recovered solution will not only be overly-smoothed/regularized (or highly-biased) but also be sparser, where as for smaller $\lambda$ the recovered solution will be noisy (or high-variance) and dense (low sparsity). Even for an appropriate parameter value ($\lambda_0$) ISTA solution ($x_{\lambda_0}$) still tends to be more biased, due the soft-thresholding, for details see the references mentioned above. ISTA also suffers from a slower convergence rate of the descent process, which can be improved by using FISTA (fast iterative soft-thresholding algorithm) \cite{FISTA_Beck_Teboulle}, but it still doesn't compensate the biasness problem. Usually, the biasness problem is handled through a two-step process: (1) either ISTA or FISTA is used for variables/features selection and then, (2) an ordinary least-square (OLS) is performed on the selected variables and the corresponding matrix.

\subsection{A new iterative algorithm for sparsity:} In this paper, we propose a new (semi-) iterative regularization method which is not only as simple (if not simpler) as ISTA or FISTA from a computational viewpoint, but also improves the sparsity of the recovered solution significantly, even removing the biasness effect shown by the mentioned methods, i.e., improving the accuracy (or relative error) of the recovery. First note that, starting from the simple Landweber iterations \eqref{steepest descent} all the generalizations and extensions \eqref{Krylov-space method} developed in the (semi-) iterative regularization literature is focused on improving the \textit{speed of the convergence} of the descent process. However, they do not motivate sparsity during the descent (or minimization) process, even when started with a zero vector ($x_0^\delta = \bar{0} \in \mathbb{R}^n$). Hence, here we provide some different descent directions that advocate sparsity during the recovery process and thus, leading to a sparse regularized solution $x^\delta$, when started with $x_0^\delta = \bar{0} \in \mathbb{R}^n$.

The paper is organized in the following manner. In \S \ref{sparsity descent directions} we provide some descent directions that will not only impart sparsity but also has a regularizing effect in the recovery process. In \S \ref{Convergence and Implementation} we prove some results corresponding to the convergence of this method, when dealing with a noisy data, and the sparsity in the regularized solution. In \S \ref{Convergence and Implementation} we also compare and contrast our algorithm with ISTA, and present the advantages of our method over ISTA. In \S \ref{Numerical Results} we present numerical results to validate the developed theory and, we also compare the numerical results obtained using our method with the results obtained using certain standard regularization methods such as total-variation regularization, FISTA (fast iterative soft-thresholding algorithm), CGLS method, hybrid LSQR and many others, see Tables \ref{satellite error comparison}, \ref{hubble telescope error comparison} and \ref{table modified DP}. Furthermore, in \S \ref{modified DP}, we also provide a very effective and efficient (practical) stopping criterion to estimate an appropriate value of the regularization parameter (here, the iteration index) so as to obtain a regularized (sparse) solution, even when no noise information is available. We show that this technique is a significant improvement over the simpler Morozov's discrepancy principle.

\section{Sparsity motivating descent directions}\label{sparsity descent directions}
In this section, first we work with noise-free data ($g$) and then extend the results to noisy data ($g_\delta$). Note that if a functional $F$ is (strictly) convex then there is a unique minimizer (say $\varphi$), which can be attained in an iterative descent process. That is, if there is a sequence of functions $\{ \psi_m \}$ (or vectors $\{x_m\}$, in discretized case) such that the corresponding sequence $\{ F(\psi_m) \}$ is strictly decreasing ($F(\psi_{m+1}) < F(\psi_m)$), then the sequence $\{\psi_m\}$ \textit{``converges"} (in some sense)\footnote{the nature of convergence, i.e., either strongly or weakly, depends on the functional $F$.} to the minimizer $\varphi$. The most common descending direction at any step $m$ is the steepest-descent direction, which (for the functional defined in \eqref{LS solution}) is given by the (negative) gradient of the functional, i.e., $\xi_m = -2T^*(T\psi_m - g)$, see below.
As aforementioned, though it's the steepest-descent direction, it doesn't advocate sparsity in the recovery process. Therefore, if sparsity is the subject of interest, we need to find some other (appropriate) descent directions that will support sparsity in the descent process.

Observe that, for any functional $F$ we have the following Taylor's expansion
\begin{equation}\label{Taylor1}
    F(\psi + \tau h) = F(\psi) + \tau F'(\psi)[h] + O(\tau^2),
\end{equation}
or equivalently, for sufficiently small $\tau > 0$, we can have
\begin{equation}\label{Taylor2}
    F(\psi + \tau h) - F(\psi) \approx \tau F'(\psi)[h],
\end{equation}
where the differential $F'(\psi)$ for the functional as defined in \eqref{LS solution} is
\begin{align}\label{F differential}
    F'(\psi)[h] &= 2\Lprod{T^*(T\psi - g)}{h},\\
                &= 2\int_\Omega T^*(T\psi - g) h \; d\mu   \notag
\end{align}
and we denote $\lgrad{\psi}F := \xi = T^*(T\psi - g)$
as the $\Ltwo$-gradient of the functional $F$. Note that, considering $h=-\xi_m := -T^*(T\psi_m - g)$ (as the descent direction at the $m^{th}$-step) we have, from \eqref{Taylor2}, $F(\psi_m + \tau \xi_m) < F(\psi_m)$, i.e., starting from an initial choice $\psi_0$ we can have a sequence of functions $\{ \psi_m \}$ such that the sequence $\{ F(\psi_m) \}$ is (strictly) decreasing, and hence, the sequence $\{ \psi_m \}$ converges to the minimizer $\varphi$, see \cite{Engl+Hanke+Neubauer, Landweber, Bakushinsky+Goncharsky, Groetsch, Baumeister,Morozov_a}.

\subsection{Sparsity inducing descent directions:}\mbox{}\\
Splitting the integral in \eqref{F differential} into its distinct sub-domains, we have, for any $\lambda > 0$
\begin{align}\label{F' integral form}
    F'(\psi)[h] &= 2\int_\Omega T^*(T\psi - g)h \; d\mu\\
                &= 2 \left[ \int_{|T^*(T\psi - g)| > \lambda} T^*(T\psi - g)h\; d\mu + \int_{|T^*(T\psi - g)| \leq \lambda} T^*(T\psi - g)h \; d\mu \right] \notag\\
                &= 2 \left[ \int_{T^*(T\psi - g) > \lambda} T^*(T\psi - g)h\; d\mu + \int_{0 < T^*(T\psi - g) \leq \lambda} T^*(T\psi - g)h \; d\mu \right. \notag\\
                & \hspace{0.5cm}+ \left. \int_{T^*(T\psi - g) < -\lambda} T^*(T\psi - g)h\; d\mu + \int_{-\lambda \leq T^*(T\psi - g) < 0} T^*(T\psi - g)h \; d\mu \right]. \notag
\end{align}
The above expression \eqref{F' integral form} inspires us to formulate some descent directions that advocates sparsity in the reconstruction during the descent process.

\subsubsection{\textbf{$\lambda$-truncated gradient:}}\mbox{}\\
From the expression \eqref{F' integral form} if we instead choose
\begin{equation}\label{lambda truncated gradient}
    h_{m, \lambda} = -T^*(T\psi_m - g)\cdot \chi_{|T^*(T\psi_m - g)| > \lambda}
\end{equation}
(for an appropriate $\lambda > 0$) as the descent direction, then not only will the \textit{``smaller"} gradient values ($|T^*(T\psi_m - g)| \leq \lambda$) be discarded but also, we have $F'(\psi_m)[h_{m, \lambda}] < 0$, i.e., $F(\psi_{m+1}) < F(\psi_m)$, where $\psi_{m+1} = \psi_m + \tau_m h_{m,\lambda}$ and $\chi_{A}$ is the standard characteristic (or indicator) function, defined as,
\begin{equation}
    \chi_{A}(x) =
    \begin{cases}
        1, \; x \in A \\
        0, \; x \notin A.
    \end{cases}
\end{equation}
Therefore, starting from a zero function ($\psi_0 \equiv 0$) and using the \textit{$\lambda$-truncated gradient}, which we denote by $\xi_{m,\lambda} := h_{m,\lambda}$, as the descent direction one may expect to recover a sparse solution, this is proved below.

First, note that here the value of $\lambda$ plays an important role in balancing between the \textit{sparsity} and the \textit{accuracy} of the recovered solution, where by ``accuracy" we mean the \textit{relative error} of the recovery. If $\lambda$ is too large (specifically, $\lambda > ||T^*g||_\infty$), then $\psi_m \equiv 0$ for all $m \geq 0$, i.e., completely sparse, but (obviously) not a very good approximated solution for \eqref{T Gen.}. Now if $\lambda$ is too small, such that there is no truncation of the gradient at any step, then it will be a simple steepest-descent algorithm, which doesn't guarantee any sparsity, but results in the most accurate estimate for the solution of \eqref{T Gen.}. Now we show that one can also obtain a non-trivial sparse (approximate) solution for \eqref{T Gen.} with an appropriate choice of $\lambda > 0$. Note that, for a sequence $\{ \psi_m \}$ if the corresponding sequence $\{ F(\psi_m) \}$ tends to zero ($F(\psi_m) \rightarrow 0$), then for every $\epsilon > 0, \exists M(\epsilon) \in \mathbb{N}$ such that $|F(\psi_m)| \leq \epsilon$, for all $m \geq M(\epsilon)$. In other words, for all $m \geq M(\epsilon)$, $||T\psi_m - g||_2^2 \leq \epsilon$, and since $T$ is a bounded linear operator (and so is $T^*$), we have $||T^*(T\psi_m - g)||_2^2 \leq ||T^*||_2^2 \; ||T\psi_m - g||_2^2 \leq ||T^*||_2^2 \; \epsilon$. That is, for any $\lambda > 0$ there is a $M(\lambda) \in \mathbb{N}$ such that the gradient-norm $||\xi_m = T^*(T\psi_m - g)||_2 \leq \lambda$ (implying $||\xi_m \cdot \chi_A||_2 \leq \lambda$, for any $A \subset \Omega_{\psi_m}$ with $\mu(A) > 0$), for all $m \geq M(\lambda)$, and hence, if the recovered solution ($\psi_{M(\lambda)}$) at $M(\lambda)$ is sparse, then so are all $\psi_m$'s for $m \geq M(\lambda)$, in fact, $\psi_m \equiv \psi_{M(\lambda)}$ for all $m \geq M(\lambda)$. Now the existence of such a $\lambda$ can be shown as below.

Note that, for the first update, without any truncation of the gradient,
\begin{equation}
    (\psi_0 \equiv 0, \xi_0 = T^*g) \xrightarrow{m=1 (no-truncation)} (\psi_1 = \tau_1 T^*g, \xi_1 = T^*( \tau_1 TT^*g - g)),
\end{equation}
where $\tau_1 > 0$ is the step-size, we have $||T^*(\tau_1 TT^* - I)g||_\infty < ||T^*g||_\infty$, since $T$ is linear and bounded. Now, for any choice of $\lambda$ satisfying $||T^*(\tau_1 TT^* - I)g||_\infty < \lambda < ||T^*g||_\infty $, leads to a non-trivial sparse $\psi_{1,\lambda} := \tau_1 (T^*g)\cdot (\chi_{|T^*g| > \lambda})$, since $\mu(x:|(T^*g)(x)| \leq \lambda) > 0$ and $\mu(x:\lambda < |(T^*g)(x)| \leq ||T^*g||_\infty) > 0$, otherwise $\lambda = ||T^*g||_\infty$. Also, we have $F(\psi_{1,\lambda}) < F(\psi_0)$, as $F(\psi_0) = ||T^*g||_2^2$ and $F(\psi_{1,\lambda}) = ||T^*(\tau_1 TT^*g\cdot (\chi_{|T^*g| > \lambda}) - g)||_2^2$, i.e., $\psi_{1,\lambda}$ is an improvement over $\psi_0$ in the approximation of the solution of \eqref{T Gen.}. Furthermore, since the gradient at $m=1$, with $\psi_{1,\lambda} = \tau_1 (T^*g)\cdot (\chi_{|T^*g| > \lambda})$, is $\xi_1 = T^*(T\psi_{1,\lambda} - g)$, we have $||\xi_1||_\infty \leq ||T^*(\tau_1 TT^* - I)g||_\infty < \lambda$, implying the $\lambda$-truncated gradient at $m=1$ will be $\xi_{1,\lambda} \equiv 0$, and hence, $\psi_m = \psi_{1,\lambda}$ for all $m \geq 2$. In other words, it says, choosing an appropriate $\lambda$ would update the initial choice $\psi_0 \equiv 0$ to get a sparse $\psi_{1,\lambda} \not \equiv 0$, such that, it's a better approximation for the solution of \eqref{T Gen.}, and $\psi_{m,\lambda} = \psi_{1,\lambda}$, for all $m \geq 2$, i.e., the sparsity is not destroyed, though there is no further improvements.

Therefore, when $\lambda$ decreases the sparsity of the recovered solution also decreases, but the accuracy of the recovery increases, and vice-versa, i.e., the parameter value $\lambda$ balances between the sparsity and the accuracy of the recovery. An algorithm to choose an appropriate value of $\lambda$ such that the constructed sequence $\{ \psi_{m,\lambda} \}$ \textit{``converges"} to a regularized (sparse) solution ($\varphi^\delta$) of \eqref{T Gen.} for a noisy data ($g_\delta$) is presented in \S \ref{Convergence and Implementation} and, its numerical implementation and results shown in \S \ref{Numerical Results}.

\begin{remark}\label{lambda regularization effect}
Note that, the parameter $\lambda$ not only serves as a sparsity parameter, but also as a regularization parameter, since (as explained above) for any $\lambda > 0$ there is a $M(\lambda) \in \mathbb{N}$ such that $\psi_m = \psi_{M(\lambda)}$ for all $m \geq M(\lambda)$. In other words, when dealing with a noisy data ($g_\delta$), where the regularization is achieved by terminating the descent process at an appropriate iteration index, the regularization effect can also be achieved by using the $\lambda$-truncated gradient, for an appropriate $\lambda > 0$.
\end{remark}

\begin{remark}\label{split responsiblities}
As stated in remark \ref{lambda regularization effect}, the value $\lambda$ can play the role of the sparsity parameter as well as the regularization parameter, in the $\lambda$-truncated gradient descent algorithm, and hence, has double the responsibilities. Now, since in a traditional (semi-) iterative regularization method the regularization effect is handled by the iteration index ($m$), we would like to distribute the responsibility of regularization and sparsity to the parameters $m$ (iteration index) and $\lambda$, respectively. This will enable us to increase the level of sparsity in the recovered regularized solution without affecting the accuracy of the recovery (i.e. the regularizing effect). One way to achieve that is through the approach stated in \ref{Convergence and Implementation}, and/or having a sequence of iteration-dependent $\lambda$ values, i.e. $\lambda_m$-truncated gradient, which is discussed below.
\end{remark}

\subsubsection{\textbf{$\lambda_m$-truncated gradient:}}\mbox{}\\
In order to split the responsibility of regularization and sparsity between the parameters $m$ (iteration index) and $\lambda$, respectively, as discussed in remark \ref{split responsiblities}, we can use a sequence $\{\lambda_m > 0\}$ of truncation values, instead of a constant $\lambda > 0$, such that in each iteration step ($m$) we can truncate the $m^{th}$ gradient ($\xi_m$) based on $\lambda_m$, i.e., our $m^{th}$ (sparsity inducing) descent direction ($h_{m,\lambda_m}$) will be now
\begin{equation}\label{lambdam truncated gradient}
    h_{m,\lambda_m} = -T^*(T\psi_m - g)\cdot \chi_{|T^*(T\psi_m - g)| > \lambda_m},
\end{equation}
and we denote $\xi_{m,\lambda_m} := h_{m,\lambda_m}$ as the $\lambda_m$-truncated gradient. One of the choices for the sequence $\{ \lambda_m \}$ can be based upon the sequence values $\{ || \xi_m = T^*(T\psi_m - g) ||_\infty \}$, such as, $\lambda_m = \alpha\% * ||\xi_m||_\infty$. In this case, one updates the $m^{th}$ iterate $\psi_{m,\lambda_{m}}$ to get $\psi_{m+1, \lambda_{m+1}}$, but only corresponding to top $(1-\alpha) \%$ of its gradient $\xi_{m} = T^*(T\psi_{m, \lambda_{m}} - g)$, that is, $\psi_{m+1, \lambda_{m+1}} = \psi_{m,\lambda_{m}} - \tau_{m+1} \xi_{m, \lambda_{m+1}}$.
We implemented this choice of $\{ \lambda_m \}$ in our numerical experiments, for different values of $\alpha$, and the results are presented in \S \ref{Numerical Results}, and also a comparison with the $\lambda$-truncated gradient approach.

\begin{remark}
Note that, in the $\lambda$-truncated gradient method one has to a-priori fix the $\lambda$-value for all the iterations and, since it also influence the regularization (or accuracy) of the recovery process, one cannot have a large $\lambda$-value, i.e., one has to compromise between the sparsity and regularization, to certain extent. Where as, in the $\lambda_m$-truncated gradient method one can have larger $\lambda_m$ values for the first few iterations (since the $||\xi_m||_\infty$ values are large in the initial stage) imposing a lot of sparsity, in addition to not compromising the accuracy of the recovery, as the regularization (or accuracy of the recovery) can be achieved through terminating the iterations independent of the $\lambda_m$-values, say discrepancy principle. That is, the stopping criterion for the descent process can be made independent of $\lambda_m$, say depending on the noise level $\delta$, leading to greater degrees of freedom for the regularization parameter ($M(\delta)$, instead of $M(\lambda)$) and the sparsity parameter ($\lambda_m$). In \S \ref{Numerical Results} we compare the sparsity and the accuracy of the results obtained by $\lambda$-truncated and $\lambda_m$-truncated gradient methods.
\end{remark}

\subsubsection{\textbf{Discretized truncated gradients:}}\mbox{}\\
Note that, all of the above developed theories and results are also valid in the discretized settings $(A,x,b)$, when the function $\psi_m$ is repalced by the vector $x_m$ and the continuous measure $\mu$i is replaced by the counting measure. In fact, it's much easier when dealing in a discretized system, for example, from the $m^{th}$ gradient $d_m = A^*(Ax_m - b) \in \mathbb{R}^n$ we can have
\begin{enumerate}
    \item  \textbf{the discrete $\lambda$-truncated gradient}, $d_{m,\lambda} = ((d_{m,\lambda})_i)_{i=1}^n$, as
    \begin{equation}
        (d_{m,\lambda})_i =
        \begin{cases}
            (d_m)_i \; &if \; |(d_m)_i| > \lambda\\
            0 \; &if \; |(d_m)_i| \leq \lambda
        \end{cases}
    \end{equation}

    \item \textbf{the discrete $\lambda_m$-truncated gradient}, $d_{m,\lambda_m} = ((d_{m,\lambda_m})_i)_{i=1}^n$, as
    \begin{equation}
        (d_{m,\lambda_m})_i =
        \begin{cases}
            (d_m)_i \; &if \; |(d_m)_i| > \lambda_m\\
            0 \; &if \; |(d_m)_i| \leq \lambda_m
        \end{cases}
    \end{equation}

    \item \textbf{the $k$-truncated gradient:}

    In addition, when working in a discretized setting where the discretized solution space is a finite high-dimensional space, say $\mathbb{R}^n$ for large $n$, but the solution ($\hat{x}$) of \eqref{Ax=b} has a sparse structure to it, i.e., $||\hat{x}||_0 \leq k << n$, then we also have the luxury to construct a \textit{$k$-truncated gradient} ($d_{m,k}$) by simply projecting the ($n-k$)-smaller values of the gradient $d_m \in \mathbb{R}^n$ to $\bar{0} \in \mathbb{R}^{(n-k)}$, i.e., $d_{m,k} = (d_{m,k})_{i=1}^n \in \mathbb{R}^n$ such that $||d_{m,k}||_0 \leq k$, where
    \begin{equation}
        (d_{m,k})_i =
        \begin{cases}
                (d_m)_i \; & if \; (d_m)_i \geq \mbox{k-th max of } d_m\\
                0 \; & otherwise.
        \end{cases}
    \end{equation}
\end{enumerate}

\begin{remark}
Note that, when using the discrete $\lambda$-truncated gradient ($d_{m,\lambda}$) or $\lambda_m$-truncated gradient ($d_{m,\lambda_m}$), the sparsity of the gradients is dynamically bounded, i.e., $||d_{m,\lambda}||_0 \leq k_m$ or $||d_{m,\lambda_m}||_0 \leq k_m$, where the sparsity bound ($k_m$) of the truncated gradients varies dynamically depending on the $\lambda$ or $\lambda_m$ values and the gradient $d_m$, i.e., $k_m = k(\lambda, d_m)$ or $k_m = k(\lambda_m, d_m)$, respectively. Where as, the bound is fixed for $k$-truncated gradient ($d_{m,k}$), i.e., $||d_{m,k}||_0 \leq k$, irrespective of the iterations. However, this doesn't imply that the sparsity of the iterates $x_{m,k}$ are $k$-bounded, i.e., $||x_{m,k}||_0 \not \leq k$, as $x_{m,k}$ are the updates of the previous iterates ($x_{m,k} = x_{m-1,k} - \tau_m d_{m,k}$) and the sparsity of $x_{m-1,k}$ may not be $k$-bounded.
\end{remark}

\begin{remark}
Also note that the $k$-truncated gradient can be derived from the $\lambda_m$-truncated gradient, by choosing $\lambda_m = \min\{k, ||d_m||_0\}$-th max of $d_m$, for the positivity of $\lambda_m$. However, we mostly determine the values of $\lambda_m$ as some percentage of $||d_m||_\infty$, i.e., $\lambda_m = \alpha\% * ||d_m||_\infty$, for some $0 \leq \alpha \leq 100$.
\end{remark}

\subsubsection{\textbf{$k\lambda_m$-truncated gradient:}}\mbox{}\\
Now one can make use of both the gradients $\lambda_m$-truncated and $k$-truncated by either having the
\begin{enumerate}
    \item \textbf{min-k$\lambda_m$-truncated gradient:} ($d_{m,k\lambda_m^-}$) which can be defined as
    \begin{equation}
        d_{m,k\lambda_m^-} =
        \begin{cases}
             d_{m,\lambda_m} \; & \; if \; ||d_{m,\lambda_m}||_0 \leq ||d_{m,k}||_0\\
             d_{m,k} \; & \; if \; ||d_{m,\lambda_m}||_0 > ||d_{m,k}||_0,
        \end{cases}
    \end{equation}
    which increases the sparsity of the recovery process, more than the individual ones, or

    \item \textbf{max-k$\lambda_m$-truncated gradient:} ($d_{m,k\lambda_m^+}$) which can be defined as
    \begin{equation}
        d_{m,k\lambda_m^+} =
        \begin{cases}
             d_{m,k} \; & \; if \; ||d_{m,\lambda_m}||_0 \leq ||d_{m,k}||_0\\
             d_{m,\lambda_m} \; & \; if \; ||d_{m,\lambda_m}||_0 > ||d_{m,k}||_0,
        \end{cases}
    \end{equation}
    which increases the accuracy of the recovery process, more than the individual ones.
\end{enumerate}
The numerical implementations, results and comparison (with the individual gradients: $\lambda_m$-truncated or $k$-truncated gradient) is presented in \S \ref{Numerical Results}.

\section{Implementation algorithms and convergence results}\label{Convergence and Implementation}
\subsection{Convergence:}
In this section we describe how to implement the above truncated gradients, when dealing with a noisy data ($b_\delta$), such that the sequence $\{ x_{m,\lambda|\lambda_m|k|k\lambda_m^-|k\lambda_m^+}^\delta \}$ (or $\{ \psi_{m,\lambda|\lambda_m}^\delta \}$) converges to a regularized (sparse) solution $x^\delta$ (or $\varphi^\delta$). First note that, as discussed in \S \ref{semi-iterative methods}, the key to obtain a regularized solution in a (semi-) iterative regularization method (when dealing with a noisy $b_\delta$) is to terminate the iterations at an appropriate index $m_0 < \infty$. A typical choice of the stopping condition is to terminate the minimization process when the residue ($Ax_m^\delta - b_\delta$) has achieved a certain level of accuracy, i.e., when $||b - b_\delta|| \leq \delta$ is known, then terminate the iteration at $M(\delta)$ such that $||Ax_{M(\delta)}^\delta - b_\delta|| \leq \eta \delta$, for $\eta > 1$, this is known as \textit{Morozov's discrepancy principle}, see \cite{Morozov_a} for the convergence analysis. In other words, it says, a reconstructed solution ($x^\delta$) is a good approximate of the true solution $(\hat{x})$ if the image of $x^\delta$, under the forward operator $A$, fits the noisy data $b_\delta$ appropriately, i.e., $Ax^\delta$ doesn't underfit or overfit $b_\delta$.

For the $\lambda$-truncated gradient, note that $F(x_{m+1,\lambda}^\delta) < F(x_{m,\lambda}^\delta)$, for all $m \leq M(\lambda)$, where $M(\lambda) \in \mathbb{N}$ is the smallest integer such that $||d_{M(\lambda)}^\delta||_\infty \leq \lambda$, and $F(x_{m,\lambda}^\delta) = F(x_{M(\lambda),\lambda}^\delta)$ for all $m > M(\lambda)$, since the $\lambda$-truncated gradient $d_{M(\lambda),\lambda}^\delta \equiv 0$. This implies the iterates $x_{m,\lambda}^\delta$ get updated for all $m \leq M(\lambda)$ and remain unchanged after $M(\lambda)$, i.e., $x_{m,\lambda}^\delta = x_{M(\lambda),\lambda}^\delta$ for all $m > M(\lambda)$. Hence, an appropriate $\lambda$ should be such that $M(\delta) \leq M(\lambda)$, implying the recovered solution $x_{M(\delta),\lambda}^\delta$ (being both sparse and regularized) as a good approximate of the true solution ($\hat{x}$), since $||Ax_{M(\delta), \lambda}^\delta - b_\delta|| \leq \delta$.

A similar argument also works for the $\lambda_m$-truncated gradient, for example, if $\lambda_m$ are chosen such that $||\{ d_m^\delta > \lambda_m \}||_0 \geq 1$, then the iterate $x_{m,\lambda_m}^\delta$ is updated to $x_{m+1, \lambda_{m+1}}^\delta = x_{x, \lambda_m}^\delta - \tau_{m+1}d_{m+1,\lambda_{m+1}}$ with $F(x_{m+1, \lambda_{m+1}}^\delta) < F(x_{m,\lambda_m}^\delta)$, i.e., an improvement over the previous iterate. Again, if we are able to satisfy the stopping criterion (i.e., $||Ax_{M(\delta),\lambda_{M(\delta)}}^\delta - b_\delta|| \leq \delta$, for some $M(\delta) < \infty$), then we are to recover a regularized (sparse) solution $x_{M(\delta), \lambda}^\delta$, which is a good approximation to the true solution $\hat{x}$. A choice for such a sequence can be $\{ \lambda_m := \alpha\% * ||d_m^\delta||_\infty \}$, for $1 \leq \alpha \leq 100$, where one can see that the sparsity increases upon increasing $\alpha$, but may compromise the accuracy for large $\alpha$, typically $\alpha < 50$.

Similarly, one can also prove the convergence to a regularized (sparse) solution, which is a good approximation to the true solution $\hat{x}$, when dealing with the $k$-truncated or $k\lambda_m$-truncated gradients.

\subsection{Implementation Algorithm:} As one can clearly see, the numerical implementation of the above developed theories is very simple, as one simply truncates the gradient obtained at any stage, based on certain criterion. Typically, first, one recovers a regularized solution without truncating the gradients, i.e., a simple steepest-descent (gradient) minimization. This gives an idea of the number of iterations needed and a base regularized solution to compare with. Then, one performs the truncated-gradients descent algorithms and compares the (sparse) regularized solution with the already obtained regularized solution, just in case. Also, it's advisable to start with smaller $\lambda$ or $\lambda_m$ values and larger $k$ values, so that it doesn't compromise the accuracy of the recovery by too much, and then, can either increase ($\lambda$ or $\lambda_m$) or decrease ($k$) the truncation values to increase the sparsity level in the reconstruction process.

\subsection{Comparison with ISTA:} Comparing with the Iterative Soft/Shrinkage Threshold Algorithm (ISTA), which is described earlier, one can see that there is some similarity between these two approaches, though they are still inherently different. For example, in ISTA the iterate $x_m^\delta$ is improved in a forward-backward splitting manner to get the new iterate $x_{m+1}^\delta$, i.e.,
\begin{equation}
    x_m^\delta \xrightarrow{forward} x_{m+1}^\delta = x_m^\delta - 2\tau A^*(Ax_m^\delta - b_\delta) \xrightarrow{backward} x_{m+1}^\delta = S_{\lambda \tau}(x_{m+1}^\delta),
\end{equation}
where the shrinkage operator $S_{\tau \lambda}$, as is defined in \eqref{proximalOP}, not only discards the smaller $x_{m+1}^\delta$ values (i.e., $|(x_{m+1}^\delta)_i| \leq \tau \lambda \Rightarrow (x_{m+1}^\delta)_i=0$), but also shrink the larger $x_{m+1}^\delta$ values (i.e., $ |(x_{m+1}^\delta)_i| > \lambda \Rightarrow sgn((x_{m+1}^\delta)_i)*(|(x_{m+1}^\delta)_i| - \tau \lambda) $).

Where as, in our method the iterate $x_m^\delta$ is also improved in a forward-backward splitting fashion to get the new iterate $x_{m+1}^\delta$, but the forward-backward splitting is applied to the gradient ($\xi_m^\delta$) instead, i.e.,
\begin{equation}
    (x_m^\delta,0) \xrightarrow{forward} (x_m^\delta, \xi_{m}^\delta) \xrightarrow{backward} (x_m^\delta, \xi_{m,\lambda}^\delta) \mapsto x_{m+1}^\delta = x_m^\delta - 2\tau \xi_{m,\lambda}^\delta,
\end{equation}
where $\xi_{m,\lambda}^\delta$, as defined in \eqref{lambda truncated gradient}, truncates the smaller $\xi_m^\delta$ values (i.e., $|(\xi_{m}^\delta)_i| \leq \lambda \Rightarrow (\xi_{m,\lambda}^\delta)_i=0$, otherwise $(\xi_{m,\lambda}^\delta)_i = (\xi_{m}^\delta)_i$), that is, thresholding the gradients instead of the iterates.

The most significant advantage of our method over ISTA or FISTA is that of the (overall) computational time. In ISTA or FISTA, to obtain a regularized solution ($x_{\lambda_0}^\delta$), one has to \textbf{completely minimize} the functional \eqref{l1-constraint} for a \textbf{certain number} of $\lambda$ values ($\lambda \in [\lambda_{min}, \lambda_{max}]$), i.e., one has to minimize the functional $F_\lambda$ such that $F_\lambda(x_{M(\epsilon)}^\delta) < \epsilon$, for some accuracy threshold parameter $\epsilon > 0$ (usually, leading to a very high $M(\epsilon)$ value for small $\epsilon$), and also for that many $\lambda$ values. In comparison, in our method one has to perform only certain number of iterations such that $F(x_{M(\delta)}^\delta) \leq (\eta \delta)^2$ (usually, $M(\delta) << M(\epsilon)$, much smaller) and also has to perform the minimization only once (or maybe for a few $\lambda$ values, just to be conservative). Furthermore, we also show that (see Tables \ref{satellite error comparison} and \ref{hubble telescope error comparison}) the sparsity level in our recovered solution far exceeds than that obtained by any other regularization methods.

\begin{remark}
In our numerical experiments (presented in \S \ref{Numerical Results}) we simply modified the steepest-descent algorithm by replacing the gradients $\xi_m^\delta$ with the truncated-gradient $\xi_{m,\lambda_m}^\delta$ the simple truncated gradients as the descent directions, i.e., in \eqref{steepest descent} we used $\xi_{m}^\delta = \xi_{m,\lambda_m}^\delta$ and $\psi_m^\delta = \psi_{m,\lambda_m}^\delta$. However, as mentioned in \S \ref{semi-iterative methods}, one can further boost the rate of convergence by opting the accelerated Landweber iterations (generalized semi-iterative methods or the Krylov-subspace methods, and even conjugate-gradient methods), instead of the simple Landweber iterative approach. Using a boosted algorithm really speeds up the descent process.
\end{remark}

\section{Numerical Results}\label{Numerical Results}
To test the numerical viability and efficiency of the developed algorithm a MATLAB code was written to solve the image deblurring inverse problem. The mechanism behind blurring an image results in a Fredholm-type integral equation of first kind, with the kernel as a Gaussian spread function, i.e., the operator describing the underlying physical process can be written as, for any $s \in \Omega_1 \subset \mathbb{R}^2$,
\begin{equation}\label{fredholm equation}
g(s) = (T\varphi)(s) := \int_\Omega K(s,t) \varphi(t) dt, \;\;\; \Omega \subset \mathbb{R}^2,
\end{equation}
where the kernel $K(s,t) = \frac{1}{\sqrt{2\pi \sigma^2}} e^{\frac{|s-t|^2}{2\sigma^2}}$, $\sigma$ determining the spread and the ill-posedness of the problem (ill-posedness increases with the increase in the $\sigma$ value, here we considered $\sigma=10$, highly ill-posed). Such integral operator equations \eqref{fredholm equation} are classical inverse problems and can be quite ill-posed. The operator equation \eqref{fredholm equation} is then discretized using either Galerkin or Nystr$\ddot{o}$m methods to yield the following linear discrete ill-posed problem
\begin{equation}
    Ax = b,
\end{equation}
where the matrix $A \in \mathbb{R}^{m\times n}$ (blurring matrix) is the discretized representation of the operator $T$ and the vectors $x \in \mathbb{R}^n$ (true image, Figure \ref{deblurring satellite and hubble telescope}), $b \in \mathbb{R}^m$ (blurry image) are the discretized source and effect functions. MATLAB functions in \cite{Hansen_MATLAB2020, Hansen_IRtools, Hansen2007} determine the discretizations $A \in \mathbb{R}^{n\times n}$, the scaled discrete approximations of $x \in \mathbb{R}^n$ and $b := Ax \in \mathbb{R}^m$, where we considered $n = m = 65536$, i.e., a $256 \times 256$ image. Specifically, we used the routine $blur()$, to generate the deblurring matrix $A$, and the routine $PRblur()$, to generate the true image, from \cite{Hansen_MATLAB2020,Hansen2007, Hansen_IRtools}. To test the stability of the method, we add a normally distributed (with mean zero) error vector $\epsilon_\delta \in \mathbb{R}^m$ to the original $b$ to get a perturbed vector $b_\delta \in \mathbb{R}^m$ (blurred noisy image, Figure \ref{deblurring satellite and hubble telescope}), such that the relative error $\left( \frac{||\epsilon_\delta||_2}{||b||_2} \right)$ in the noisy data is around $10\%$, unless otherwise stated. In particular, when using the discrepancy principle to terminate the descent process we assume $\delta = ||\epsilon_\delta||_2$ to be known.

We also compare the results obtained using our method (with different truncated gradients) with the results obtained using certain standard iterative regularization methods, such as (i) Least squares minimization with $\mathcal{L}^1$-norm (IRell1), (ii) (heuristic) total variation penalization (IRhtv), (iii) Conjugate Gradient (CGLS), (iv) FISTA for constrained LS (IRfista), (v) modified residual norm steepest descent (IRmrnsd), (vi) enriched CGLS (IRenrich), (vii) $\mathcal{L}^1$-norm penalization term (IRirn), (viii) hybrid FGMRES for enforcing $\mathcal{L}^1$-norm penalization (IRhybrid\_fgmres), (ix) hybrid GMRES for square systems (IRhybrid\_gmres), (x) hybrid LSQR (IRhybrid\_lsqr), (xi) restarted Krylov subspace method (IRrestart),  where the routines for the above mentioned methods are taken from \cite{Hansen_IRtools}, with the stopping rule as the discrepancy principle (i.e., in the `options' settings for terminating the iterations the `stoprule' is set to be `DP'). The comparisons can be seen in Tables \ref{satellite error comparison} and \ref{hubble telescope error comparison}.

\begin{example}\label{example xmin=0 constraint}
In the first case, we perform a constraint minimization of the functional $F$ defined in \eqref{LS solution}, i.e., $\min_{x \geq 0} F(x)$, since we consider the pixel values to be non-negative. This is achieved by simply projecting the negative values of the iterates to zero during the descent process, i.e., $x_m^\delta = ((x_m^\delta)_i)_{i=1}^n$ such that $(x_m^\delta)_i = \max \{ (x_m^\delta)_i, 0 \}$. We do the same for the above mentioned other iterative methods by turning the `nonnegativity', `on' and having `xMin', `0' in the `options' settings. The recovered deblurred image is shown in Figure \ref{satellite_recovery} and, the relative errors and sparsity comparisons are presented in Table \ref{satellite error comparison}. Note that, $\hat{x} \in \mathbb{R}^{n}$, for $n=65536$, and the sparsity of the true image ($\hat{x}$) is 58858. We can see from Table \ref{satellite error comparison}, that using truncated gradients, we are able to attain sparsity level much better than most of the other methods. Also, note that the truncated gradient $d_{m,k\lambda_m^-}^\delta$ always yields better sparsity level than using the corresponding individual ones ($d_{m,\lambda_m}^\delta$ or $d_{m,k}^\delta$). Hence, if some a-priori knowledge about the sparsity of the true solution ($\hat{x}$) is known, then one can combine it with the $\lambda_m$-truncated gradients to significantly improve the sparsity of the recovery.
\end{example}

\begin{example}\label{example satellite, xmin<0}
Note that, in Example \ref{example xmin=0 constraint}, the sparsity in the recovery is also heavily influenced by the non-negativity constraint ($x \geq 0$), since the negative co-efficients of any iterate $x_m^\delta$ is projected to zero, i.e., $(x_{m}^\delta)_i = \max\{ (x_m^\delta)_i, 0 \}$. Hence, to observe the real sparsity effect of these methods we instead have the following constraint minimization: $\min_{x \geq -10^{-100}} F(x)$, i.e., the co-efficients of the iterate $x_m^\delta = ((x_m^\delta)_i)_{i=1}^n$ are modified as $(x_m^\delta)_i = \max \{ (x_m^\delta)_i, -10^{-100} \}$. We do the same for the other mentioned iterative methods by turning the `nonnegativity', `off' and having `xMin', `$-10^{-100}$' in the `options' settings. From a computational viewpoint it does not effect the accuracy (or relative error) of the recovery (since $-10^{-100} \approx 0$), as can be seen in Table \ref{satellite error comparison}, but we do see a significant lapse in the sparsity level, implying that the sparsity in those methods arose from the constraint. However, we notice that using our method, though the sparsity level decreases, the recovered solution is still quite sparse, especially for high truncation values.
\end{example}

\begin{example}\label{example hubble telescope, xmin=0}
In Example \ref{example xmin=0 constraint} the true image $\hat{x} \in \mathbb{R}^{65536}$ is too sparse (sparsity of $\hat{x}$ is 58858), and hence, are the recovered regularized sparse solutions. In this example, we repeat the same procedure, constraint minimization: $\min_{x\geq 0}F(x)$ as performed in example \ref{example xmin=0 constraint}, but on a non-sparse true image ($\hat{x} \in \mathbb{R}^{65536}$, with \# of zeroes in $\hat{x}$ as 1) with lot of smaller values, i.e., we consider the Hubble space telescope, see Figure \ref{hubble_true}. The recovered deblurred image is shown in Figure \ref{hubble_recovery} and, the relative errors and sparsity comparisons are presented in Table \ref{hubble telescope error comparison}. Note that, unlike the truncated gradient $d_{m,k\lambda_m^-}^\delta$ which supports sparsity, the truncated gradient $d_{m,k\lambda_m^+}^\delta$ supports the accuracy of the recovery and yields better results than using the corresponding individual ones ($d_{m,\lambda_m}^\delta$ or $d_{m,k}^\delta$).
\end{example}

\begin{example}\label{example hubble telescope, xmin<0}
In this example, we consider the non-sparse true image from Example \ref{example hubble telescope, xmin=0} and perform the modified constraint minimization: $\min_{x \geq -10^{-100}}$, as done in Example \ref{example satellite, xmin<0}. The recovered deblurred image and, relative errors and sparsity comparisons are shown in Table \ref{hubble telescope error comparison}.
\end{example}

\begin{remark}
Note that, in all of the above examples, we presented the results corresponding to only two truncation values that are widely separated. Now, one take even finer truncation values to either increase the sparsity of the recovery or the accuracy of it, depending on the scenario.
\end{remark}

\begin{figure}[ht]
    \centering
    \begin{subfigure}{0.4\textwidth}
        \includegraphics[width=\textwidth]{satellite_true}
        \caption{True Satellite Image}
        \label{satellite_true}
    \end{subfigure}
    \begin{subfigure}{0.4\textwidth}
        \includegraphics[width=\textwidth]{hubble_true}
        \caption{True Hubble Telescope}
        \label{hubble_true}
    \end{subfigure}
    \begin{subfigure}{0.4\textwidth}
        \includegraphics[width=\textwidth]{satellite_noisy.eps}
        \caption{Satellite noisy blurred}
        \label{satellite_noisy}
    \end{subfigure}
    \begin{subfigure}{0.4\textwidth}
        \includegraphics[width=\textwidth]{hubble_noisy.eps}
        \caption{Hubble Telescope noisy}
        \label{hubble_noisy}
    \end{subfigure}
    \begin{subfigure}{0.4\textwidth}
        \includegraphics[width=\textwidth]{satellite_recovery.eps}
        \caption{Satellite deblurred}
        \label{satellite_recovery}
    \end{subfigure}
    \begin{subfigure}{0.4\textwidth}
        \includegraphics[width=\textwidth]{hubble_recovery.eps}
        \caption{Hubble deblurred}
        \label{hubble_recovery}
    \end{subfigure}
    \caption{Deblurring of the blurred and noisy Satellite and Hubble telescope images, Examples \ref{example xmin=0 constraint}, \ref{example satellite, xmin<0}, \ref{example hubble telescope, xmin=0} and \ref{example hubble telescope, xmin<0}.}
    \label{deblurring satellite and hubble telescope}
\end{figure}

\begin{table}[ht]
    \centering
    \begin{tabular}{|p{2.6cm}||p{2.1cm}|p{2.3cm}||p{2.1cm}|p{2cm}|}
    \hline
    \multicolumn{5}{|c|}{Relative errors using different regularization methods, for $10\%$ data noise.}\\
    \hline
     & Constraint: & $x \geq 0$ & Constraint & $x \geq -10^{-100}$ \\
    \hline
    Methods & Relative error & Sparsity (\# 0s) & Relative error & Sparsity\\
    \hline
    IRell1 & 0.2937 & 0 & 0.2937 & 0 \\
    \hline
    IRhtv & 0.2138 & 53634 & 0.3769 & 0 \\
    \hline
    IRcgls & 0.2250 & 0 & 0.2250 & 0 \\
    \hline
    IRfista & 0.2096 & 41189 & 0.2096 & 0 \\
    \hline
    IRmrnsd & 0.1972 & 0 & 0.1972 & 0 \\
    \hline
    IRenrich & 0.2249 & 0 & 0.2249 & 0 \\
    \hline
    IRirn & 0.2375 & 7389 & 0.2466 & 0 \\
    \hline
    IRhy.\_fgmres & 0.2937 & 0 & 0.2937 & 0 \\
    \hline
    IRhy.\_gmres & 0.2334 & 0 & 0.2334 & 0 \\
    \hline
    IRhy.\_lsqr & 0.2736 & 0 & 0.2736 & 0 \\
    \hline

    {\color{red} Truncated-gradients:}  \\
    \hline
    $\xi_{m,\lambda} = \xi_m$, $\lambda = 0$, no truncation & 0.2105 & 34698 & 0.2105 & 0 \\
    \hline
    $\xi_{m,\lambda}$=$\xi_m\chi_{|\xi_m| > |\lambda|}$ $\lambda = 10^{-5}$ & 0.2082 & 43096 & 0.2082 & 25975\\
    \hline
    $\xi_{m,\lambda}$=$\xi_m\chi_{|\xi_m| > |\lambda|}$ $\lambda = 2*10^{-5}$ & 0.2002 & 52670 & 0.2002 & 45751\\
    \hline
    $\lambda_m$=$\alpha\%||\xi_m||_\infty$, $\alpha =10$ & 0.2107 & 47050 & 0.2107 & 28446\\
    \hline
    $\lambda_m$=$\alpha\%||\xi_m||_\infty$, $\alpha=40$ & 0.2043 & 57149 & 0.2043 & 55556\\
    \hline
    $k = 32768 (\approx 50\%*n)$ & 0.2094 & 40584 & 0.2094 & 15302\\
    \hline
    $k = 6554 (\approx 10\%*n)$ & 0.2023 & 54725 & 0.2023 & 48526\\
    \hline
    $k\lambda_m^-$, k=32768, $\lambda_m$=$40\% ||\xi_m||_\infty$ & 0.2043 & 57149 & 0.2043 & 55556\\
    \hline
    $k\lambda_m^-$, k=6554, $\lambda_m$=$10\%||\xi_m||_\infty$ & 0.2023 & 54725 & 0.2023 & 48526\\
    \hline
    \end{tabular}
    \caption{Image deblurring, Examples \ref{example xmin=0 constraint} and \ref{example satellite, xmin<0}}
    \label{satellite error comparison}
\end{table}

\begin{table}[ht]
    \centering
    \begin{tabular}{|p{2.6cm}||p{2.1cm}|p{2.3cm}||p{2.1cm}|p{2cm}|}
    \hline
    \multicolumn{5}{|c|}{Relative errors using different regularization methods, for $10\%$ data noise.}\\
    \hline
     & Constraint: & $x \geq 0$ & Constraint & $x \geq -10^{-100}$ \\
    \hline
    Methods & Relative error & Sparsity (\# 0s) & Relative error & Sparsity\\
    \hline
    IRell1 & 0.2285 & 0 & 0.2285 & 0 \\
    \hline
    IRhtv & 0.3462 & 23454 & 0.4236 & 0 \\
    \hline
    IRcgls & 0.1717 & 0 & 0.1717 & 0 \\
    \hline
    IRfista & 0.1767 & 13590 & 0.1767 & 0 \\
    \hline
    IRmrnsd & 0.1618 & 1 & 0.1618 & 1 \\
    \hline
    IRenrich & 0.1717 & 0 & 0.1717 & 0 \\
    \hline
    IRirn & 0.1756 & 15699 & 0.3840 & 0 \\
    \hline
    IRhy.\_fgmres & 0.2285 & 0 & 0.2285 & 0 \\
    \hline
    IRhy.\_gmres & 0.1918 & 0 & 0.1918 & 0 \\
    \hline
    IRhy.\_lsqr & 0.2539 & 0 & 0.2539 & 0 \\
    \hline

    {\color{red} Truncated-gradients:}  \\
    \hline
    $\xi_{m,\lambda} = \xi_m$, $\lambda = 0$, no truncation & 0.1700 & 12644 & 0.1700 & 0 \\
    \hline
    $\xi_{m,\lambda}$=$\xi_m\chi_{|\xi_m| > |\lambda|}$ $\lambda = 10^{-5}$ & 0.1695 & 11460 & 0.1695 & 3747\\
    \hline
    $\xi_{m,\lambda}$=$\xi_m\chi_{|\xi_m| > |\lambda|}$ $\lambda = 3*10^{-5}$ & 0.1673 & 27440 & 0.1673 & 25961 \\
    \hline
    $\lambda_m$=$\alpha\%||\xi_m||_\infty$, $\alpha =10$ & 0.1739 & 21765 & 0.1739 & 19343 \\
    \hline
    $\lambda_m$=$\alpha\%||\xi_m||_\infty$, $\alpha=40$ & 0.1927 & 39592 & 0.1927 & 39015\\
    \hline
    $k = 32768 (\approx 50\%*n)$ & 0.1709 & 23945 & 0.1709 & 14491 \\
    \hline
    $k = 6554 (\approx 10\%*n)$ & 0.2265 & 39393 & 0.2265 & 39128\\
    \hline
    $k\lambda_m^+$, k=32768, $\lambda_m$=$40\% ||\xi_m||_\infty$ & 0.1709 & 23945 & 0.1709 & 14491\\
    \hline
    $k\lambda_m^+$, k=6554, $\lambda_m$=$10\%||\xi_m||_\infty$ & 0.1739 & 21765 & 0.1739 & 19343\\
    \hline
    \end{tabular}
    \caption{Image deblurring, Examples \ref{example hubble telescope, xmin=0} and \ref{example hubble telescope, xmin<0}}
    \label{hubble telescope error comparison}
\end{table}

\section{Practical stopping criterion}\label{modified DP}
As described above, terminating the iterations at an ``\textit{appropriate instance}" is the key to the regularizing effect in a (semi-) iterative regularization method, i.e., the iteration index ($m$) serves as a regularizing parameter. Where an appropriate instance is usually achieved by some \textit{stopping criteria} for the descent or minimization process. Typically, when the error norm $ ||b - b_\delta||_2 \leq \delta$ is known, one terminates the iterations when the residue norm is less than the error norm ($||Ax_{M(\delta)}^\delta - b_\delta||_2 \leq \eta \delta$, for some $\eta > 1$), known as \textit{Morozov's discrepancy principle} (DP) \cite{Morozov_a}. However, in a practical situation the exact knowledge of $\delta$ is rarely known, only an estimate of the noise present in the data can be made. Hence, the estimation of the noise present in data is also noisy, i.e., instead of the error-norm $\delta$ we have $\delta_{\delta_1}$, such that $|\delta_{\delta_1} - \delta| \leq \delta_1$, where $\delta_1$ is the error in estimating the error-norm $\delta$. This is crucial, since (now) the noisy DP will result in a recovered solution ($x_{M(\delta_{\delta_1})}^{\delta_1}$) which may be either overly-smoothed (or under-fitted) if $\delta_1 > \delta$, or noisy (over-fitted) if $\delta_1 < \delta$. Even if the exact $\delta$ can be estimated, it's known that the DP condition can sometimes lead to (slightly) under-fit (or over-smoothed) solutions, especially for small $\delta$, see Example \ref{example modified DP}.

In this section we provide a modified discrepancy principle which is a generalization of the discrepancy principle and serves as an improvement, especially, in the absence of error-norm ($\delta$) or when a noisy approximate ($\delta_{\delta_1}$) of the original $\delta$ is known. Moreover, it also improves the interpretability of the stopping criterion. First, note that the DP condition is equivalent to $\frac{||Ax_M^\delta - b_\delta||_2}{||b||_2} \leq \eta \frac{\delta}{||b||_2}$, i.e., terminate the iterations when the \textit{relative error} in the residue is less than the \textit{relative error} in the data, instead of the \textit{residue norm} less than the \textit{error-norm}. The advantage of this transformation is that for $||b||_2 >> 1$ (which is usually the case) we have $\frac{\delta_1}{||b||_2} << 1$, and hence, $\frac{\delta}{||b||_2} \approx \frac{\delta_{\delta_1}}{||b||_2}$. Now, although, we don't have the exact data $b$ and hence $||b||_2$, we can always approximate $||b||_2$ with $||b_\delta||_2$ (norm of the noisy data), since the error-norm $\delta$ is usually small. The following describes the stopping criterion algorithm:
\begin{enumerate}
    \item First, get an estimate ($\delta_{\delta_1}$) of the exact noise-norm ($\delta$), and hence, an estimate of the noise relative-error $\left( \frac{\delta_{\delta_1}}{||b_\delta||_2}, \mbox{ where $b_\delta$ is the noisy data} \right)$.
    \item Construct some critical points, such as $\gamma_1 = \ceil*{\frac{\delta_{\delta_1}}{||b||_2}}\%$, where $\ceil{x}$ is the smallest integer greater than $x$, $\gamma_2 = \gamma_1 - 0.5\%$, $\gamma_3 = \gamma_2 - 0.5\%$ and so on. Typically, $\gamma_1$, $\gamma_2$ and $\gamma_3$ are enough to get a good recovery, but one can opt for more $\gamma_i$'s and/or with finer spacings.
    \item Perform the minimization for certain number of iterations and captures the recoveries at those critical points, i.e., recoveries ($x_{M(\gamma_i)}^\delta$) at $M(\gamma_i)$, corresponding to $\gamma_i$, such that $\frac{||Ax_{M(\gamma_i)}^\delta - b_\delta||_2}{||b_\delta||_2}\% \leq \eta \gamma_i$.
    \item Now, setting $x_{M(\gamma_1)}^\delta$ as the base recovery, one can compare (heuristically) the other recoveries with it and select appropriately, such as the most sparse recovery close to $x_{M(\gamma_1)}^\delta$ etc., see Example \ref{example modified DP}.
\end{enumerate}

\begin{example}\label{example modified DP}
In this example we compare the results obtained by using the modified DP (as discussed in \S \ref{modified DP}) over the simple discrepancy principle. We considered the Hubble telescope image (from Example \ref{example hubble telescope, xmin=0}) but with a smaller relative error in the added noise, i.e., $\frac{\delta}{||b||_2}\% = 3.4242\%$, for $\delta = 0.098591$. The corresponding $\gamma_i$ values are $\gamma_1 = 4\%$, $\gamma_2 = 3.5\%$, $\gamma_3 = 3\%$ and $\gamma_4 = 2.5\%$. We performed two sets of minimization: (1) Constrained minimization ($x_i \geq 0$) for truncated and not-truncated gradients, (2) Unconstrained minimization ($x_i > -\infty$), for 100 iterations each. One can notice, from Table \ref{table modified DP}, that using the modified DP one can recover a regularized sparse solution in the absence of the exact error-norm ($\delta$), even with better sparsity and/or accuracy than the simple DP-solution.
\end{example}

\begin{table}[ht]
    \centering
    \begin{tabular}{|p{3cm}||p{1cm}|p{1.5cm}|p{1.8cm}|p{1.8cm}|p{1.1cm}|}
    \hline
    \multicolumn{6}{|c|}{Discrepancy principle vs. Modified discrepancy principle.}\\
    \hline
    {\cl{red}Unconstrained, $x_i > -\infty$} & index (m) & $\gamma_i$\% & $\frac{||Ax_m^\delta - b_\delta||_2}{||b_\delta||_2}\%$ & $\frac{||\hat{x} - x_m^\delta||_2}{||\hat{x}||_2}\%$ & sparsity \\
    \hline
    $\lambda_m=0$ & 5 & 4 & 3.8453 & 15.9971 & 0 \\
     & 7 & 3.5 & 3.4889 & 15.3097\ & 0\\
    $||Ax_m^\delta - b_\delta||_2 \leq \delta \mapsto$ & 8 & 3.4242 & 3.3836 & 15.0990 & 0 \\
     & 16 & 3 & 2.9915 & 14.6378 & 0\\
     & 64 & 2.5 & 2.4998 & 18.1620 & 0\\
     \hline
    \hline
    $\lambda_m= 10\%||\xi_m||_\infty$ & 5 & 4 & 3.8490 & 16.1104 & 2280 \\
     & 7 & 3.5 & 3.4713 & 15.3842 & 1640\\
    $||Ax_m^\delta - b_\delta||_2 \leq \delta \mapsto$ & 8 & 3.4242 & 3.3693 & 15.2295 & 1539 \\
     & 15 & 3 & 2.9908 & 14.6846 & 958\\
     & 53 & 2.5 & 2.4951 & 18.3960 & 36\\
     \hline
    \hline
    $\lambda_m= 40\%||\xi_m||_\infty$ & 9 & 4 & 3.9741 & 17.9874 & 18697 \\
     & 14 & 3.5 & 3.4152 & 16.9031 & 16884\\
    $||Ax_m^\delta - b_\delta||_2 \leq \delta \mapsto$ & 14 & 3.4242 & 3.4152 & 16.9030 & 16884 \\
     & 43 & 3 & 2.9714 & 16.2881 & 9390\\
     \hline
     \hline
     {\cl{red}Constrained, $x_i \geq 0$} \\
     \hline
     $\lambda_m = 0$ & 5 & 4 & 3.8607 & 15.9490 & 5791 \\
     & 8 & 3.5 & 3.4160 & 14.8190 & 10601 \\
     $||Ax_m^\delta - b_\delta||_2 \leq \delta \mapsto$ & 8 & 3.4242 & 3.4160 & 14.8190 & 10601 \\
     & 26 & 3 & 2.9991 & 13.7672 & 18355\\
     \hline
     \hline
     $\lambda_m = 10\%||\xi_m||_\infty$ & 5 & 4 & 3.8767 & 15.9377 & 11887 \\
     & 8 & 3.5 & 3.4915 & 15.0107 & 10732 \\
     $||Ax_m^\delta - b_\delta||_2 \leq \delta \mapsto$ & 10 & 3.4242 & 3.3720 & 14.6331 & 11177 \\
     & 31 & 3 & 2.9990 & 13.7393 & 19756\\
     \hline
     \hline
     $\lambda_m = 40\%||\xi_m||_\infty$ & 10 & 4 & 3.8407 & 17.6140 & 31385 \\
     & 19 & 3.5 & 3.4318 & 16.8342 & 29986 \\
     $||Ax_m^\delta - b_\delta||_2 \leq \delta \mapsto$ & 21 & 3.4242 & 3.3981 & 16.6935 & 29877 \\
     \hline
    \end{tabular}
    \caption{Relative errors using the modified DP, Example \ref{example modified DP}.}
    \label{table modified DP}
\end{table}

\begin{remark}
Note that, in our iterative algorithm both the parameters: (1) the iteration index ($m$) and (ii) the truncation values ($\lambda, \lambda_m$ or $k$) serve as the regularization parameter as well as the sparsity parameter. That is, in a regularized setting ($m \leq M_0$, for some appropriate $M_0 \in \mathbb{N}$), the sparsity of the recovery decreases with increasing $m$ values (and/or decreasing truncation values) and vice-versa, where as, the accuracy (relative error) of the recovery increases with increasing $m$ values (and/or decreasing truncation values) and vice-versa.
\end{remark}

\section{Conclusion}
In this paper we presented a very simple iterative regularization method that motivates sparsity in the reconstruction. A sparse regularized solution of an inverse problem is achieved by opting a descent direction, during the minimization process, that motivates sparsity. We provided many such descent directions that serve the purpose. We prove the convergence of this algorithm to a regularized (sparse) solution when dealing with a noisy data. We also provided an efficient and effective (practical) stopping criterion to terminate the iterations during the descent process, when dealing with a noisy data, without the knowledge of the exact error-norm.

\bibliography{thesisref}
\bibliographystyle{ieeetr}

\end{document}


TITLE: Recurrence relation paths traced on a tetrahedron
QUESTION [2 upvotes]: I am attempting a recurrence relation question in which the number of possible length $n$ paths $a_n$ is being calculated. A path is defined as a traversal along a tetrahedron between connected vertices. All valid paths must return to the original vertex at the end of the path. The example given labels the vertices $ABCD$ and shows that the value of $a_2$ is $3$ - $ABA$, $ACA$ and $ADA$ are the only valid paths.
The question is what is the formula that defines this recurrence relation? The answer stated is $a_n=2a_{n-1}+3a_{n-2}$ but I don't understand where the $2a_{n-1}$ comes from. The $3a_{n-2}$ term makes sense as one could traverse any length two path - $3$ of them - and then do a path of length $n-2$ and this would be a valid path of length $n$.

REPLY [2 votes]: The $2 a_{n-1}$ term can be described by taking any path, which must start with $AX$ where $X \in \{B,C,D\}$. Now, which ever one $X$ is, there are two others $Y$ and $Z$, for example if $X=B$ then $Y=C$ and $Y=D$. So you can replace $AX$ with either $AYX$ or $AZX$, which gives $2 a_{n-1}$ new paths.
But there's still an issue. The reason why this method counts every length $n$ path exactly once is:

If a path starts $AXA$ then it is counted in the $3 a_{n-2}$ term --- just remove the initial $AXA$;
If it starts $AXY$ then it is counted in the $2 a_{n-1}$ term --- just replace the initial $AXY$ by $AY$.


TITLE: How do I calculate the surface area of a cone with a segment removed? (see picture)
QUESTION [1 upvotes]: Suppose I have a cone with a segment removed, as shown in the following image:
https://i.stack.imgur.com/lIziU.png
I can assume all the cone dimensions are known, and that the plane intersecting the cone cutting it into the two sections is parallel to the axis of the cone. Then, I'd like to calculate the surface area of each of the curved sections (I don't need the surface areas of the flat sections). However, I'm not sure how to go about it, as the radial symmetry has been broken.
Apologies if this has been asked before, but I couldn't find the answer anywhere. Any advice would be greatly appreciated, thank you!

REPLY [1 votes]: If we look at the cone from a point of view parallel to the base.  There is a linear relationship between the slant length and the base radius.
And when we look as some section of the wall of the cone, there is a corresponding projection on the base.  And the difference in the area is proportional to this $\frac {\text {slant length}}{\text {radius}} = \sec \theta$ ratio.
One factor of distance has been dilated in a linear way, so area has been dilated by the same factor.
If you know multivariate calculus I can show you that way, too.


TITLE: Branch cut and $\log(z)$ derivative
QUESTION [5 upvotes]: I'm a little confused about the branch cut thing. Given an entire functions $f(z),g(z),h(z)$, $z\in \mathbb C$, such that $f(x)=g(x)+h(x)$ for all $x\in \mathbb R$, $f$ and $g+h$ doesn't vanish on $\mathbb  R$ . I take the $\log$ for both sides then differentiate and get:
$$\log f(x)=\log(g(x)+h(x))$$
$$\frac{f'(x)}{f(x)}=\frac{g'(x)+h'(x)}{g(x)+h(x)}$$
I thought this is correct, but I was asked about "what is the branch cut I used here?", and I don't know what does this mean, and what is the answer for this question! Any help!?
Edit:
I also have the same problem with the following case:
If $f(x)=e^{g(x)}$ then $\log f(x)=g(x)$, and $\frac{f'(x)}{f(x)}=g'(x)$, also a branch cut issue!

REPLY [15 votes]: Sometimes, problems arise when trying to extend continuous functions from $\mathbb{R}$ to $\mathbb{C}$. For example, consider the complex-valued function $f(z) = \sqrt{z}$. Consider the following (poorly-drawn) diagram:

If your domain $D$ sweeps out a quarter-circle of radius $r$, the the image $f(D)$ sweeps out half the angle (an eighth-circle) with radius $\sqrt{r}$. Suppose then that your domain $D$ was an entire circle in the domain.

Consider the blue-red interface in the domain. As the image suggests, $\sqrt{z}$ is discontinuous there, because as you transition across the line $x = 0$ in the domain, you jump half-way across the circle in the image. To deal with this, we extend $\sqrt{z}$ to a new domain: two complex planes "glued" together. It looks something like this.

Now the function is continuous. Instead of transitioning directly from (1) to (2) in the image plane, we "move through the cut" to the blue (3) on the second complex plane, so that $\sqrt{z}$ is continuous. Then we travel around the second plane, jumping from the red (4) to the red (2) in the first complex. We have gone around twice in the domain, and once in the image, and our issues with continuity are resolved. To be more precise about the gluing, the domain actually looks like this:

Now, as you can expect, you need to glue three complex planes together to make $\sqrt[3]{z}$ continuous, and so on, etc. Unfortunately, for the complex logarithm, no number of complex planes in the domain will make the image continuous. Thus, you need to define the complex logarithm on an infinite collection of complex planes, "spiraling" away in either direction. It looks like this.

Recall the example with the square root earlier. We had two copies of the complex plane, and their images under $\sqrt{z}$ were different. In the first copy, the image was the right semi-circle, and in the second copy, it was the left semicircle. Thus, when working with the complex square root, one needs to specify which copy of the complex plane one wishes to work with. Each such copy is called a branch (symmetrically, you can think of it as two copies of $\sqrt{z}$ defined on the same complex plane). Similarly, when working with the complex log, you need to talk about which of the infinitely many complex planes in the domain you wish to work with, and so you must specify which branch you are using. The popular choice is the so-called "principal branch":
Log$(z) = ln|z| + i$Arg$(z)$.

REPLY [2 votes]: The problem here is that $\log(z)$ is a multi-valued function on the complex plane. How can you see this? First of all recall that $\log(z)$ is usually defined to be an inverse of $e^{z}$ (i.e. $log(e^{z})=z$). Working out what this means shows that one way of defining the $\log$ is to take
$$\log(z)=\ln|z| +iArg(z)$$
but then you run into the problem that $Arg(z)$ has an infinite number of values, just add or subtract a multiple of $2 \pi$. Now consider $\log(1)$. We want $\log$ to be continuous (and even analytic), so it should satisfy
$$\lim_{y \to 0^{+}} \log(1+iy) =\lim_{\theta \to 2\pi^{-}} \log(1e^{i\theta}) = 0$$
The first one appears to be true, but the above definition of log shows that in the second limit $Arg(z)$ approaches $2\pi$ so that
$$\lim_{\theta \to 2\pi^{-}} \log(1e^{i\theta}) = 2\pi i$$
The problem here is that the domain of the log function is not simply connected. A way to fix this problem is to remove all of the negative real axis and zero so that you can't go all the way around $0$ and force 2 different values of $Arg(z)$. Removing these points gives you an analytic function on $\mathbb{C}\setminus \{x\in \mathbb{R}|x\leq 0\}$ and this is called a branch cut because it singles out a single valued function of $\log$ from an infinite number of choices for the range of $Arg(z)$. One can also make a branch cut by removing any ray that connects zero and infinity.
In your problem, we are not sure of what domain your final statement is true on as we dont know what branch of the $\log$ function you have chosen.


TITLE: graphing and functions
QUESTION [1 upvotes]: one of my problem says that the average slope formula is: $\dfrac{f(x_2)-f(x_1)}{x_2-x_1}$
It tells us the average slope over the interval from $x_1$ to $x_2$
Then it says to take the average slope of the function $x^2+3x+1$ over one unit intervals from $0$ to $10$.
can i use a particular function function to solve this problem?

REPLY [0 votes]: The naive approach is to make yourself a table, with $10$ rows (one for each $x_i$ with: $x_1 = 1, x_2 = 2, \dots, x_{10} = 10$). The first column would be $f(x_i)$, the second column would be $x_i - x_{i-1}$ (you should notice an obvious pattern in this column), the third column would be $f(x_i) - f(x_{i-1})$, and the fourth column would be:
$\dfrac{f(x_i) - f(x_{i-1})}{x_i - x_{i-1}}$
(and if you noticed the pattern in the second column, you might not need this last one....why?)
To get this started, set $x_0 = 0$, with $f(x_0) = f(0) = 1$.
Note that Henry's answer above gives a short-cut, this is "the long way".
The first row (for $x_1 = 1$) would look like this:
$f(1) = 5;\ \ x_1 - x_0 = 1 -0 = 1;\ \ f(1) - f(0) = 5-1 = 4;$ slope $= 4/1 = 4$.
Nine more to go....


TITLE: Dry adiabatic lifting
QUESTION [3 upvotes]: If an unsaturated air parcel ascends adiabatically in the atmosphere with an initial temperature, would the final temperature of the air parcel, when it reaches it's destination, be the same at the initial temperature since it does not exchange heat with it's surroundings?

REPLY [3 votes]: No. The parcel of gas expands adiabatically as it rises. This makes it cool.
Here's one way of looking at what happens: No heat is transferred to or from the parcel of gas (so, adiabatic). At the same time, the parcel of gas expands as it rises. That means the parcel is doing work on the external environment. The temperature of the parcel must drop to conserve energy.


TITLE: If a chain of distinct irreducible closed subsets of a quasi-affine variety $Y$ is maximal, a chain of their closures is maximal in $\overline{Y}$?
QUESTION [5 upvotes]: The following is the Proposition 1.10 of Hartshorne's Algebraic Geometry:

If $Y$ is a quasi-affine variety, then $\dim Y =\dim \overline{Y}$.

In the proof, there is a statement saying that if
$$
Z_0 \subset \cdots \subset Z_n \;\;\; (1)
$$
is a maximal chain of distinct irreducible closed subsets of $Y$, then
$$
\overline{Z}_0 \subset \cdots \subset \overline{Z}_n   \;\;\; (2)
$$
is also a maximal chain of distinct irreducible closed subsets of $\overline{Y}$
refering to the example (1.1.3) which is

Any nonempty open subset of an irreducible space is irreducible and dense.

I can't understand this. Yes, Y is an open subset of an affine variety. Then, what ?
Help me!

REPLY [2 votes]: We are given a maximal chain in $Y$ given by $(1)$. Each $Z_i$ is irreducible in $Y$ and hence is irreducible in $X=\bar{Y}$. We know that if $Z_i$ is irreducible in $X$ then it's closure is also irreducible (and of course $\bar{Z_i}$ is closed in $X$). What remains to show is that your chain in $(2)$ is maximal. Suppose we could extend it to $\bar{Z_n}\subset W$ with proper containment and $W$ is closed and irreducible.
Let $W' = Y\cap W$. Since $W$ is closed, $W'$ is closed in $Y$. Now we use your quoted theorem. $W'$ is irreducible since it is a non-empty open subset of an irreducible set.
Finally, $Z_n\subset W'$ since intersection respects containment, and is strict since their closures are distinct (We are using the fact that $W'$ is dense in $W$ so that $\bar{W'}=W$, which again makes use of your quoted theorem). Thus we have extended the maximal chain $(1)$ which is a contradiction.


TITLE: Given $a+b+c=abc$, find $\tan^{-1}a+\tan^{-1}b+\tan^{-1}c$.
QUESTION [0 upvotes]: I need help with the following problem - a hint of how to approach it would be fine:

Given $a$, $b$, and $c$ are real numbers and $a+b+c=abc$, find
  $$\tan^{-1}(a)+\tan^{-1}(b)+\tan^{-1}(c).$$

REPLY [1 votes]: By the trigonometric identity
$$
\tan(\alpha+\beta+\gamma) = \frac{\tan\alpha+\tan\beta+\tan\gamma-\tan\alpha\tan\beta\tan\gamma}{1-\tan\alpha\tan\beta-\tan\alpha\tan\gamma-\tan\beta\tan\gamma}.
$$
let $\alpha$, $\beta$, $\gamma$ be respectively the arctangents of $a$, $b$, and $c$ so that
$$
\frac{a+b+c-abc}{1-ab-ac-bc} = \tan\Big(\arctan(a)+\arctan(b)+\arctan(c)\Big)
$$
then since $a+b+c=abc$ we have that the numerator of the above equation is zero
$$\tan\Big(\arctan(a)+\arctan(b)+\arctan(c)\Big)=0$$
from which we know the tangent function is zero at integer multiples of $\pi$. However, we shouldn't accept all integer values of $\pi$ since $\arctan(\theta):\mathbb R \to (-\frac{\pi}{2},\frac{\pi}{2})$ and therefore the range of $\arctan(\theta)$ is restricted to $(-\frac{\pi}{2},\frac{\pi}{2})$. Following the technique shown inside this question, we may conclude
$$\arctan(a)+\arctan(b)+\arctan(c) = \begin{cases}
\pi, & 1 < ab+bc+ca ~~\text{and}~~ a+b+c > 0\\
0,   & 1 > ab+bc+ca\\
-\pi & 1 < ab+bc+ca ~~\text{and}~~ a+b+c < 0
\end{cases}
$$


TITLE: Show that simple bipartite graph with n vertices doesn't have more than $\frac{n^2}4$ edges.
QUESTION [0 upvotes]: I am stuck at finding the solution because the lack understanding of bipartite graph theory.
Bipartite graph is defined as a graph which can be separated into two sets of vertices such as that each vertices in one set has edges connecting the vertices in the another set.
Please, could you show me an answer for the question.

REPLY [3 votes]: In a bipartite graph your $n$ vertices can be partitioned into two subsets of size $i$ and $(n-i)$, $0 \leq i \leq n$, and you can only have edges between vertices of different subsets, so you have a maximum of $i(n-i)$ edges if every member of one subset is connected to every member of the other subset.
$f(i) = i(n-i), 0 \leq i \leq n$ is maximized by $i = \frac{n}{2}$ which leads to $\frac{n}{2} \frac{n}{2} = \frac{n^2}{4}$ being the maximum number of edges.

REPLY [1 votes]: Hint:
A bipartite graph is divided into two pieces, say of size $p$ and $q$, where $p+q=n$. Then the maximum number of edges is $pq$.
Using calculus, we can deduce that this product is maximal when $p=q$, in which case it is equal to $n^2/4$.

REPLY [0 votes]: Suppose a graph $G$ is partitioned into two disjoint subsets $U$ and $V$ then $G$ is bipartite if every edge connects a vertex in $U$ to a vertex in $V$. Therefore, for each vertex in $U$ there are at most $|V|$ many possible edges origininating from it . So we have $|U||V|$ possible edges in the whole graph at most. Now with the constraint  that $|U|+|V|=n$ we maximize $|U||V|$ and this is indeed the case when $|U|=|V|=\frac{n}{2}$ that is we have at most $|U||V|=\frac{n^2}{4}$ edges.


TITLE: Kalman Filter...Denoising measurement data to track objects
QUESTION [1 upvotes]: Hi Everyone,
I am about to implement a Kalman Filter in a software.
I found this very helpful article here:
http://bilgin.esme.org/BitsBytes/KalmanFilterforDummies.aspx
The example helps a lot, but in my case I have measurement data and that data tracks an object (if available) so I do not expect the same measurement data each time I take a measurement.
So my measurement data is not constant, it changes depending if a trackable object appears.
In the article above we have the constant called A, and in my case I guess it should be calculated somehow... But I am not sure how...
Any help would be really great!
Thanks a lot!

REPLY [2 votes]: This book describes such models at extreme length and with almost all possible variations.

West and Harrison (1997). Bayesian Forecasting and Dynamic Models

Edit: This an excellent additional reference that is more empirical in nature.

Petris, Petrone and Campagnoli (2009). Dynamic Linear Models with R


TITLE: Elliptic functions (in particular, the Weierstrass P-function) to solve ODE
QUESTION [1 upvotes]: I have the following ODE
$$-\frac{vf^2}{2} = \frac{(f')^2}{2} + 2f^3 + Cf + D,$$
where $v, C, D$ are constants.
I've been given the hint to use elliptic functions (in particular, the Weierstrass P-function). I've looked at the wikipedia page for it, but really can't see how it can be of use.
Any help/hints?

REPLY [3 votes]: HINT.-Clearly first you should know well the $\wp_L(z)$ for the lattice $L$. It is defined by
$$\wp(z)=\frac{1}{z^2}+\sum_{\omega\in L - \{0\}}\big [\frac{1}{(z-\omega)^2}-\frac{1}{\omega^2}\big ]$$ and an important algebraic relation between the trascendental funtions $\wp$ and $\wp'$ is $$(\wp')^2=4(\wp)^3-g_2\wp-g_3\qquad (*)$$ where the called invariants of $\wp$, $g_2$ and $g_3$ are defined by
$$g_2=60\sum_{\omega\in L - \{0\}}\frac{1}{\omega^4}\\g_3=140\sum_{\omega\in L - \{0\}}\frac{1}{\omega^6}\\$$ Now your equation has the form
$$(f')^2=4(-f)^3-vf^2-2Cf-2D\qquad (**)$$
Try to associate the relations $(*)$ and $(**)$ in such a way that $f(x)=\wp'(t)$ and $x=\wp(t)$. You have to construct a lattice with $v,C,D$ and satisfying $C^3-27D^2\ne 0$. You can discard the term $vf^2$ in $(**)$ like as usual for solving the cubic equation (for a beginner all this it is not easy at all).