∎

11institutetext: Zhong-Feng Sun 22institutetext: School of Mathematics and
Statistics, Shandong University of Technology, Zibo, Shandong, China. 22email:
<EMAIL_ADDRESS>33institutetext: Jin-Chuan Zhou 44institutetext: School of
Mathematics and Statistics, Shandong University of Technology, Zibo, Shandong,
China. 44email<EMAIL_ADDRESS>55institutetext: Yun-Bin Zhao
66institutetext: Corresponding author. Shenzhen Research Institute of Big
Data, Chinese University of Hong Kong, Shenzhen, Guangdong, China. 66email:
<EMAIL_ADDRESS>
# Heavy-ball-based optimal thresholding algorithms for sparse linear inverse
problems ††thanks: The work was founded by the National Natural Science
Foundation of China (NSFC 12071307 and 11771255), Young Innovation Teams of
Shandong Province (2019KJI013), and Shandong Province Natural Science
Foundation (ZR2021MA066).

Zhong-Feng Sun  Jin-Chuan Zhou  Yun-Bin Zhao

(Received: date / Accepted: date)

###### Abstract

Linear inverse problems arise in diverse engineering fields especially in
signal and image reconstruction. The development of computational methods for
linear inverse problems with sparsity tool is one of the recent trends in this
area. The so-called optimal $k$-thresholding is a newly introduced method for
sparse optimization and linear inverse problems. Compared to other sparsity-
aware algorithms, the advantage of optimal $k$-thresholding method lies in
that it performs thresholding and error metric reduction simultaneously and
thus works stably and robustly for solving medium-sized linear inverse
problems. However, the runtime of this method remains high when the problem
size is relatively large. The purpose of this paper is to propose an
acceleration strategy for this method. Specifically, we propose a heavy-ball-
based optimal $k$-thresholding (HBOT) algorithm and its relaxed variants for
sparse linear inverse problems. The convergence of these algorithms is shown
under the restricted isometry property. In addition, the numerical performance
of the heavy-ball-based relaxed optimal $k$-thresholding pursuit (HBROTP) has
been studied, and simulations indicate that HBROTP admits robust capability
for signal and image reconstruction even in noisy environments.

###### Keywords:

Sparse linear inverse problems optimal $k$-thresholding heavy-ball method
restricted isometry property phase transition peak signal-to-noise ratio

###### MSC:

94A12 15A29 90C25 90C20 49M20

##  1 Introduction

In recent years, the linear inverse problem has gained much attention in
various fields such as wireless communication BSR17 ; CZJAZ21 , signal and
image processing OS18 ; BT09 ; TR20 ; BD10 ; XDY21 ; LAZS20 . The main task of
the linear inverse problem is to recover the unknown data $z\in\mathbb{R}^{n}$
from the linear measurements

$y=\Phi z+\nu,$ (1.1)

where $y\in\mathbb{R}^{m}$ are the acquired measurements of $z$,
$\nu\in\mathbb{R}^{m}$ is a noisy vector, and $\Phi\in\mathbb{R}^{m\times n}$
is a given measurement matrix. In this paper, we only consider the
underdetermined case $m\ll n.$ In this case, it is generally impossible to
reconstruct data $z$ from the linear system (1.1) unless $z$ possesses a
certain special structure such as sparsity. When $z$ in its current domain is
not sparse enough, a practical approach to achieve sparsity is to consider a
certain transformed domain in which $z$ can be sparsely represented. For
instance, many natural signal and image can be sparsely represented under the
fourier transform or some wavelet transform. Suppose that $z$ can be sparsely
represented on an orthonormal basis $\Psi\in\mathbb{R}^{n\times n}$, i.e.,
$z=\Psi x$ where the coefficient vector $x\in\mathbb{R}^{n}$ is $k$-sparse or
$k$-compressible. That is, $\|x\|_{0}\leq k$ or $x$ can be well approximated
by a $k$-sparse vector, where $\|\cdot\|_{0}$ denotes the number of nonzero
entries of a vector. Via such a sparse representation, the non-sparse linear
inverse model (1.1) can be transformed to the sparse linear inverse (SLI)
problem

$y=Ax+\nu,$ (1.2)

where the product $A=\Phi\Psi\in\mathbb{R}^{m\times n}$ is still referred to
as the measurement (or sensing) matrix. The purpose of the SLI problem is to
reconstruct the sparse signal $x$ based on the observed (inaccurate) signal
$y$ and the adopted measurement matrix $A$ and then reconstruct the original
signal $z=\Psi x.$ The SLI problem can be formulated as the sparse
optimization problem BD10 ; CZJAZ21 ; OS18 ; TW10

$\underset{u}{\rm min}\\{{\|y-Au\|}_{2}^{2}:\|u\|_{0}\leq k\\}$ (1.3)

or other related models such as the LASSO-type models BT09 ; BSR17 . The
problem (1.3) has been widely used in signal and image processing via
compressed samplings.

The existing methods for (1.3) includes thresholding E10 ; E12 ; FR13 , greedy
DM09 ; NT09 ; TG07 , convex optimization CT05 ; CWB08 ; CDS98 ; ZL17 ,
nonconvex optimization C07  and Bayesian methods SPZ08 ; WB04 . In this paper,
we focus our attention on the thresholding method which was first proposed by
Donoho and Johnstone DJ94 . It has experienced a significant development since
1994 and has evolved into a large family of algorithms including the hard
thresholding BTW15 ; BD08 ; BD10 ; F11 ; KC14 ; MZ20 , soft thresholding DDM04
; D95 ; E06  and optimal $k$-thresholding MZ22 ; MZKS22 ; Z20 ; ZL21 . It is
worth stressing that the main advantage of thresholding methods is that the
algorithms can always guarantee the points generated by the algorithms being
feasible or nearly feasible to the problem (1.3), since the iterates are
$k$-sparse or $k$-compressible. The simplest thresholding iterative method
might be the so-called iterative hard thresholding BD08  which is a gradient-
based thresholding method. The combination of IHT and orthogonal projection
yields the well known hard thresholding pursuit (HTP) method F11 . Due to
their low computational complexity, IHT, HTP and their modifications have been
widely used in numerous areas such as compressive sampling and image
processing BTW15 ; B12 ; BD10 ; KC14 . However, performing hard thresholding
on an iterate does not necessarily leads to the reduction of the objective
function value of (1.3) and might also cause numerical oscillation during the
course of iterations Z20 ; ZL21 .

To alleviate the above-mentioned weakness of hard thresholding, the optimal
$k$-thresholding (OT) operator was first introduced in Z20  (see also ZL21 ).
Such an operator may reduce the objective value of (1.3) while performing
thresholding of an iterate. Based on this operator, the basic framework of
optimal $k$-thresholding (OT) and optimal $k$-thresholding pursuit (OTP) is
proposed in Z20 . The optimal $k$-thresholding of a vector $v$ is defined as

${\rm min}\\{{\|y-A(v\otimes{w})\|}_{2}^{2}:~{}{\bf
e}^{T}w=k,~{}w\in\\{0,1\\}^{n}\\},$ (1.4)

where ${\bf e}=(1,1,\ldots,1)^{T}\in\mathbb{R}^{n}$, $v\in\mathbb{R}^{n}$ is a
given vector, $\\{0,1\\}^{n}$ is the set of $n$-dimensional binary vectors and
$v\otimes w:=(v_{1}w_{1},\ldots,v_{n}w_{n})^{T}$ denotes the Hadamard product
of two vectors. From a computational point of view, it is generally more
convenient to solve the convex relaxation of (1.4) instead, which is the
following quadratic convex optimization

${\rm min}\\{{\|y-A(v\otimes{w})\|}_{2}^{2}:\ \ {\bf e}^{T}w=k,\ \ 0\leq
w\leq\bf e\\}.$ (1.5)

This problem is referred to as data compressing problem in Z20 ; ZL21 . Based
on (1.5), the relaxed optimal $k$-thresholding (ROT$\omega$) and relaxed
optimal $k$-thresholding pursuit (ROTP$\omega$) algorithms were proposed in
Z20 ; ZL21 , where $\omega$ represents times of data compression that are
performed in the algorithms. When $\omega=1$, the algorithm is named as ROTP.
Using partial gradient and Newton-type search direction, some modifications of
ROTP were also studied in MZ22 ; MZKS22  recently. However, solving the above-
mentioned quadratic optimization problem (1.5) remains time-consuming for
large-scale problems. It is important to further investigate how the
computational cost of ROTP-type methods might be reduced or how these methods
can be accelerated by incorporating such acceleration techniques as heavy-ball
(HB) and Nesterov’s method. By using linearization together with certain
regularization method, the so-called nature thresholding algorithm is newly
developed in ZL22  which has significantly lower the computational complexity
than ROTP${\omega}$. In this paper, we investigate the ROTP-type algorithm
from acceleration perspective by showing that the HB technique does be able to
improve the efficiency of the ROTP-type algorithms.

The HB method introduced by Polyak P64  can be seen as a two-step method which
combines the momentum term and gradient descent direction. In recent years, HB
has gained significant attention and has found wide applications in image
processing, data analysis, distributed optimization and undirected networks
ADR22 ; GOP17 ; KBGY22 ; LRP16 ; LCWW21 ; MRJ21 ; UPS22 ; XK20 . The
theoretical analysis (global convergence and local convergence rate) has been
studied by several researchers. For example, the linear convergence rate of HB
for unconstrained convex optimization problem was established by Aujol et. al
ADR22 ; Mohammadi et. al MRJ21  analyzed the relation between the convergence
rate of HB and its variance amplification when the objective functions of
problems are strongly convex and quadratic; Xin and Khan XK20  showed that the
distributed HB method with appropriate parameters attains a global $R$-linear
rate, and it has potential acceleration compared with some first-order methods
for the ill-conditioned problem. Besides the HB acceleration technique, the
Nesterov’s acceleration method are also studied in KBGY22 ; LRP16 ; LCWW21 ;
MRJ21 .

In this paper, we study the acceleration issue of the optimal $k$-thresholding
algorithms by incorporating the HB acceleration technique. Merging the optimal
$k$-thresholding and heavy-ball techniques, we propose the following novel
algorithms for the linear inverse problems reformulated as (1.3):

  * •

Heavy-ball-based optimal $k$-thresholding (HBOT),

  * •

Heavy-ball-based optimal $k$-thresholding pursuit (HBOTP),

  * •

Heavy-ball-based relaxed optimal $k$-thresholding ($\textrm{HBROT}\omega$),

  * •

Heavy-ball-based relaxed optimal $k$-thresholding pursuit
($\textrm{HBROTP}\omega$),

where the integer parameter $\omega$ denotes the number of times for solving
(1.5) at every iteration. The global convergence of these algorithms is
established in this paper under the restricted isometry property (RIP)
introduced by Candès and Tao CT05 , and the main results are summarized in
Theorems 3.1 and 4.1. The performance of HBROTP (i.e., $\textrm{HBROTP}\omega$
with $\omega=1$) and several existing algorithms such as ROTP2 Z20 , PGROTP
MZKS22 , $\ell_{1}$-minimization CDS98  and OMP E10  are compared through
numerical experiments. The phase transition with Gaussian random data is
adopted to illustrate the performance of the proposed algorithms for solving
sparse signal reconstruction problems. Also, a few real images are used to
demonstrate the reconstruction performance of the algorithms in noiseless
situations.

This paper is organized as follows. In Section 2, we introduce some notations,
definitions, basic inequalities and the algorithms. In Section 3, we discuss
the error bounds and convergences of HBOT and HBOTP under the RIP. Similarly,
the error bounds for HBROT$\omega$ and HBROTP$\omega$ are given in Section 4.
Numerical results observed from synthetic signals and real images are
estimated in Section 5. Conclusions are drawn in Section 6.

##  2 Preliminary and algorithms

###  2.1 Notations

Denote by $N:=\\{1,2,\ldots,n\\}.$ Given a subset $\Omega\subseteq N,$
$\overline{\Omega}:=N\setminus\Omega$ denotes the complement set of $\Omega$
and $|\Omega|$ denotes its cardinality. For a vector $z\in\mathbb{R}^{n}$, the
support of $z$ is represented as $\textrm{supp}(z):=\\{i\in N:z_{i}\neq 0\\}$,
and the vector $z_{\Omega}\in\mathbb{R}^{n}$ is obtained by zeroing out the
elements of $z$ supported on $\overline{\Omega}$ and retaining other elements.
Given the sparse level $k$, $\mathcal{L}_{k}(z)$ denotes the index set of the
$k$ largest absolute entries of $z.$ As usual,
$\mathcal{H}_{k}(z):=z_{\mathcal{L}_{k}(z)}$ is called the hard thresholding
of $z.$ The symbols $\|\cdot\|_{1}$ and $\|\cdot\|_{2}$ represent
$\ell_{1}$-norm and $\ell_{2}$-norm of a vector, respectively. Let
$\mathcal{W}^{k}:=\\{w:~{}{\bf e}^{T}w=k,~{}w\in\\{0,1\\}^{n}\\}$ and
$\mathcal{P}^{k}:=\\{w\in\mathbb{R}^{n}:{\bf e}^{T}w=k,\ \ 0\leq w\leq\bf
e\\}$ be two subsets of $\mathbb{R}^{n}$.

###  2.2 Definitions and basic inequalities

Let us first recall the restricted isometry property (RIP) of a matrix $A$ and
the optimal $k$-thresholding operator $Z^{\\#}_{k}(\cdot)$ in the following
two definitions respectively.

###### Definition 1

CT05  Given a matrix $A\in\mathbb{R}^{m\times n}$ with $m<n$. The $k$th order
restricted isometry constant (RIC) of $A$, denoted by $\delta_{k},$ is the
smallest nonnegative number $\delta$ such that

$(1-\delta){\left\|u\right\|}^{2}_{2}\leq{\left\|Au\right\|}^{2}_{2}\leq(1+\delta){\left\|u\right\|}^{2}_{2}$
(2.1)

for all $k$-sparse vectors $u\in\mathbb{R}^{n}$. The matrix $A$ is said to
satisfy the RIP of order $k$ if $\delta_{k}<1$.

###### Definition 2

Z20 ; ZL21  Given a vector $u\in\mathbb{R}^{n}$, let $w^{*}(u)$ be the
solution of the binary optimization problem

$\underset{w}{\min}\left\\{\|y-A(u\otimes
w)\|^{2}_{2}:~{}w\in\mathcal{W}^{k}\right\\}.$

The optimal $k$-thresholding of $u$ is denoted by $Z^{\\#}_{k}(u):=u\otimes
w^{*}(u)$ which is a $k$-sparse vector. $Z^{\\#}_{k}(\cdot)$ is called the
optimal $k$-thresholding operator.

The two lemmas below will be used for the analysis in Sections 3 and 4.

###### Lemma 1

_F11_ Let $u,v\in\mathbb{R}^{n}$, $W\subseteq N$ be an index set and $s,t\in
N$ be two integers.

  * (i)

If $|W\cup\textrm{supp}(v)|\leq t$, then
$\|\left[(I-A^{T}A)v\right]_{W}\|_{2}\leq\delta_{t}{\|v\|}_{2}.$

  * (ii)

If $|W|\leq t$, then
$\|\left(A^{T}v\right)_{W}\|_{2}\leq\sqrt{1+\delta_{t}}{\|v\|}_{2}.$

  * (iii)

If $\textrm{supp}(u)\cap\textrm{supp}(v)=\emptyset,$ $|\textrm{supp}(u)|\leq
s$ and $|\textrm{supp}(v)|\leq t$, then
$|u^{T}A^{T}Av|\leq\delta_{s+t}{\|u\|}_{2}{\|v\|}_{2}.$

###### Lemma 2

Let $\\{a^{p}\\}\subseteq\mathbb{R}~{}(p=0,1,\dots)$ be a nonnegative sequence
satisfying

$a^{p+1}\leq b_{1}a^{p}+b_{2}a^{p-1}+b_{3}$

for $p\geq 1,$ where $b_{1},b_{2}$ and $b_{3}\geq 0$ are constants and
$b_{1}+b_{2}<1$. Then

$a^{p}\leq\theta^{p-1}\left(a^{1}+(\theta-
b_{1})a^{0})\right)+\frac{b_{3}}{1-\theta}$

for $p\geq 2,$ where $0\leq\theta<1$ is a constant given by
$\theta=(b_{1}+\sqrt{b_{1}^{2}+4b_{2}})/2<1.$

The proof of Lemma 2 is straightforward and hence omitted.

###  2.3 Algorithms

Given iterates $x^{p-1}$ and $x^{p}$, the heavy-ball search direction is
defined as

$d^{p}=\alpha A^{T}(y-Ax^{p})+\beta(x^{p}-x^{p-1}),$

where $\alpha>0$ and $\beta\geq 0$ are two parameters. We use the optimal
$k$-thresholding operator $Z^{\\#}_{k}(\cdot)$ to generate the new iterate
$x^{p+1}$, i.e.,

$x^{p+1}=Z^{\\#}_{k}(x^{p}+d^{p}),$

which is called the heavy-ball-based optimal $k$-thresholding (HBOT) of the
iterate $x^{p}+d^{p}.$ Combining HBOT with the orthogonal projection
technique, we propose the heavy-ball-based optimal $k$-thresholding pursuit
(HBOTP) algorithm. The two algorithms are described as follows. It is evident
that HBOT and HBOTP are the multi-step extension of existing OT and OTP
algorithms in Z20 ; ZL21 , respectively.

HBOT and HBOTP algorithms.  Input the data $(A,y,k)$ and two initial points
$x^{0}$ and $x^{1}$. Choose the parameters $\alpha>0$ and $\beta\geq 0$.

  * S1.

At $x^{p}$, set

$u^{p}=x^{p}-\alpha A^{T}(Ax^{p}-y)+\beta(x^{p}-x^{p-1}).$ (2.2)

  * S2.

Solve the optimization problem

$w^{*}=\arg\min_{w}\\{{\|y-A(u^{p}\otimes{w})\|}_{2}^{2}:~{}{\bf
e}^{T}w=k,~{}w\in\\{0,1\\}^{n}\\}.$ (2.3)

  * S3.

Generate the next iterate $x^{p+1}$ as follows:

    * For HBOT, let $x^{p+1}=u^{p}\otimes{w}^{*}$.

    * For HBOTP, let $S^{p+1}=\textrm{supp}(u^{p}\otimes{w}^{*})$ and $x^{p+1}$ be the solution to the least squares problem

$x^{p+1}=\arg\min_{x\in\mathbb{R}^{n}}\\{{\|y-Ax\|}_{2}^{2}:\textrm{supp}(x)\subseteq
S^{p+1}\\}.$ (2.4)

Repeat S1-S3 until a certain stopping criterion is met.

In general, the computational cost for solving the binary optimization problem
(2.3) is usually high BT18 ; CAP08 . Replacing (2.3) by its convex relaxation

$\arg\min_{w}\left\\{\|y-A(u\otimes
w)\|^{2}_{2}:~{}w\in\mathcal{P}^{k}\right\\}$

yields the next heavy-ball-based relaxed optimal $k$-thresholding
(HBROT$\omega$) and the heavy-ball-based relaxed optimal $k$-thresholding
pursuit (HBROTP$\omega$) algorithms, where the parameter $\omega$ represents
the times of solving such a convex relaxation problem at each iteration
(which, as pointed out in Z20 , can be interpreted as the times of data
compression within each iteration). As $\omega=1$, we simply use HBROT and
HBROTP to denote the algorithms HBROT$1$ and HBROTP$1$, respectively. Clearly,
as $\alpha=1$ and $\beta=0$, HBROT$\omega$ and HBROTP$\omega$ reduce to
ROT$\omega$ and ROTP$\omega$ in ZL21 .

$\textrm{HBROT}\omega$ and $\textrm{HBROTP}\omega$ algorithms. Input the data
$(A,y,k)$, two initial points $x^{0},x^{1}$ and $\omega$. Choose the
parameters $\alpha>0$ and $\beta\geq 0$.

  * S1.

At $x^{p}$, calculate $u^{p}$ by (2.2).

  * S2.

Set $v\leftarrow u^{p}$. Perform the following loops to produce the vectors
${w}^{(j)}(j=1,\ldots,\omega)$:

for $j=1:\omega$ do

${w}^{(j)}=\arg\min_{w}\\{{\|y-A(v\otimes{w})\|}_{2}^{2}:\ \ {\bf e}^{T}w=k,\
\ 0\leq w\leq\bf e\\},$ (2.5)

and set $v\leftarrow v\otimes{w}^{(j)}$.

end

  * S3.

Let
$x^{\sharp}=\mathcal{H}_{k}(u^{p}\otimes{w}^{(1)}\otimes\cdots\otimes{w}^{(\omega)})$.
Generate the next iterate $x^{p+1}$ as follows:

    * For $\textrm{HBROT}\omega$, let $x^{p+1}=x^{\sharp}$.

    * For $\textrm{HBROTP}\omega$, let $S^{p+1}=\textrm{supp}(x^{\sharp})$, and $x^{p+1}$ be the solution to the least squares problem

$x^{p+1}=\arg\min_{x\in\mathbb{R}^{n}}\\{{\|y-Ax\|}_{2}^{2}:\textrm{supp}(x)\subseteq
S^{p+1}\\}.$ (2.6)

Repeat S1-S3 until a certain stopping criterion is met.

The choice of stopping criterions depends on the application scenarios. For
instance, one can simply prescribe the maximum number of iterations,
$p_{\max},$ which allows the algorithm to perform a total number of $p_{\max}$
iterations. One can also terminate the algorithm when
$\|y-Ax^{p}\|_{2}\leq\varepsilon,$ where $\varepsilon>0$ is a prescribed
tolerance depending on the noise level.

##  3 Analysis of HBOT and HBOTP

In this section, we establish the error bounds for HBOT and HBOTP under the
RIP of order $k$ or $k+1$. Taking into account of noise influence, the error
bound provides an estimation for the distance between the problem solution and
the iterates generated by the algorithms. Thus the error bound is an important
measurement of the quality of iterates as the approximation to the sparse
solution of the linear inverse problem. In noiseless situations, the error
bound implies the global convergence of the algorithms under the RIP
assumption. Let us first introduce the following property, which is a
combination of Lemmas 3.3 and 3.6 in ZL21 .

###### Lemma 3

_ZL21_ Let $z$ be a $(2k)$-sparse vector. Then
$\|Az\|_{2}^{2}\geq(1-2\delta_{k}-\delta_{k+s(k)})\|z\|_{2}^{2}$ with

$\displaystyle s(k)=\left\\{\begin{array}[]{ll}1,{}&\textrm{ if
}~{}k~{}\textrm{is an odd number},\\\ 0,{}&\textrm{ if }~{}k~{}\textrm{is an
even number}.\\\ \end{array}\right.$ (3.3)

Note that Lemma 3.4 in ZL21  was established for the sparsity level $k$ being
an even number. We now establish the similar result even when $k$ is odd.

###### Lemma 4

Let $h,z\in\mathbb{R}^{n}$ be two $k$-sparse vectors, and let
$\widehat{w}\in\mathcal{W}^{k}$ be any $k$-sparse binary vector satisfied
$\textrm{supp}(h)\subseteq\textrm{supp}(\widehat{w})$, then

$\|\left[(I-A^{T}A)(h-z)\right]\otimes\widehat{w}\|_{2}\leq\sqrt{5}\delta_{k+s(k)}{\|h-z\|}_{2},$
(3.4)

where $s(k)$ is given by (3.3).

###### Proof

For given vectors $h,z,\widehat{w}$ satisfying the conditions of the lemma,
from (ZL21, , Lemma 3.4), we get

$\|\left[(I-A^{T}A)(h-z)\right]\otimes\widehat{w}\|_{2}\leq\sqrt{5}\delta_{k}{\|h-z\|}_{2}$
(3.5)

for even number $k$. Therefore, we just need to show that (3.4) also holds
when $k$ is an odd number.

Indeed, assume that $k$ is an odd number. Taking a $(k+1)$-sparse binary
vector $\overline{w}\in\mathcal{W}^{k+1}$ such that ${\rm
supp}(\widehat{w})\subseteq{\rm supp}(\overline{w})$, we obtain

$\displaystyle\|\left[(I-A^{T}A)(h-z)\right]\otimes\widehat{w}\|_{2}$
$\displaystyle=\|\left[(I-A^{T}A)(h-z)\right]_{\textrm{supp}(\widehat{w})}\|_{2}$
$\displaystyle\leq\|\left[(I-A^{T}A)(h-z)\right]_{\textrm{supp}(\overline{w})}\|_{2}$
$\displaystyle=\|\left[(I-A^{T}A)(h-z)\right]\otimes\overline{w}\|_{2}.$ (3.6)

As $h$ and $z$ are two $k$-sparse vectors, they are also $(k+1)$-sparse
vectors. Since
$\textrm{supp}(h)\subseteq\textrm{supp}(\widehat{w})\subseteq\textrm{supp}(\overline{w})$,
and since $k+1$ is even (when $k$ is odd), applying (3.5) to this case yields

$\|\left[(I-A^{T}A)(h-z)\right]\otimes\overline{w}\|_{2}\leq\sqrt{5}\delta_{k+1}{\|h-z\|}_{2}.$
(3.7)

Combining (3) and (3.7), we obtain

$\|\left[(I-A^{T}A)(h-z)\right]\otimes\widehat{w}\|_{2}\leq\sqrt{5}\delta_{k+1}{\|h-z\|}_{2}$

for odd number $k$. Hence, we conclude that (3.4) holds for any positive
integers $k$. $\Box$

The main results for HBOT and HBOTP are summarized as follows.

###### Theorem 3.1

Let $x\in\mathbb{R}^{n}$ be a solution to the system $y=Ax+\nu$ where $\nu$ is
a noise vector. Assume that the RIC, $\delta_{k+s(k)},$ of $A$ and the
parameters $\alpha,\beta$ in HBOT and HBOTP satisfy that
$\delta_{k+s(k)}<\gamma^{*}$ and

$0\leq\beta<\frac{1+1/\eta}{1+\sqrt{5}\delta_{k+s(k)}}-1,~{}\frac{1+2\beta-1/\eta}{1-\sqrt{5}\delta_{k+s(k)}}<\alpha<\frac{1+1/\eta}{1+\sqrt{5}\delta_{k+s(k)}},$
(3.8)

where $\gamma^{*}(\approx 0.2274)$ is the unique root of the equation
$5\gamma^{3}+5\gamma^{2}+3\gamma-1=0$ in the interval $(0,1)$, $s(k)$ is given
by (3.3) and
$\eta:=\sqrt{\frac{1+\delta_{k}}{1-2\delta_{k}-\delta_{k+s(k)}}}$. Then the
sequence $\\{x^{p}\\}$ generated by HBOT or HBOTP obeys

${\|x_{S}-x^{p}\|}_{2}\leq C_{1}\theta^{p-1}+C_{2}\|\nu^{\prime}\|_{2},$ (3.9)

where $S:=\mathcal{L}_{k}(x)$, $\nu^{\prime}:=\nu+Ax_{\overline{S}}$, and the
quantities $C_{1},C_{2}$ are defined as

$C_{1}={\|x_{S}-x^{1}\|}_{2}+({\theta}-{b}){\|x_{S}-x^{0}\|}_{2},~{}C_{2}=\frac{2+(1+\delta_{k})\alpha}{(1-\theta)\sqrt{1-2\delta_{k}-\delta_{k+s(k)}}},$
(3.10)

and $\theta:=(b+\sqrt{b^{2}+4\eta\beta})/2<1$ is ensured under the conditions
(3.8) and the constant $b$ is given by

$b:=\eta\left(|1+\beta-\alpha|+\sqrt{5}\alpha\delta_{k+s(k)}\right).$ (3.11)

###### Proof

From (2.2), we have

$u^{p}-x_{S}=(1-\alpha+\beta)(x^{p}-x_{S})+\alpha(I-A^{T}A)(x^{p}-x_{S})-\beta(x^{p-1}-x_{S})+\alpha
A^{T}\nu^{\prime},$ (3.12)

where $S=\mathcal{L}_{k}(x)$ and $\nu^{\prime}=\nu+Ax_{\overline{S}}$. Let
$\widehat{w}\in\mathcal{W}^{k}$ be a $k$-sparse binary vector such that
$\textrm{supp}(x_{S})\subseteq\textrm{supp}(\widehat{w})$. Then
$x_{S}=x_{S}\otimes\widehat{w}$. Since $(x_{S}-u^{p})\otimes\widehat{w}$ is a
$k$-sparse vector and $y=Ax_{S}+\nu^{\prime}$, we have

$\displaystyle{\|y-A(u^{p}\otimes\widehat{w})\|}_{2}=$
$\displaystyle{\|A[x_{S}-u^{p}\otimes\widehat{w}]+\nu^{\prime}\|}_{2}$
$\displaystyle\leq$
$\displaystyle{\|A[(x_{S}-u^{p})\otimes\widehat{w}]\|_{2}+\|\nu^{\prime}\|}_{2}$
$\displaystyle\leq$
$\displaystyle\sqrt{1+\delta_{k}}{\|(x_{S}-u^{p})\otimes\widehat{w}\|_{2}+\|\nu^{\prime}\|}_{2},$
(3.13)

where the last inequality is obtained by using (2.1). From (3.12), one has

$\displaystyle\|(x_{S}-u^{p})\otimes\widehat{w}\|_{2}$
$\displaystyle\leq|1-\alpha+\beta|\cdot\|(x^{p}-x_{S})\otimes\widehat{w}\|_{2}+\alpha\|[(I-A^{T}A)(x^{p}-x_{S})]\otimes\widehat{w}\|_{2}$
$\displaystyle~{}~{}~{}+\beta\|(x^{p-1}-x_{S})\otimes\widehat{w}\|_{2}+\alpha\|(A^{T}\nu^{\prime})\otimes\widehat{w}\|_{2}.$
(3.14)

Since $x_{S},x^{p},\widehat{w}$ are $k$-sparse vectors and
$\textrm{supp}(x_{S})\subseteq\textrm{supp}(\widehat{w})$, by using Lemmas 4
and 1 (ii), we obtain

$\|[(I-A^{T}A)(x^{p}-x_{S})]\otimes\widehat{w}\|_{2}\leq\sqrt{5}\delta_{k+s(k)}\|x^{p}-x_{S}\|_{2}$
(3.15)

and

$\|(A^{T}\nu^{\prime})\otimes\widehat{w}\|_{2}=\|(A^{T}\nu^{\prime})_{\textrm{supp}(\widehat{w})}\|_{2}\leq\sqrt{1+\delta_{k}}\|\nu^{\prime}\|_{2}.$
(3.16)

Substituting (3.15) and (3.16) into (3) yields

$\displaystyle\|(x_{S}-u^{p})\otimes\widehat{w}\|_{2}$
$\displaystyle\leq|1-\alpha+\beta|\cdot\|x^{p}-x_{S}\|_{2}+\alpha\sqrt{5}\delta_{k+s(k)}\|x^{p}-x_{S}\|_{2}$
$\displaystyle~{}~{}+\beta\|x^{p-1}-x_{S}\|_{2}+\alpha\sqrt{1+\delta_{k}}\|\nu^{\prime}\|_{2}$
$\displaystyle=(|1+\beta-\alpha|+\sqrt{5}\alpha\delta_{k+s(k)})\|x^{p}-x_{S}\|_{2}+\beta\|x^{p-1}-x_{S}\|_{2}$
$\displaystyle~{}~{}+\alpha\sqrt{1+\delta_{k}}\|\nu^{\prime}\|_{2}.$

It follows from (3) that

$\displaystyle{\|y-A(u^{p}\otimes\widehat{w})\|}_{2}\leq$
$\displaystyle\sqrt{1+\delta_{k}}\left(|1+\beta-\alpha|+\sqrt{5}\alpha\delta_{k+s(k)}\right)\|x^{p}-x_{S}\|_{2}$
$\displaystyle+\beta\sqrt{1+\delta_{k}}\|x^{p-1}-x_{S}\|_{2}+[1+{(1+\delta_{k})}\alpha]\|\nu^{\prime}\|_{2}.$
(3.17)

Since $x^{p+1}=u^{p}\otimes{w}^{*}$ in HBOT or $x^{p+1}$ is the optimal
solution of (2.4) in HBOTP, the iterates $\\{x^{p}\\}$ generated by HBOT or
HBOTP satisfies

${\|y-Ax^{p+1}\|}_{2}\leq{\|y-A(u^{p}\otimes{w^{*}})\|}_{2}\leq{\|y-A(u^{p}\otimes{w})\|}_{2}$
(3.18)

for all $w\in\mathcal{W}^{k},$ where the second inequality follows from (2.3).
For $\widehat{w}\in\mathcal{W}^{k}$, it follows from (3.18) that

${\|y-Ax^{p+1}\|}_{2}\leq{\|y-A(u^{p}\otimes\widehat{w})\|}_{2}.$ (3.19)

As $x_{S}-x^{p+1}$ is a $(2k)$-sparse vector, by using Lemma 3, one has

$\displaystyle{\|y-Ax^{p+1}\|}_{2}=$
$\displaystyle{\|A(x_{S}-x^{p+1})+\nu^{\prime}\|}_{2}$ $\displaystyle\geq$
$\displaystyle{\|A(x_{S}-x^{p+1})\|_{2}-\|\nu^{\prime}\|}_{2}$
$\displaystyle\geq$
$\displaystyle\sqrt{1-2\delta_{k}-\delta_{k+s(k)}}{\|x_{S}-x^{p+1}\|_{2}-\|\nu^{\prime}\|}_{2}.$
(3.20)

Combining (3), (3.19) and (3) yields

$\displaystyle\|x^{p+1}-x_{S}\|_{2}$
$\displaystyle\leq\eta(|1+\beta-\alpha|+\sqrt{5}\alpha\delta_{k+s(k)})\|x^{p}-x_{S}\|_{2}+\eta\beta\|x^{p-1}-x_{S}\|_{2}$
$\displaystyle~{}~{}~{}+\frac{2+{(1+\delta_{k})}\alpha}{\sqrt{1-2\delta_{k}-\delta_{k+s(k)}}}\|\nu^{\prime}\|_{2}$
$\displaystyle=b\|x^{p}-x_{S}\|_{2}+\eta\beta\|x^{p-1}-x_{S}\|_{2}+(1-\theta)C_{2}\|\nu^{\prime}\|_{2},$
(3.21)

where $\eta,b,\theta,C_{2}$ are given exactly as in Theorem 3.1. Since
$\delta_{k}\leq\delta_{k+s(k)}<\gamma^{*}$, we have

$\eta\sqrt{5}\delta_{k+s(k)}=\sqrt{5}\delta_{k+s(k)}\sqrt{\frac{1+\delta_{k}}{1-2\delta_{k}-\delta_{k+s(k)}}}<\sqrt{5}\gamma^{*}\sqrt{\frac{1+\gamma^{*}}{1-3\gamma^{*}}}=1,$

where the last equality follows from the fact that $\gamma^{*}$ is the root of
$5\gamma^{3}+5\gamma^{2}+3\gamma=1$ in $(0,1)$. It implies that
$0<\frac{1+1/\eta}{1+\sqrt{5}\delta_{k+s(k)}}-1$, which shows that the range
of $\beta$ in (3.8) is well defined. Furthermore, the first inequality in
(3.8) implies that

$\frac{1+2\beta-1/\eta}{1-\sqrt{5}\delta_{k+s(k)}}<1+\beta<\frac{1+1/\eta}{1+\sqrt{5}\delta_{k+s(k)}},$

which indicates that the range for $\alpha$ in (3.8) is also well defined.
Combining (3.11) with (3.8), we deduce that

$\displaystyle b=$
$\displaystyle\eta\left(|1+\beta-\alpha|+\sqrt{5}\alpha\delta_{k+s(k)}\right)$
$\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}\eta\left[1+\beta-\alpha(1-\sqrt{5}\delta_{k+s(k)})\right],&\textrm{
if }\frac{1+2\beta-1/\eta}{1-\sqrt{5}\delta_{k+s(k)}}<\alpha\leq 1+\beta,\\\
\eta\left[-1-\beta+\alpha(1+\sqrt{5}\delta_{k+s(k)})\right],&\textrm{ if
}1+\beta<\alpha<\frac{1+1/\eta}{1+\sqrt{5}\delta_{k+s(k)}},\end{array}\right.$
$\displaystyle<$ $\displaystyle 1-\eta\beta,$

which means that the relation (3) obeys the conditions of Lemma 2. It follows
from Lemma 2 that (3.9) holds with
$\theta=\frac{b+\sqrt{b^{2}+4\eta\beta}}{2}<1$ and $C_{1},C_{2}$ are given by
(3.10). $\Box$

The error bound (3.9) indicates that the iterates $x^{p}$ generated by the
algorithms can approximate the significant components $x_{S}$ of the solution
of the linear inverse problem. In particular, we immediately obtain the
following convergence result for the algorithms.

###### Corollary 1

Let $x\in\mathbb{R}^{n}$ be a $k$-sparse solution to the system $y=Ax.$ Assume
that the RIC, $\delta_{k+s(k)},$ of $A$ and the parameters $\alpha,\beta$ in
HBOT and HBOTP satisfy $\delta_{k+s(k)}<\gamma^{*},$ where $\gamma^{*}$ is the
same constant given in Theorem 3.1. Then the sequence $\\{x^{p}\\}$ generated
by HBOT or HBOTP obeys that ${\|x-x^{p}\|}_{2}\leq C_{1}\theta^{p-1},$ where
the constant $C_{1}$ is defined as in Theorem 3.1. Thus the iterates
$\\{x^{p}\\}$ generated by HBOT or HBOTP converges to $x.$

##  4 Analysis of HBROT$\omega$ and HBROTP$\omega$

In this section, we establish the error bounds for the two algorithms
HBROT$\omega$ and HBROTP$\omega.$ The analysis is far from being trivial. We
need a few useful technical results before we actually establish the error
bounds. We first recall a helpful lemma concerning the polytope
$\mathcal{P}^{k}$, which is the special case of Lemma 4.2 with $\tau=k$ in
ZL21 .

###### Lemma 5

_ZL21_ Given an index set $\Lambda\subseteq N$ and a vector
$w\in\mathcal{P}^{k}$. Decompose $w_{\Lambda}$ as the sum of $k$-sparse
vectors: $w_{\Lambda}=\sum_{j=1}^{q}w_{\Lambda_{j}}$, where
$q:=\lceil\frac{|\Lambda|}{k}\rceil$, $\Lambda=\bigcup_{j=1}^{q}{\Lambda_{j}}$
and $\Lambda_{1}:=\mathcal{L}_{k}(\omega_{\Lambda})$,
$\Lambda_{2}:=\mathcal{L}_{k}(\omega_{\Lambda\setminus\Lambda_{1}})$, and so
on. Then

$\sum_{j=1}^{q}\|w_{\Lambda_{j}}\|_{\infty}<2.$

We now give an inequality concerning the norms
$\|\cdot\|_{2},~{}\|\cdot\|_{1}$ and $\|\cdot\|_{\infty}.$ This inequality is
an modification of Lemma 6.14 in FR13 , but tailored to the need of the later
analysis in this paper.

###### Lemma 6

Let $h\in\mathbb{R}^{r}\setminus\\{0\\}$ be a vector with positive integer
$r\geq 2$, and $\zeta_{1}>\zeta_{2}$ be two positive numbers such that
$\|h\|_{1}\leq\zeta_{1}$ and $\|h\|_{\infty}\leq\zeta_{2}$. Then

$\|h\|_{2}\leq\left\\{\begin{array}[]{cl}g(r),&~{}\textrm{ if }~{}r\leq
t_{0},\\\ \min\\{g(t_{0}),g(t_{0}+1)\\},&~{}\textrm{ if }~{}r\geq t_{0}+1,\\\
\end{array}\right.$ (4.1)

where $t_{0}:=\lfloor\frac{4\zeta_{1}}{\zeta_{2}}\rfloor$ and

$g(j):=\frac{1}{\sqrt{j}}\zeta_{1}+\frac{\sqrt{j}}{4}\zeta_{2},~{}~{}j\in(0,+\infty),$
(4.2)

is strictly decreasing in the interval $(0,\frac{4\zeta_{1}}{\zeta_{2}}]$ and
strictly increasing in the interval $[\frac{4\zeta_{1}}{\zeta_{2}},+\infty)$.

###### Proof

Without loss of generality, we assume that $h$ is a nonnegative vector since
the norms of $h$ are independent of the sign of every elements of it. Sort the
components of $h$ into descending order, and denote such ordered components by
$z_{1}\geq z_{2}\geq\cdots\geq z_{r}\geq 0$ and $z=(z_{1},\ldots,z_{r})^{T}$.
Thus, $\|z\|_{q}=\|h\|_{q}$ for $q\geq 1.$ For a given positive integer $s$
and $a_{1}\geq a_{2}\geq\cdots\geq a_{s}\geq 0$, from (FR13, , Lemma 6.14),
one has

$\sqrt{a_{1}^{2}+\cdots+a_{s}^{2}}\leq\frac{a_{1}+\cdots+a_{s}}{\sqrt{s}}+\frac{\sqrt{s}}{4}(a_{1}-a_{s}).$
(4.3)

The further discussion is divided into two cases according to the relation
between $r$ and $t_{0}$.

Case 1. For $r\leq t_{0}$, by using (4.2) and (4.3), we have

$\|z\|_{2}\leq\frac{\|z\|_{1}}{\sqrt{r}}+\frac{\sqrt{r}}{4}(z_{1}-z_{r})\leq\frac{\|z\|_{1}}{\sqrt{r}}+\frac{\sqrt{r}}{4}\|z\|_{\infty}\leq\frac{1}{\sqrt{r}}\zeta_{1}+\frac{\sqrt{r}}{4}\zeta_{2}=g(r).$
(4.4)

Case 2. For $r\geq t_{0}+1$, denote by
$t:=\arg\min_{j}\\{g(j):j=t_{0},t_{0}+1\\}$ and let $r_{1},r_{2}$ be
nonnegative integers such that $r=r_{1}t+r_{2}(0\leq r_{2}<t)$. Decompose $z$
as the sum of $t$-sparse vectors: $z=\sum_{j=1}^{r_{1}+1}z_{Q_{j}}$, where
$Q_{j}:=\\{(j-1)t+1,\ldots,jt\\}(j=1,\ldots,r_{1})$ and
$Q_{r_{1}+1}:=\\{r_{1}t+1,\ldots,r_{1}t+r_{2}\\}$.

Firstly, we consider the case when $r_{2}>0$. With the aid of (4.3), we see
that

$\|z_{Q_{j}}\|_{2}\leq\frac{\|z_{Q_{j}}\|_{1}}{\sqrt{t}}+\frac{\sqrt{t}}{4}\left(z_{(j-1)t+1}-z_{jt}\right),\
\ j=1,\ldots,r_{1},$ (4.5)

and

$\|z_{Q_{r_{1}+1}}\|_{2}\leq\frac{\|z_{Q_{r_{1}+1}}\|_{1}}{\sqrt{t}}+\frac{\sqrt{t}}{4}z_{r_{1}t+1},$
(4.6)

which is ensured under the condition
$a_{1}=z_{r_{1}t+1},\ldots,a_{r_{2}}=z_{r_{1}t+r_{2}},a_{r_{2}+1}=\ldots=a_{t}=0$.
Merging (4.5) with (4.6), one has

$\|z\|_{2}=\|\sum_{j=1}^{r_{1}+1}z_{Q_{j}}\|_{2}\leq\sum_{j=1}^{r_{1}+1}\|z_{Q_{j}}\|_{2}\leq\frac{1}{\sqrt{t}}\sum_{j=1}^{r_{1}+1}\|z_{Q_{j}}\|_{1}+\frac{\sqrt{t}}{4}\mu$

with

$\mu:=\sum_{j=1}^{r_{1}}\left(z_{(j-1)t+1}-z_{jt}\right)+z_{r_{1}t+1}=z_{1}-\sum_{j=1}^{r_{1}}\left(z_{jt}-z_{jt+1}\right)\leq
z_{1},$

where the inequality is resulted from $z_{1}\geq z_{2}\geq\cdots\geq z_{r}$.
It follows that

$\|z\|_{2}\leq\frac{1}{\sqrt{t}}\|z\|_{1}+\frac{\sqrt{t}}{4}z_{1}\leq\frac{1}{\sqrt{t}}\zeta_{1}+\frac{\sqrt{t}}{4}\zeta_{2}=g(t),$
(4.7)

where the second inequality is ensured by $\|z\|_{1}=\|h\|_{1}\leq\zeta_{1}$
and $z_{1}=\|h\|_{\infty}\leq\zeta_{2}$.

Secondly, we now consider the case $r_{2}=0$, which means that
$Q_{r_{1}+1}=\emptyset$ and $z=\sum_{j=1}^{r_{1}}z_{Q_{j}}$. Hence, by using
(4.3), we obtain

$\|z\|_{2}\leq\sum_{j=1}^{r_{1}}\|z_{Q_{j}}\|_{2}\leq\frac{1}{\sqrt{t}}\sum_{j=1}^{r_{1}}\|z_{Q_{j}}\|_{1}+\frac{\sqrt{t}}{4}\sum_{j=1}^{r_{1}}\left(z_{(j-1)t+1}-z_{jt}\right).$
(4.8)

For $z_{1}\geq z_{2}\geq\cdots\geq z_{r}\geq 0$, we have

$\sum_{j=1}^{r_{1}}\left(z_{(j-1)t+1}-z_{jt}\right)=z_{1}-\sum_{j=1}^{r_{1}-1}\left(z_{jt}-z_{jt+1}\right)-z_{r_{1}}\leq
z_{1}=\|z\|_{\infty}.$ (4.9)

Merging (4.8) with (4.9) leads to

$\|z\|_{2}\leq\frac{1}{\sqrt{t}}\|z\|_{1}+\frac{\sqrt{t}}{4}\|z\|_{\infty}\leq\frac{1}{\sqrt{t}}\zeta_{1}+\frac{\sqrt{t}}{4}\zeta_{2}=g(t).$
(4.10)

Combining (4.4), (4.7), (4.10) with $\|z\|_{q}=\|h\|_{q}(q\geq 1)$, we obtain
the relation (4.1) directly. $\Box$

Now, we use an example to show that the upper bound of $\|h\|_{2}$ in Lemma 6
is tighter than that of Lemma 6.14 in FR13  in some situations.

###### Example 1

Let
$h=(1,\epsilon_{1},\ldots,\epsilon_{14},\epsilon_{0})^{T}\in\mathbb{R}^{16}$,
where $\epsilon_{j}\geq\epsilon_{0}~{}(j=1,\ldots,14),$
$\sum_{j=1}^{14}\epsilon_{j}=1-\epsilon_{0}$ and $\epsilon_{0}\in(0,1/15]$.
Hence, $\|h\|_{1}=2$ and $\|h\|_{\infty}=1$. Set $\zeta_{1}=2$ and
$\zeta_{2}=1$, then $t_{0}=\frac{4\zeta_{1}}{\zeta_{2}}=8$. The upper bound of
$\|h\|_{2}$ can be given by $\|h\|_{2}\leq 1.5-\epsilon_{0}$ in (4.3) and
$\|h\|_{2}\leq g(8)=\sqrt{2}$ in (4.1), respectively. Since
$1.5-\epsilon_{0}>\sqrt{2}$ for $\epsilon_{0}\in(0,1/15]$, we conclude that
the upper bound of $\|h\|_{2}$ given by (4.1) is tighter than that of (4.3) if
$r>t_{0}+1$ and $\epsilon_{0}=\min_{1\leq i\leq r}|h_{i}|$ is small enough.

Taking
$h=(\|w_{\Lambda_{1}}\|_{\infty},\ldots,\|w_{\Lambda_{q}}\|_{\infty})^{T}$
with $q=\lceil\frac{|\Lambda|}{k}\rceil$, by using Lemma 5, we have
$\|h\|_{1}<2$ and $\|h\|_{\infty}=\|w_{\Lambda_{1}}\|_{\infty}\leq 1$. Hence,
by setting $\zeta_{1}=2$ and $\zeta_{2}=1$ , we get $t_{0}=8$ in Lemma 6,
which results in the following corollary.

###### Corollary 2

Under the conditions of Lemma 5, one has
$\left(\sum_{j=1}^{q}\|w_{\Lambda_{j}}\|_{\infty}^{2}\right)^{1/2}\leq\xi_{q},$
where

$\displaystyle\xi_{q}=\left\\{\begin{array}[]{cl}1,&~{}\textrm{ if }~{}q=1,\\\
\frac{2}{\sqrt{q}}+\frac{\sqrt{q}}{4},&~{}\textrm{ if }~{}2\leq q<8,\\\
\sqrt{2},&~{}\textrm{ if }~{}q\geq 8,\end{array}\right.$ (4.14)

which is strictly decreasing in the interval $[2,8]$ and $\underset{q\geq
1}{max}~{}\xi_{q}=\xi_{2}=\frac{5}{4}\sqrt{2}$.

Using Lemmas 5 and 6, we now establish a helpful inequality in the following
lemma.

###### Lemma 7

Let $x\in\mathbb{R}^{n}$ be a vector satisfying $y=Ax+\nu$ where $\nu$ is a
noise vector. Let $S=\mathcal{L}_{k}(x)$ and $V\subseteq N$ be any given index
set obeys $S\subseteq V$. At the iterate $x^{p}$, the vectors
$w^{(j)}(j=1,\cdots,\omega)$ are generated by HBROT$\omega$ or HBROTP$\omega$.
For every $i\in\\{1,\ldots,\omega\\}$, we have

$\displaystyle\Theta^{i}:=$
$\displaystyle\|A[(u^{p}-x_{S})\otimes(\bigotimes_{j=1}^{i}w^{(j)})]_{\overline{V}}\|_{2}$
$\displaystyle\leq$
$\displaystyle\sqrt{1+\delta_{k}}\Big{[}\big{(}\xi_{q}|1-\alpha+\beta|+2\alpha\delta_{3k}\big{)}\|x^{p}-x_{S}\|_{2}+\beta\xi_{q}\|x^{p-1}-x_{S}\|_{2}$
$\displaystyle~{}~{}~{}+2\alpha\sqrt{1+\delta_{k}}\|\nu^{\prime}\|_{2}\Big{]},$
(4.15)

where $q=\lceil\frac{n-|V|}{k}\rceil$ and $\xi_{q}$ is given by (4.14).

###### Proof

Taking $w=w^{(1)}$ and $\Lambda=\overline{V}$ in Lemma 5 and Corollary 2, one
has

$\sum_{j=1}^{q}\|(w^{(1)})_{\Lambda_{j}}\|_{\infty}<2,~{}\sqrt{\sum_{j=1}^{q}\|(w^{(1)})_{\Lambda_{j}}\|_{\infty}^{2}}\leq\xi_{q},$
(4.16)

where $q=\lceil\frac{n-|V|}{k}\rceil$ and the definition of
${\Lambda_{j}}(j=1,\ldots,q)$ can be found in Lemma 5. Define the $k$-sparse
vectors
$z^{(l)}:=[(u^{p}-x_{S})\otimes(\bigotimes_{j=1}^{i}w^{(j)})]_{\Lambda_{l}}(l=1,\ldots,q)$,
where $u^{p}$ is given by (2.2). Since
$w^{(j)}\in\mathcal{P}^{k}(j=1,\ldots,\omega)$, we have

$\|z^{(l)}\|_{2}\leq\|(\bigotimes_{j=1}^{i}w^{(j)})_{\Lambda_{l}}\|_{\infty}\cdot\|(u^{p}-x_{S})_{\Lambda_{l}}\|_{2}\leq\|(w^{(1)})_{\Lambda_{l}}\|_{\infty}\cdot\|(u^{p}-x_{S})_{\Lambda_{l}}\|_{2},$
(4.17)

where the second inequality follows from $0\leq w^{(j)}\leq\bf e$ for
$j=1,\ldots,\omega$. Since $z^{(l)}(l=1,\ldots,q)$ are $k$-sparse vectors,
from the definition of $\Theta^{i}$ in (7) and (4.17), we obtain

$\displaystyle\Theta^{i}$
$\displaystyle=\|A\sum_{l=1}^{q}z^{(l)}\|_{2}\leq\sum_{l=1}^{q}\|Az^{(l)}\|_{2}\leq\sqrt{1+\delta_{k}}\sum_{l=1}^{q}\|z^{(l)}\|_{2}$
$\displaystyle\leq\sqrt{1+\delta_{k}}\sum_{l=1}^{q}\|(w^{(1)})_{\Lambda_{l}}\|_{\infty}\cdot\|(u^{p}-x_{S})_{\Lambda_{l}}\|_{2},$
(4.18)

where the second inequality is given by (2.1). Since $|\Lambda_{l}|\leq k$ and
$|\textrm{supp}(x^{p}-x_{S})\cup\Lambda_{l}|\leq 3k(l=1,\dots,q)$, by using
(3.12) and Lemma 1, we have

$\displaystyle\|(u^{p}-x_{S})_{\Lambda_{l}}\|_{2}\leq$
$\displaystyle|1-\alpha+\beta|\cdot\|(x^{p}-x_{S})_{\Lambda_{l}}\|_{2}+\alpha\|[(I-A^{T}A)(x^{p}-x_{S})]_{\Lambda_{l}}\|_{2}$
$\displaystyle+\beta\|(x^{p-1}-x_{S})_{\Lambda_{l}}\|_{2}+\alpha\|(A^{T}\nu^{\prime})_{\Lambda_{l}}\|_{2}$
$\displaystyle\leq$
$\displaystyle|1-\alpha+\beta|\cdot\|(x^{p}-x_{S})_{\Lambda_{l}}\|_{2}+\alpha\delta_{3k}\|x^{p}-x_{S}\|_{2}$
$\displaystyle+\beta\|(x^{p-1}-x_{S})_{\Lambda_{l}}\|_{2}+\alpha\sqrt{1+\delta_{k}}\|\nu^{\prime}\|_{2}.$
(4.19)

Substituting (4) into (4) yields

$\displaystyle\frac{\Theta^{i}}{\sqrt{1+\delta_{k}}}\leq$
$\displaystyle|1-\alpha+\beta|\cdot\sum_{l=1}^{q}\|(w^{(1)})_{\Lambda_{l}}\|_{\infty}\cdot\|(x^{p}-x_{S})_{\Lambda_{l}}\|_{2}$
$\displaystyle+\alpha\delta_{3k}\sum_{l=1}^{q}\|(w^{(1)})_{\Lambda_{l}}\|_{\infty}\cdot\|x^{p}-x_{S}\|_{2}$
$\displaystyle+\beta\sum_{l=1}^{q}\|(w^{(1)})_{\Lambda_{l}}\|_{\infty}\cdot\|(x^{p-1}-x_{S})_{\Lambda_{l}}\|_{2}$
$\displaystyle+\alpha\sqrt{1+\delta_{k}}\sum_{l=1}^{q}\|(w^{(1)})_{\Lambda_{l}}\|_{\infty}\cdot\|\nu^{\prime}\|_{2}.$

It follows from Cauchy-Schwarz inequality and (4.16) that

$\displaystyle\frac{\Theta^{i}}{\sqrt{1+\delta_{k}}}$
$\displaystyle\leq|1-\alpha+\beta|\sqrt{\sum_{l=1}^{q}\|(w^{(1)})_{\Lambda_{l}}\|_{\infty}^{2}}\sqrt{\sum_{l=1}^{q}\|(x^{p}-x_{S})_{\Lambda_{l}}\|_{2}^{2}}+2\alpha\delta_{3k}\|x^{p}-x_{S}\|_{2}$
$\displaystyle~{}~{}~{}+\beta\sqrt{\sum_{l=1}^{q}\|(w^{(1)})_{\Lambda_{l}}\|_{\infty}^{2}}\sqrt{\sum_{l=1}^{q}\|(x^{p-1}-x_{S})_{\Lambda_{l}}\|_{2}^{2}}+2\alpha\sqrt{1+\delta_{k}}\|\nu^{\prime}\|_{2}$
$\displaystyle\leq|1-\alpha+\beta|\xi_{q}\|x^{p}-x_{S}\|_{2}+2\alpha\delta_{3k}\|x^{p}-x_{S}\|_{2}+\beta\xi_{q}\|x^{p-1}-x_{S}\|_{2}$
$\displaystyle~{}~{}~{}+2\alpha\sqrt{1+\delta_{k}}\|\nu^{\prime}\|_{2},$

where the last inequality can be obtained from the following relation
$\sum_{j=1}^{q}\|z_{\Lambda_{l}}\|_{2}^{2}=\|z_{\overline{V}}\|_{2}^{2}\leq\|z\|_{2}^{2}$
for any $z\in\mathbb{R}^{n}$ due to
$\overline{V}=\bigcup_{j=1}^{q}{\Lambda_{j}}$ and
$\Lambda_{j}\bigcap\Lambda_{l}=\emptyset(j\neq l)$. Thus (7) holds. $\Box$

We now establish the estimation of
$\|y-A[u^{p}\otimes(\bigotimes_{j=1}^{\omega}w^{(j)})]\|_{2}$ by using Lemma
7.

###### Lemma 8

Let $x\in\mathbb{R}^{n}$ be a vector satisfying $y=Ax+\nu$ where $\nu$ is a
noise vector. At the iterate $x^{p}$, the vectors $u^{p}$ and
$w^{(j)}(j=1,\cdots,\omega)$ are generated by HBROT$\omega$ or HBROTP$\omega$.
Then

$\displaystyle\|y-A[u^{p}\otimes(\bigotimes_{j=1}^{\omega}w^{(j)})]\|_{2}$
$\displaystyle\leq
c_{1,q}\sqrt{1+\delta_{k}}\|x^{p}-x_{S}\|_{2}+\sqrt{1+\delta_{k}}\beta\big{[}\xi_{q}(\omega-1)+1\big{]}\|x^{p-1}-x_{S}\|_{2}$
$\displaystyle~{}~{}~{}+\big{[}\alpha(2\omega-1)(1+\delta_{k})+1\big{]}\|\nu^{\prime}\|_{2},$
(4.20)

where $S:=\mathcal{L}_{k}(x)$, $q=\lceil\frac{n-k}{k}\rceil$, $\xi_{q}$ is
given by (4.14) and $c_{1,q}$ is given as

$\begin{array}[]{rl}c_{1,q}:=\big{(}\xi_{q}(\omega-1)+1\big{)}|1-\alpha+\beta|+\alpha\big{(}2(\omega-1)\delta_{3k}+\delta_{2k}\big{)}.\end{array}$
(4.21)

###### Proof

Let $\widehat{w}\in\mathcal{W}^{k}$ be a binary vector satisfying
$\textrm{supp}(x_{S})\subseteq\textrm{supp}(\widehat{w})$ and
$V=\textrm{supp}(\widehat{w})$. From Lemma 4.3 in ZL21 , we get

$\displaystyle\|y-A[u^{p}\otimes(\bigotimes_{j=1}^{\omega}w^{(j)})]\|_{2}$
$\displaystyle\leq\|y-A(u^{p}\otimes{\widehat{w}})\|_{2}+\sum_{i=1}^{\omega-1}\|A[(u^{p}-x_{S})\otimes(\bigotimes_{j=1}^{i}w^{(j)})\otimes({\bf
e}-\widehat{w})]\|_{2}.$ (4.22)

As $V=\textrm{supp}(\widehat{w})$ and $|V|=k$, it follows from (7) that

$\displaystyle\|A[(u^{p}-x_{S})\otimes(\bigotimes_{j=1}^{i}w^{(j)})\otimes({\bf
e}-\widehat{w})]\|_{2}=\|A[(u^{p}-x_{S})\otimes(\bigotimes_{j=1}^{i}w^{(j)})]_{\overline{\textrm{supp}(\widehat{w})}}\|_{2}$
$\displaystyle\leq\sqrt{1+\delta_{k}}\Big{[}(\xi_{q}|1-\alpha+\beta|+2\alpha\delta_{3k})\|x^{p}-x_{S}\|_{2}+\beta\xi_{q}\|x^{p-1}-x_{S}\|_{2}$
$\displaystyle~{}~{}~{}+2\alpha\sqrt{1+\delta_{k}}\|\nu^{\prime}\|_{2}\Big{]},$
(4.23)

where $q=\lceil\frac{n-k}{k}\rceil$ and $i=1,\ldots,\omega-1$.

We now estimate the term ${\|y-A(u^{p}\otimes\widehat{w})\|}_{2}$ in (4).
Because $|\textrm{supp}(x^{p}-x_{S})\cup\textrm{supp}(\widehat{w})|\leq 2k$,
by using (3) and Lemma 1, we obtain

$\displaystyle\|(x_{S}-u^{p})\otimes\widehat{w}\|_{2}$
$\displaystyle\leq|1-\alpha+\beta|\cdot\|x^{p}-x_{S}\|_{2}+\alpha\|[(I-A^{T}A)(x^{p}-x_{S})]_{\textrm{supp}(\widehat{w})}\|_{2}$
$\displaystyle~{}~{}~{}+\beta\|x^{p-1}-x_{S}\|_{2}+\alpha\|(A^{T}\nu^{\prime})_{\textrm{supp}(\widehat{w})}\|_{2}$
$\displaystyle\leq(|1-\alpha+\beta|+\alpha\delta_{2k})\|x^{p}-x_{S}\|_{2}+\beta\|x^{p-1}-x_{S}\|_{2}+\alpha\sqrt{1+\delta_{k}}\|\nu^{\prime}\|_{2}.$

It follows from (3) that

$\displaystyle{\|y-A(u^{p}\otimes\widehat{w})\|}_{2}$
$\displaystyle\leq\sqrt{1+\delta_{k}}\Big{[}\big{(}|1-\alpha+\beta|+\alpha\delta_{2k}\big{)}\|x^{p}-x_{S}\|_{2}$
$\displaystyle~{}~{}~{}+\beta\|x^{p-1}-x_{S}\|_{2}+\alpha\sqrt{1+\delta_{k}}\|\nu^{\prime}\|_{2}\Big{]}+\|\nu^{\prime}\|_{2}.$
(4.24)

Combining (4), (4) with (4) yields (8). $\Box$

The following property of the hard thresholding operator
$\mathcal{H}_{k}(\cdot)$ is shown in (ZL21, , Lemma 4.1).

###### Lemma 9

_ZL21_ Let $z,h\in\mathbb{R}^{n}$ be two vectors and $\|h\|_{0}\leq k$. Then

$\|h-\mathcal{H}_{k}(z)\|_{2}\leq\|(z-h)_{S\cup
S^{*}}\|_{2}+\|(z-h)_{S^{*}\setminus S}\|_{2},$

where $S:=\textrm{supp}(h)$ and $S^{*}:=\textrm{supp}(\mathcal{H}_{k}(z))$.

Let us state a fundamental property of the orthogonal projection in the
following lemma, which can be found in (F11, , Eq.(3.21)) and (Z20, , p.49),
and was extended to the general case in (ZL20, , Lemma 4.2).

###### Lemma 10

_F11; Z20 ; ZL20 _ Let $x\in\mathbb{R}^{n}$ be a vector satisfying $y=Ax+\nu$
where $\nu$ is a noise vector. Let $S^{*}\subseteq N$ be an index set
satisfying $|S^{*}|\leq k$ and

$z^{*}=\arg\min_{z\in\mathbb{R}^{n}}\\{{\|y-Az\|}_{2}^{2}:\textrm{supp}(z)\subseteq
S^{*}\\}.$

Then

$\|z^{*}-x_{S}\|_{2}\leq\frac{1}{\sqrt{1-(\delta_{2k})^{2}}}\|(z^{*}-x_{S})_{\overline{S^{*}}}\|_{2}+\frac{\sqrt{1+\delta_{k}}}{1-\delta_{2k}}\|\nu^{\prime}\|_{2},$

where $S:=\mathcal{L}_{k}(x)$ and $\nu^{\prime}:=\nu+Ax_{\overline{S}}.$

We now establish the error bounds for HBROT$\omega$ and HBROTP$\omega$.

###### Theorem 4.1

Suppose that $n>3k$ and denote $\sigma:=\lceil\frac{n-2k}{k}\rceil$. Let
$x\in\mathbb{R}^{n}$ be a vector satisfying $y=Ax+\nu$ where $\nu$ is a noise
vector. Denote

$t_{k}:=\frac{\sqrt{1+\delta_{k}}}{\sqrt{1-\delta_{2k}}},\ \
z_{k}:=\sqrt{1-\delta^{2}_{2k}}.$ (4.25)

  * _(i)_

Assume that the $(3k)$-th order RIC, $\delta_{3k}$, of the matrix $A$ and the
pair $(\alpha,\beta)$ satisfy $\delta_{3k}<\gamma^{*}(\omega)$ and

$0\leq\beta<\frac{1-d_{1}}{1+d_{1}+d_{2}},~{}\frac{(d_{0}+d_{2}+2)\beta+d_{0}}{d_{0}-d_{1}+1}<\alpha<\frac{d_{0}+2-(d_{2}-d_{0})\beta}{d_{0}+d_{1}+1},$
(4.26)

where $\gamma^{*}(\omega)$ is the unique root of the equation
$G_{\omega}(\gamma)=1$ in the interval $(0,1)$, where

$G_{\omega}(\gamma):=(2\omega+1)\gamma\sqrt{\frac{1+\gamma}{1-\gamma}}+\gamma,$
(4.27)

and the constants $d_{0},d_{1},d_{2}$ are given as

$\left\\{\begin{array}[]{l}d_{0}:=t_{k}(\omega\xi_{\sigma}+1),\\\
d_{1}:=t_{k}(2\omega\delta_{3k}+\delta_{2k})+\delta_{3k},\\\
d_{2}:=t_{k}[\xi_{\sigma}(\omega-1)+1]\frac{2\omega\delta_{3k}+\delta_{2k}}{2(\omega-1)\delta_{3k}+\delta_{2k}}.\end{array}\right.$
(4.28)

Then, the sequence $\\{x^{p}\\}$ produced by HBROT$\omega$ obeys

$\|x^{p}-x_{S}\|_{2}\leq\theta_{1}^{p-1}\left[\|x^{1}-x_{S}\|_{2}+(\theta_{1}-b_{1})\|x^{0}-x_{S}\|_{2}\right]+\frac{b_{3}}{1-\theta_{1}}\|\nu^{\prime}\|_{2}$
(4.29)

with $\theta_{1}:=\frac{b_{1}+\sqrt{b_{1}^{2}+4b_{2}}}{2}$. The fact
$\theta_{1}<1$ is ensured under (4.26) and $b_{1},b_{2},b_{3}$ are given as

$\displaystyle b_{1}:=$ $\displaystyle
t_{k}c_{\sigma}+\left(|1+\beta-\alpha|+\alpha\delta_{3k}\right),~{}~{}b_{2}:=\beta
t_{k}[\xi_{\sigma}(\omega-1)+1]\frac{c_{\sigma}}{c_{1,\sigma}}+\beta,$
$\displaystyle b_{3}:=$
$\displaystyle\frac{\alpha(2\omega-1)(1+\delta_{k})+2}{\sqrt{1-\delta_{2k}}}\cdot\frac{2\delta_{3k}}{2(\omega-1)\delta_{3k}+\delta_{2k}}\cdot\frac{c_{\sigma}}{c_{\sigma}-c_{1,\sigma}}+\alpha\sqrt{1+\delta_{k}},$
(4.30)

where $c_{1,\sigma}$ and $\xi_{\sigma}$ are given by (4.21) and (4.14),
respectively, and

$c_{\sigma}:=(\omega\xi_{\sigma}+1)|1-\alpha+\beta|+\alpha(2\omega\delta_{3k}+\delta_{2k}).$
(4.31)

  * _(ii)_

Suppose that the $(3k)$-th order RIC, $\delta_{3k}$, of the matrix $A$ and the
pair $(\alpha,\beta)$ satisfy $\delta_{3k}<\gamma^{\sharp}(\omega)$ and

$\small
0\leq\beta<\frac{z_{k}-d_{1}}{1+d_{1}+d_{2}},~{}\frac{(d_{0}+d_{2}+2)\beta+d_{0}+1-z_{k}}{d_{0}-d_{1}+1}<\alpha<\frac{d_{0}+1+z_{k}-(d_{2}-d_{0})\beta}{d_{0}+d_{1}+1},$
(4.32)

where the constants $d_{0},d_{1},d_{2}$ are given by (4.28) and
$\gamma^{\sharp}(\omega)$ is the unique root of the equation
$\frac{1}{\sqrt{1-\gamma^{2}}}G_{\omega}(\gamma)=1$ in the interval $(0,1)$,
where $G_{\omega}(\gamma)$ is given by (4.27). Then, the sequence
$\\{x^{p}\\}$ produced by HBROTP$\omega$ obeys

$\displaystyle\|x^{p}-x_{S}\|_{2}$
$\displaystyle\leq\theta_{2}^{p-1}\left[\|x^{1}-x_{S}\|_{2}+(\theta_{2}-\frac{b_{1}}{z_{k}})\|x^{0}-x_{S}\|_{2}\right]$
$\displaystyle~{}~{}~{}+\frac{1}{1-\theta_{2}}\left(\frac{b_{3}}{z_{k}}+\frac{\sqrt{1+\delta_{k}}}{1-\delta_{2k}}\right)\|\nu^{\prime}\|_{2}$
(4.33)

with $\theta_{2}:=\frac{b_{1}+\sqrt{b_{1}^{2}+4b_{2}z_{k}}}{2z_{k}}$. The fact
$\theta_{2}<1$ is ensured under (4.32) and the constants $b_{i}(i=1,2,3)$ and
$z_{k}$ are given by (_(i)_) and (4.25), respectively.

###### Proof

For
$x^{\sharp}=\mathcal{H}_{k}(u^{p}\otimes(\bigotimes_{j=1}^{\omega}w^{(j)}))$
generated by the Algorithms, by using Lemma 9, we have

$\|x_{S}-x^{\sharp}\|_{2}\leq\|[u^{p}\otimes(\bigotimes_{j=1}^{\omega}w^{(j)})-x_{S}]_{X\cup
S}\|_{2}+\|[u^{p}\otimes(\bigotimes_{j=1}^{\omega}w^{(j)})-x_{S}]_{X\setminus
S}\|_{2},$ (4.34)

where $X=\textrm{supp}(x^{\sharp})$ and $w^{(j)}(j=1,\ldots,\omega)$ are given
by (2.5). Using (3.12) and the triangle inequality, we have that

$\displaystyle\|[u^{p}\otimes(\bigotimes_{j=1}^{\omega}w^{(j)})-x_{S}]_{X\setminus
S}\|_{2}$
$\displaystyle=\|[(u^{p}-x_{S})\otimes(\bigotimes_{j=1}^{\omega}w^{(j)})]_{X\setminus
S}\|_{2}\leq\|(u^{p}-x_{S})_{X\setminus S}\|_{2}$
$\displaystyle\leq|1-\alpha+\beta|\cdot\|(x^{p}-x_{S})_{X\setminus
S}\|_{2}+\alpha\|[(I-A^{T}A)(x^{p}-x_{S})]_{X\setminus S}\|_{2}$
$\displaystyle~{}~{}~{}+\beta\|(x^{p-1}-x_{S})_{X\setminus
S}\|_{2}+\alpha\|(A^{T}\nu^{\prime})_{X\setminus S}\|_{2},$

where the first equality is ensured by $(x_{S})_{X\backslash S}=0$ and the
first inequality is due to $0\leq w^{(j)}\leq{\bf e}(j=1,\ldots,\omega)$.
Since $|X\setminus S|\leq k$ and $|\textrm{supp}(x^{p}-x_{S})\cup(X\setminus
S)|\leq 3k$, by using Lemma 1, we see that

$\displaystyle\|[u^{p}\otimes(\bigotimes_{j=1}^{\omega}w^{(j)})-x_{S}]_{X\setminus
S}\|_{2}$
$\displaystyle\leq\left(|1+\beta-\alpha|+\alpha\delta_{3k}\right)\|x^{p}-x_{S}\|_{2}$
$\displaystyle~{}~{}+\beta\|x^{p-1}-x_{S}\|_{2}+\alpha\sqrt{1+\delta_{k}}\|\nu^{\prime}\|_{2}.$
(4.35)

Denote

$\Theta_{1}:=\|A[u^{p}\otimes(\bigotimes_{j=1}^{\omega}w^{(j)})-x_{S}]_{X\cup
S}\|_{2},\ \
\Theta_{2}:=\|A[u^{p}\otimes(\bigotimes_{j=1}^{\omega}w^{(j)})-x_{S}]_{\overline{X\cup
S}}\|_{2}.\ \ $ (4.36)

As $|X\cup S|\leq 2k$, by using (2.1), we obtain

$\Theta_{1}\geq\sqrt{1-\delta_{2k}}\|(u^{p}\otimes(\bigotimes_{j=1}^{\omega}w^{(j)})-x_{S})_{X\cup
S}\|_{2}.$ (4.37)

For any given $\zeta\in(0,1)$, we consider the following two cases associated
with $\Theta_{1}$ and $\Theta_{2}$.

Case 1.  $\Theta_{2}\leq\zeta\Theta_{1}$. Since $y=Ax_{S}+\nu^{\prime}$, by
the triangle inequality and (4.36), we have

$\displaystyle\|y-A[u^{p}\otimes(\bigotimes_{j=1}^{\omega}w^{(j)})]\|_{2}=\|A[u^{p}\otimes(\bigotimes_{j=1}^{\omega}w^{(j)})-x_{S}]-\nu^{\prime}\|_{2}$
$\displaystyle=\|A[u^{p}\otimes(\bigotimes_{j=1}^{\omega}w^{(j)})-x_{S}]_{X\cup
S}+A[u^{p}\otimes(\bigotimes_{j=1}^{\omega}w^{(j)})-x_{S}]_{\overline{X\cup
S}}-\nu^{\prime}\|_{2}$
$\displaystyle\geq\Theta_{1}-\Theta_{2}-\|\nu^{\prime}\|_{2}$
$\displaystyle\geq(1-\zeta)\Theta_{1}-\|\nu^{\prime}\|_{2}.$ (4.38)

Merging (4.37), (4) with (8) yields

$\displaystyle\|[u^{p}\otimes(\bigotimes_{j=1}^{\omega}w^{(j)})-x_{S}]_{X\cup
S}\|_{2}$
$\displaystyle\leq\frac{1}{(1-\zeta)\sqrt{1-\delta_{2k}}}(\|y-A[u^{p}\otimes\bigotimes_{j=1}^{\omega}w^{(j)}]\|_{2}+\|\nu^{\prime}\|_{2})$
$\displaystyle\leq\frac{t_{k}c_{1,{q_{1}}}}{1-\zeta}\|x^{p}-x_{S}\|_{2}+\frac{\beta
t_{k}}{1-\zeta}\big{[}\xi_{q_{1}}(\omega-1)+1\big{]}\|x^{p-1}-x_{S}\|_{2}$
$\displaystyle~{}~{}+\frac{\alpha(2\omega-1)(1+\delta_{k})+2}{(1-\zeta)\sqrt{1-\delta_{2k}}}\|\nu^{\prime}\|_{2},$
(4.39)

where $q_{1}=\lceil\frac{n-k}{k}\rceil=\sigma+1$ and $t_{k},c_{1,{q_{1}}}$ are
given in (4.25) and (4.21), respectively.

Case 2.  $\Theta_{2}>\zeta\Theta_{1}$. From (4.36) and (4.37), we obtain

$\|[u^{p}\otimes(\bigotimes_{j=1}^{\omega}w^{(j)})-x_{S}]_{X\cup
S}\|_{2}\leq\frac{1}{\zeta\sqrt{1-\delta_{2k}}}\|A[(u^{p}-x_{S})\otimes(\bigotimes_{j=1}^{\omega}w^{(j)})]_{\overline{X\cup
S}}\|_{2}.$ (4.40)

Taking $V=X\cup S$ and $i=\omega$ in (7), one has

$\displaystyle\|A[(u^{p}-x_{S})\otimes(\bigotimes_{j=1}^{\omega}w^{(j)})]_{\overline{X\cup
S}}\|_{2}$
$\displaystyle\leq\sqrt{1+\delta_{k}}\left[c_{2,{q_{2}}}\|x^{p}-x_{S}\|_{2}+\beta\xi_{q_{2}}\|x^{p-1}-x_{S}\|_{2}+2\alpha\sqrt{1+\delta_{k}}\|\nu^{\prime}\|_{2}\right],$
(4.41)

with $q_{2}=\lceil\frac{n-|X\cup S|}{k}\rceil\geq\sigma$ which is due to
$|X\cup S|\leq 2k$, and $c_{2,{q_{2}}}$ is defined as

$c_{2,{q_{2}}}:=\xi_{q_{2}}|1-\alpha+\beta|+2\alpha\delta_{3k}.$ (4.42)

Substituting (4) into (4.40), we get

$\displaystyle\|[u^{p}\otimes(\bigotimes_{j=1}^{\omega}w^{(j)})-x_{S}]_{X\cup
S}\|_{2}$
$\displaystyle\leq\frac{t_{k}}{\zeta}[c_{2,{q_{2}}}\|x^{p}-x_{S}\|_{2}+\beta\xi_{q_{2}}\|x^{p-1}-x_{S}\|_{2}$
$\displaystyle~{}~{}~{}+2\alpha\sqrt{1+\delta_{k}}\|\nu^{\prime}\|_{2}].$
(4.43)

From (4.14), we see that $\xi_{q}$ is decreasing in $[2,n]$. For
$q_{1}=\sigma+1$ and $q_{2}\geq\sigma\geq 2$, we have
$\xi_{q_{1}},\xi_{q_{2}}\leq\xi_{\sigma}$. It follows from (4.21) and (4.42)
that $c_{1,{q_{1}}}\leq c_{1,\sigma}$ and $c_{2,{q_{2}}}\leq c_{2,\sigma}$.
Combining (4) and (4) leads to

$\displaystyle\|(u^{p}\otimes(\bigotimes_{j=1}^{\omega}w^{(j)})-x_{S})_{X\cup
S}\|_{2}$ $\displaystyle\leq
t_{k}\max\left\\{\frac{c_{1,\sigma}}{1-\zeta},\frac{c_{2,\sigma}}{\zeta}\right\\}\|x^{p}-x_{S}\|_{2}$
$\displaystyle~{}~{}~{}+\beta
t_{k}\max\left\\{\frac{\xi_{\sigma}(\omega-1)+1}{1-\zeta},\frac{\xi_{\sigma}}{\zeta}\right\\}\|x^{p-1}-x_{S}\|_{2}$
$\displaystyle~{}~{}~{}+\frac{1}{\sqrt{1-\delta_{2k}}}\max\left\\{\frac{\alpha(2\omega-1)(1+\delta_{k})+2}{1-\zeta},\frac{2\alpha(1+\delta_{k})}{\zeta}\right\\}\|\nu^{\prime}\|_{2}$
(4.44)

for any $\zeta\in(0,1)$.

Next, we select a suitable parameter $\zeta\in(0,1)$ such that the right hand
of (4) is as small as possible. For $\delta_{2k}\leq\delta_{3k}$ and
$\xi_{\sigma}<2$ in (4.14), we have

$\frac{c_{2,\sigma}}{c_{1,\sigma}}=\frac{\xi_{\sigma}|1-\alpha+\beta|+2\alpha\delta_{3k}}{(\xi_{\sigma}(\omega-1)+1)|1-\alpha+\beta|+\alpha[2(\omega-1)\delta_{3k}+\delta_{2k}]}\leq\frac{2\delta_{3k}}{2(\omega-1)\delta_{3k}+\delta_{2k}}.$
(4.45)

It is easy to check that

$\min_{\zeta\in(0,1)}\max\left\\{\frac{c_{1,\sigma}}{1-\zeta},\frac{c_{2,\sigma}}{\zeta}\right\\}=c_{1,\sigma}+c_{2,\sigma}=c_{\sigma},$
(4.46)

where $c_{\sigma}$ is given by (4.31) and its minimum attains at

$\zeta^{*}=\frac{c_{2,\sigma}}{c_{1,\sigma}+c_{2,\sigma}}=\frac{\xi_{\sigma}|1-\alpha+\beta|+2\alpha\delta_{3k}}{(\omega\xi_{\sigma}+1\big{)}|1-\alpha+\beta|+\alpha(2\omega\delta_{3k}+\delta_{2k})}.$
(4.47)

That is,

$\max\left\\{\frac{c_{1,\sigma}}{1-\zeta^{*}},\frac{c_{2,\sigma}}{\zeta^{*}}\right\\}=c_{\sigma}.$
(4.48)

Moreover, noting that $\xi_{\sigma}<2$ and $\delta_{2k}\leq\delta_{3k}$, we
have $\zeta^{*}\geq\frac{\xi_{\sigma}}{\omega\xi_{\sigma}+1}$. In particular,
by taking $\zeta=\zeta^{*}$ in (4), we deduce that

$\max\left\\{\frac{\xi_{\sigma}(\omega-1)+1}{1-\zeta^{*}},\frac{\xi_{\sigma}}{\zeta^{*}}\right\\}=\frac{\xi_{\sigma}(\omega-1)+1}{1-\zeta^{*}}=\frac{\xi_{\sigma}(\omega-1)+1}{c_{1,\sigma}}c_{\sigma},$

and

$\displaystyle\max\left\\{\frac{\alpha(2\omega-1)(1+\delta_{k})+2}{1-\zeta^{*}},~{}\frac{2\alpha(1+\delta_{k})}{\zeta^{*}}\right\\}$
$\displaystyle=\frac{1}{\zeta^{*}}\max\left\\{[\alpha(2\omega-1)(1+\delta_{k})+2]\frac{\zeta^{*}}{1-\zeta^{*}},~{}2\alpha(1+\delta_{k})\right\\}$
$\displaystyle=\frac{c_{\sigma}}{c_{2,\sigma}}\max\left\\{[\alpha(2\omega-1)(1+\delta_{k})+2]\frac{c_{2,\sigma}}{c_{1,\sigma}},~{}2\alpha(1+\delta_{k})\right\\}$
$\displaystyle\leq\frac{c_{\sigma}}{c_{2,\sigma}}\max\left\\{[\alpha(2\omega-1)(1+\delta_{k})+2]\frac{2\delta_{3k}}{2(\omega-1)\delta_{3k}+\delta_{2k}},~{}2\alpha(1+\delta_{k})\right\\}$
$\displaystyle=\Big{[}\alpha(2\omega-1)(1+\delta_{k})+2\Big{]}\frac{2\delta_{3k}}{2(\omega-1)\delta_{3k}+\delta_{2k}}\cdot\frac{c_{\sigma}}{c_{2,\sigma}},$
(4.49)

where the second equality is given by (4.47), the inequality above follows
from (4.45), and the last equality holds owing to
$\delta_{2k}\leq\delta_{3k}$. Merging (4) with (4.48)-(4), we obtain

$\displaystyle\|[u^{p}\otimes(\bigotimes_{j=1}^{\omega}w^{(j)})-x_{S}]_{X\cup
S}\|_{2}$ $\displaystyle\leq t_{k}c_{\sigma}\|x^{p}-x_{S}\|_{2}+\beta
t_{k}\frac{\xi_{\sigma}(\omega-1)+1}{c_{1,\sigma}}c_{\sigma}\|x^{p-1}-x_{S}\|_{2}$
$\displaystyle~{}~{}~{}+\frac{\alpha(2\omega-1)(1+\delta_{k})+2}{\sqrt{1-\delta_{2k}}}\cdot\frac{2\delta_{3k}}{2(\omega-1)\delta_{3k}+\delta_{2k}}\cdot\frac{c_{\sigma}}{c_{2,\sigma}}\|\nu^{\prime}\|_{2}.$
(4.50)

Combining (4), (4) with (4.34), we have

$\|x_{S}-x^{\sharp}\|_{2}\leq
b_{1}\|x^{p}-x_{S}\|_{2}+b_{2}\|x^{p-1}-x_{S}\|_{2}+b_{3}\|\nu^{\prime}\|_{2},$
(4.51)

where the constants $b_{1},b_{2},b_{3}$ are given by (_(i)_).

Next, we estimate $\|x^{p+1}-x_{S}\|_{2}$ for HBROT$\omega$ and HBROTP$\omega$
based on the relation (4.51).

(i) Since $x^{p+1}=x^{\sharp}$ in HBROT$\omega$, (4.51) becomes

$\|x^{p+1}-x_{S}\|_{2}\leq
b_{1}\|x^{p}-x_{S}\|_{2}+b_{2}\|x^{p-1}-x_{S}\|_{2}+b_{3}\|\nu^{\prime}\|_{2}.$
(4.52)

Now, we consider the conditions of Lemma 2. Merging (4.45) with (4.46)
produces

$\frac{c_{\sigma}}{c_{1,\sigma}}\leq\frac{2\omega\delta_{3k}+\delta_{2k}}{2(\omega-1)\delta_{3k}+\delta_{2k}}.$

It follows from (_(i)_) and (4.31) that

$\displaystyle b_{1}+b_{2}$ $\displaystyle\leq
t_{k}c_{\sigma}+\left(|1+\beta-\alpha|+\alpha\delta_{3k}\right)$
$\displaystyle~{}~{}~{}+\left\\{t_{k}\frac{2\omega\delta_{3k}+\delta_{2k}}{2(\omega-1)\delta_{3k}+\delta_{2k}}[\xi_{\sigma}(\omega-1)+1]+1\right\\}\beta=F(\alpha,\beta),$
(4.53)

where

$\displaystyle F(\alpha,\beta):=$
$\displaystyle(d_{0}+1)|1-\alpha+\beta|+d_{1}\alpha+(d_{2}+1)\beta,$
$\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}-(d_{0}-d_{1}+1)\alpha+(d_{0}+d_{2}+2)\beta+d_{0}+1,&\textrm{
if }\ \alpha\leq 1+\beta,\\\
(d_{0}+d_{1}+1)\alpha+(d_{2}-d_{0})\beta-(d_{0}+1),&\textrm{ if }\
\alpha>1+\beta,\end{array}\right.$ (4.56)

with the constants $d_{0},d_{1},d_{2}$ are given by (4.28).

Based on the fact
$\delta_{k}\leq\delta_{2k}\leq\delta_{3k}<\gamma^{*}(\omega)$ and the function
$G_{\omega}(\gamma)$ in (4.27) is strictly increasing in the interval $(0,1)$,
by using (4.28), we have

$d_{1}\leq[(2\omega+1)t_{k}+1]\delta_{3k}\leq
G_{\omega}(\delta_{3k})<G_{\omega}(\gamma^{*}(\omega))=1,$ (4.57)

which shows that the range of $\beta$ in (4.26) is well defined. From the
first inequality in (4.26), we see that

$\frac{(d_{0}+d_{2}+2)\beta+d_{0}}{d_{0}-d_{1}+1}<1+\beta<\frac{d_{0}+2-(d_{2}-d_{0})\beta}{d_{0}+d_{1}+1},$
(4.58)

which implies that the range of $\alpha$ in (4.26) is also well defined.
Merging (4)-(4.58) with the second inequality in (4.26), we see that if
$\frac{(d_{0}+d_{2}+2)\beta+d_{0}}{d_{0}-d_{1}+1}<\alpha\leq 1+\beta,$ then

$F(\alpha,\beta)<-(d_{0}-d_{1}+1)\frac{(d_{0}+d_{2}+2)\beta+d_{0}}{d_{0}-d_{1}+1}+(d_{0}+d_{2}+2)\beta+d_{0}+1=1,$

and if $1+\beta<\alpha<\frac{d_{0}+2-(d_{2}-d_{0})\beta}{d_{0}+d_{1}+1},$ then

$F(\alpha,\beta)<(d_{0}+d_{1}+1)\frac{d_{0}+2-(d_{2}-d_{0})\beta}{d_{0}+d_{1}+1}+(d_{2}-d_{0})\beta-(d_{0}+1)=1.$

It follows from (4) that $b_{1}+b_{2}<1$. Hence, applying Lemma 2 to the
relation (4.52), we conclude that (4.29) holds with
$\theta_{1}=\frac{b_{1}+\sqrt{b_{1}^{2}+4b_{2}}}{2}<1$.

(ii) Since $x^{p+1}$ is given by (2.6) and $S^{p+1}=\textrm{supp}(x^{\sharp})$
in HBROTP$\omega$, by setting $S^{*}=S^{p+1}$ and $z^{*}=x^{p+1}$ in Lemma 10,
we have

$\displaystyle\|x^{p+1}-x_{S}\|_{2}$
$\displaystyle\leq\frac{1}{z_{k}}\|(x^{\sharp}-x_{S})_{\overline{S^{p+1}}}\|_{2}+\frac{\sqrt{1+\delta_{k}}}{1-\delta_{2k}}\|\nu^{\prime}\|_{2}$
$\displaystyle~{}~{}~{}\leq\frac{1}{z_{k}}\|x^{\sharp}-x_{S}\|_{2}+\frac{\sqrt{1+\delta_{k}}}{1-\delta_{2k}}\|\nu^{\prime}\|_{2},$
(4.59)

where $z_{k}$ is given in (4.25) and the first inequality follows from the
fact $(x^{p+1})_{\overline{S^{p+1}}}=(x^{\sharp})_{\overline{S^{p+1}}}=0$.
Combining (4) with (4.51), we have

$\displaystyle\|x^{p+1}-x_{S}\|_{2}$
$\displaystyle\leq\frac{b_{1}}{z_{k}}\|x^{p}-x_{S}\|_{2}+\frac{b_{2}}{z_{k}}\|x^{p-1}-x_{S}\|_{2}+\left(\frac{b_{3}}{z_{k}}+\frac{\sqrt{1+\delta_{k}}}{1-\delta_{2k}}\right)\|\nu^{\prime}\|_{2}.$
(4.60)

Similar to the analysis in Part (i), we need to show
$\frac{b_{1}}{z_{k}}+\frac{b_{2}}{z_{k}}<1$.

From (4.32), we have $\delta_{2k}\leq\delta_{3k}<\gamma^{\sharp}(\omega)$.
Since the function $G_{\omega}(\gamma)$ in (4.27) is strictly increasing in
$(0,1)$, one has

$d_{1}\leq
G_{\omega}(\delta_{3k})<G_{\omega}(\gamma^{\sharp}(\omega))=\sqrt{1-(\gamma^{\sharp}(\omega))^{2}}<\sqrt{1-(\delta_{2k})^{2}}=z_{k},$

where the first inequality is given by (4.57), the first equality follows from
the fact that $\gamma^{\sharp}(\omega)$ is the root of
$\frac{1}{\sqrt{1-\gamma^{2}}}G_{\omega}(\gamma)=1$ in $(0,1)$ and the last
equality is given by (4.25). It follows that the range of $\beta$ in (4.32) is
well defined. From the first inequality in (4.32), we derive

$\frac{(d_{0}+d_{2}+2)\beta+d_{0}+1-z_{k}}{d_{0}-d_{1}+1}<1+\beta<\frac{d_{0}+1+z_{k}-(d_{2}-d_{0})\beta}{d_{0}+d_{1}+1},$
(4.61)

which means that the range of $\alpha$ in (4.32) is well defined. Combining
(4), (4.61) with the second inequality in (4.32) leads to

$\displaystyle F(\alpha,\beta)<$
$\displaystyle\left\\{\begin{array}[]{ll}-(d_{0}-d_{1}+1)\frac{(d_{0}+d_{2}+2)\beta+d_{0}+1-z_{k}}{d_{0}-d_{1}+1}+(d_{0}+d_{2}+2)\beta+d_{0}+1,\\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\textrm{ if
}~{}\frac{(d_{0}+d_{2}+2)\beta+d_{0}+1-z_{k}}{d_{0}-d_{1}+1}<\alpha\leq
1+\beta,\\\ \\\
(d_{0}+d_{1}+1)\frac{d_{0}+1+z_{k}-(d_{2}-d_{0})\beta}{d_{0}+d_{1}+1}+(d_{2}-d_{0})\beta-(d_{0}+1),\\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\textrm{ if
}~{}1+\beta<\alpha<\frac{d_{0}+1+z_{k}-(d_{2}-d_{0})\beta}{d_{0}+d_{1}+1},\end{array}\right.$
$\displaystyle=$ $\displaystyle z_{k}.$

It follows from (4) that $\frac{b_{1}}{z_{k}}+\frac{b_{2}}{z_{k}}<1$.
Therefore, by Lemma 2, it follows from (4.60) that (_(ii)_) holds with
$\theta_{2}=\frac{b_{1}+\sqrt{b_{1}^{2}+4b_{2}z_{k}}}{2z_{k}}<1$. $\Box$

###### Remark 1

  * _(i)_

When $\nu=0$ and $x$ is a $k$-sparse vector, from (4.29) and (_(ii)_), we
observe that the iterates $\\{x^{p}\\}$ generated by HBROT$\omega$ or
HBROTP$\omega$ converges to $x$.

  * _(ii)_

The condition $n>3k$ in Theorem 4.1 can be removed. If so, the constant
$\xi_{\sigma}$ will be replaced by $\underset{q\geq
1}{max}~{}\xi_{q}=\frac{5}{4}\sqrt{2}$ (see Corollary 2). In addition, if
$n>9k$, then $\sigma=\lceil\frac{n-2k}{k}\rceil\geq 8$. In this case, we see
from (4.14) that $\xi_{\sigma}$ in Theorem 4.1 can be replaced by
$\underset{q\geq 2}{min}~{}\xi_{q}=\sqrt{2}$.

  * _(iii)_

In the case $\omega=1$, HBROT$\omega$ and HBROTP$\omega$ reduce to HBROT and
HBROTP, respectively. In this case, the RIP bounds in Theorem 4.1 are reduced
to $\delta_{3k}<\gamma^{*}(1)\approx 0.2118$ for HBROT and
$\delta_{3k}<\gamma^{\sharp}(1)\approx 0.2079$ for HBROTP.

  * _(iv)_

It is not convenient to calculate the RIC of the matrix A and (4.32) is just a
sufficient condition for the theoretical performance of HBROTP$\omega$. In
practical implementation, the parameters $(\alpha,\beta)$ in HBROTP may be set
simply as $0\leq\beta<1/4$ and $\alpha\geq 1+\beta$ to roughly meet the
conditions (4.32).

##  5 Numerical experiments

Sparse signal and image recovery/reconstrcution through measurements
$y=Ax+\nu,$ where $x$ denotes the original signal or image and $\nu$ denotes
the measurement errors, are typical linear inverse problems. In this section,
we provide some experiment results for the HBROTP algorithm proposed in this
paper and compare its performance with several existing methods. The
experiments are performed on a server with the processor Intel(R) Xeon(R) CPU
E5-2680 v3@ 2.50GHz and 256GB memory. All involved convex optimization
problems are solved by CVX GB17  with solver _’Mosek’_ AA20 . The comparison
of five algorithms including HBROTP, ROTP2, PGROTP, $\ell_{1}$-min and OMP are
mainly made via the phase transition with synthetic data and the
reconstruction of a few images are shown in sections 5.1 and 5.2,
respectively.

###  5.1 Phase transition

The first experiment is carried out to compare the performances of the
algorithms through phase transition curves (PTC) BT15 ; BTW15  and average
recovery time. All sparse vectors $x^{*}\in\mathbb{R}^{n}$ and matrices
$A\in\mathbb{R}^{m\times n}$ are randomly generated, and the position of
nonzero elements of $x^{*}$ follows the uniform distribution. In addition, all
columns of $A$ are normalized and the entries of $A$ and the nonzeros of
$x^{*}$ are independent and identically distributed random variables following
${\mathcal{N}}(0,1)$. In this experiment, we consider both accurate
measurements $y=Ax^{*}$ and inaccurate measurements $y=Ax^{*}+\epsilon h$ with
fixed $n=1000$, where $\epsilon=5\times 10^{-3}$ is the noise level and
$h\in\mathbb{R}^{m}$ is a normalized standard Gaussian noise. We let HBROTP
start from $x^{1}=x^{0}=0$ with the fixed parameters $\alpha=5$ and
$\beta=0.2$, while other algorithms start from $x^{0}=0.$ The maximum number
of iterations of HBROTP, ROTP2 and PGROTP are set as 50, while OMP is
performed exactly $k$ iterations and $\ell_{1}$-min is performed by the solver
_’Mosek’_ directly. Given the random data $(A,x^{*})$ or $(A,x^{*},h),$ the
recovery is counted as _’success’_ when the criterion

$\|x^{p}-x^{*}\|_{2}/\|x^{*}\|_{2}\leq 10^{-3}$

is satisfied, in which $x^{p}$ is the approximated solution generated by
algorithms.

(a) Accurate measurements

(b) Inaccurate measurements ($\epsilon=5\times 10^{-3}$)

Figure 1: The $50\%$ success rate phase transition curves for algorithms.

Denote by $\kappa=m/n$ and $\rho=k/m$, where $\kappa$ is often called the
sampling rate or the compression ratio. In the $(\kappa,\rho)$-space, the
region below the PTC is called the _’success’_ recovery region, where the
solution of the problem instances can be exactly or approximately recovered,
while the region above the PTC corresponds to the _’failure’_ region. Thus if
the region below the PTC is wider, the performance of an algorithm would be
better. We now briefly describe the mechanism for plotting the PTC which is
taken as the classical 50% logistic regression curves in this paper, and more
detailed information can be found in BT15 ; BTW15 . To generate the PTC, 13
groups of $m=\lceil\kappa\cdot n\rceil$ are considered, where the sampling
rate $\kappa$ is ranged from 0.1 to 0.7 with stepsize 0.05. For any given $m$,
by using the bisection method, the approximated recovery phase transition
region $[k_{\min},k_{\max}]$ is produced for each algorithm, in which the
success rate of recovery is at least 90% as $k<k_{\min}$ and at most 10% as
$k>k_{\max}$. The interval $[k_{\min},k_{\max}]$ will be equally divided into
$\min\\{k_{\max}-k_{\min},50\\}$ parts, and 10 problem instances are tested
for each $k$ to produce the recovery success rate for given algorithm. Thus
PTC can be obtained from the logistic regression model in BT15 ; BTW15
directly.

(a) Time ratio of ROTP2 and HBROTP

(b) Time ratio of PGROTP and HBROTP

(c) Time ratio of $\ell_{1}$-min and HBROTP

(d) Time ratio of OMP and HBROTP

Figure 2: The ratios of average CPU time of the algorithms.

The PTC for the experimented algorithms are shown in Fig. 1(a) and (b), which
correspond to the accurate measurements and inaccurate measurements with the
noise level $\epsilon=5\times 10^{-3}$, respectively. The results indicate
that HBROTP has the highest PTC as $\kappa\leq 0.5$, which shows that the
recovery capability of HBROTP is superior to the mentioned a few existing
algorithms both in noiseless and the above-described noisy environments.
However, the PTC of ROTP2, PGROTP and OMP may approach and might surpass that
of HBROTP with a larger value of $\kappa$. The comparison between Fig. 1(a)
and (b) demonstrates that all algorithms are robust for signal recovery when
the measurements are slightly inaccurate except $\ell_{1}$-min. The comparison
indicates that the overall performance of HBROTP is very comparable to that of
mentioned existing methods.

In the intersection of the recovery regions of multiple algorithms, we compare
the average CPU time for recovery of the algorithms. Specifically, for each
above-mentioned $\kappa$, we test 10 problem instances for each algorithm with
the mesh $(\kappa,\rho)$, wherein $\rho$ is ranged from 0.02 to 1 with
stepsize 0.02 until the success rate of recovery is less than 90%. The ratios
of the average computational time of ROTP2, PGROTP, $\ell_{1}$-min and OMP
against that of HBROTP are displayed in Fig. 2 (a)-(d), respectively. Fig. 2
(a) and (b) show that HBROTP is at least 1.6 times faster than ROTP2 in most
areas and slower than PGROTP except in the region $[0.1,0.2]\times[0.02,0.1]$.
On the other hand, from Fig. 2 (a)-(d), we observe that the ROT-type
algorithms such as HBROTP, ROTP2 and PGROTP are always time-consuming than
$\ell_{1}$-min and OMP, due to solving and resolving quadratic convex
optimization problems.

###  5.2 Image reconstruction

In this section, we compare the performance of algorithms on the
reconstruction of several images (Lena, Peppers and Baboon) of size $512\times
512$. Only accurate measurements are used in the experiment, and the
measurement matrices are $m\times n$ normalized standard Gaussian matrices
with $n=512$ and $m=\lceil\kappa\cdot n\rceil,$ where $\kappa$ is the sampling
rate. The discrete wavelet transform with the _‘sym8’_ wavelet is used to
establish the sparse representation of the images. The input sparsity level is
set as $k=\lceil n/10\rceil$ for HBROTP, ROTP2 and PGROTP. The peak signal-to-
noise ratio (PSNR) is used to compare the reconstruction quality of images,
which is defined by

$PSNR:=10\cdot log_{10}(V^{2}/MSE),$

where $MSE$ denotes the mean-squared error between the reconstructed and
original image, and $V$ represents the maximum fluctuation in the original
image data type ($V=255$ is used in our experiments). Clearly, the larger the
value of PSNR, the higher the reconstruction quality.

Table 1: Comparison of PSNR (dB) for algorithms with different sampling rates | $\kappa$ | HBROTP | ROTP2 | PGROTP |  $\ell_{1}$-min | OMP
---|---|---|---|---|---|---
| 0.3 | 32.60 | 32.34 | 33.12 | 33.63 | 31.37
Lena | 0.4 | 34.37 | 32.49 | 31.75 | 35.10 | 32.95
| 0.5 | 35.63 | 33.11 | 31.93 | 37.04 | 34.34
| 0.3 | 31.31 | 32.33 | 33.27 | 33.03 | 30.17
Peppers | 0.4 | 33.10 | 31.78 | 31.60 | 34.08 | 31.66
| 0.5 | 34.23 | 32.04 | 31.10 | 35.90 | 33.38
| 0.3 | 28.70 | 31.35 | 32.33 | 29.90 | 28.35
Baboon | 0.4 | 29.12 | 30.06 | 30.00 | 30.05 | 28.53
| 0.5 | 29.37 | 30.06 | 30.07 | 30.20 | 28.78

The results in terms of PSNR with the sampling rate $\kappa=0.3,0.4,0.5$ are
summarized in Tab. 1, from which we see that HBROTP is always superior to OMP
and inferior to $\ell_{1}$-min in reconstruction quality. For ROTP-type
algorithms with $\kappa=0.4,0.5$, the PSNR of HBROTP exceed that of ROTP2,
PGROTP at least 1.88 dB for Lena and 1.32 dB for Peppers, respectively. In
other cases, ROTP2 and PGROTP obtained better results than HBROTP in
reconstruction quality. In particular, the performance of ROTP2 and PGROTP are
always equivalent or superior to $\ell_{1}$-min for Baboon. In the meantime,
the performance of HBROTP in the visual qualities with $\kappa=0.3,0.4,0.5$
are displayed in Fig. 3. It can be seen that the reconstruction quality has
been significantly improved for three images as the sampling rate $\kappa$ is
ranged from $0.3$ to $0.5$, and the best visual results have been achieved
around $\kappa=0.5$.

(a) $\begin{array}[]{c}\textrm{Lena}\\\ \textrm{(Original)}\end{array}$

(b) $\kappa=0.3$

(c) $\kappa=0.4$

(d) $\kappa=0.5$

(e) $\begin{array}[]{c}\textrm{Peppers}\\\ \textrm{(Original)}\end{array}$

(f) $\kappa=0.3$

(g) $\kappa=0.4$

(h) $\kappa=0.5$

(i) $\begin{array}[]{c}\textrm{Baboon}\\\ \textrm{(Original)}\end{array}$

(j) $\kappa=0.3$

(k) $\kappa=0.4$

(l) $\kappa=0.5$

Figure 3: Performance of HBROTP for three images with different sampling
rates.

##  6 Conclusions

The new algorithms that combine the optimal $k$-thresholding and heavy-ball
method are proposed in this paper. Such algorithms can be seen as the
acceleration versions of the optimal $k$-thresholding method. The solution
error bounds and convergence of the proposed algorithms have been shown mainly
under the RIP of the matrices. The numerical performance of the proposed
HBROTP algorithm has been evaluated through phase transition, average runtime
and image reconstruction. The experiment results indicate that HBROTP is
robust signal recovery method, especially as the sampling rates $\kappa\leq
0.5.$ The proposed algorithm is generally faster than the standard ROTP method
due to integrating the heavy-ball acceleration technique.

## Declarations and statements

The authors have no relevant financial or non-financial interests to disclose
and have no competing interests to declare that are relevant to the content of
this article. All authors contributed to the conception and analysis in this
work. Specifically, Zhong-Feng Sun contributed in conceptualization,
methodology, formal analysis, the first draft writing, coding and experiment
design; Yun-Bin Zhao contributed in conceptualization, methodology, analysis,
funding acquisition, resources, supervision, experiment design and editing;
Jin-Chuan Zhou contributed to conception, analysis, methodology and reviewing.
The first draft of the manuscript was written by Zhong-Feng Sun supported by
Jin-Chuan Zhou. All authors read, discussed and approved the final manuscript.

## References

  * (1) E. D. Andersen and K. D. Andersen. The MOSEK interior point optimizer for linear programming: an implementation of the homogeneous algorithm. _High performance optimization_. Springer, Boston, MA, 33:197–232, 2000.
  * (2) J. F. Aujol, C. Dossal, and A. Rondepierre. Convergence rates of the heavy-ball method under the Łojasiewicz property. _Math. Program._ , 2022. https://doi.org/10.1007/s10107-022-01770-2
  * (3) A. Beck and M. Teboulle. A fast iterative shrinkage-thresholding algorithm for linear inverse problems. _SIAM J. Imaging Sci._ , 2(1):183–202, 2009.
  * (4) J. D. Blanchard and J. Tanner. Performance comparisons of greedy algorithms in compressed sensing. _Numer. Linear Algebra Appl._ , 22(2):254–282, 2015.
  * (5) J. D. Blanchard, J. Tanner, and K. Wei. CGIHT: conjugate gradient iterative hard thresholding for compressed sensing and matrix completion. _Information and Inference: A Journal of the IMA_ , 4(4):289–327, 2015.
  * (6) T. Blumensath. Accelerated iterative hard thresholding. _Signal Process._ , 92(3):752–756, 2012.
  * (7) T. Blumensath and M. E. Davies. Iterative thresholding for sparse approximations. _J. Fourier Anal. Appl._ , 14(5-6):629–654, 2008.
  * (8) T. Blumensath and M. E. Davies. Normalized iterative hard thresholding: Guaranteed stability and performance. _IEEE J. Sel. Top. Signal Process._ , 4(2):298–309, 2010.
  * (9) M. Borgerding, P. Schniter, and S. Rangan. AMP-inspired deep networks for sparse linear inverse problems. _IEEE Trans. Signal Process._ , 65(16):4293–4308, 2017.
  * (10) C. Buchheim and E. Traversi. Quadratic combinatorial optimization using separable underestimators, _INFORMS J. Comput._ , 30(3):424–437, 2018.
  * (11) E. J. Candès and T. Tao. Decoding by linear programming. _IEEE Trans. Inform. Theory_ , 51(12):4203–4215, 2005.
  * (12) E. J. Candès, M. B. Wakin, and S. P. Boyd. Enhancing sparsity by reweighted $\ell_{1}$-minimization. _J. Fourier Anal. Appl._ , 14:877–905, 2008.
  * (13) W. A. Chaovalitwongse, I. P. Androulakis, and P. M. Pardalos. Quadratic integer programming: Complexity and equivalent forms. In: C. Floudas, P. Pardalos (Eds.) _Encyclopedia of Optimization_ , Springer, Boston, MA, 2008.
  * (14) R. Chartrand. Exact reconstruction of sparse signals via nonconvex minimization. _IEEE Signal Process. Lett._ , 14(10):707–710, 2007.
  * (15) S. S. Chen, D. L. Donoho, and M. A. Saunders. Atomic decomposition by basis pursuit. _SIAM J. Sci. Comput._ , 20(1):33–61, 1998.
  * (16) W. Chen, B. Zhang, S. Jin, B. Ai, and Z. Zhong. Solving sparse linear inverse problems in communication systems: A deep learning approach with adaptive depth. _IEEE J. Sel. Areas Commun._ , 39(1):4–17, 2021.
  * (17) W. Dai and O. Milenkovic. Subspace pursuit for compressive sensing signal reconstruction. _IEEE Trans. Inform. Theory_ , 55(5):2230–2249, 2009.
  * (18) I. Daubechies, M. Defrise, and C. De Mol. An iterative thresholding algorithm for linear inverse problems with a sparsity constraint. _Comm. Pure Appl. Math._ , 57(11):1413–1457, 2004.
  * (19) D. L. Donoho. De-noising by soft-thresholding. _IEEE Trans. Inform. Theory_ , 41(3):613–627, 1995.
  * (20) D. L. Donoho and J. M. Johnstone. Ideal spatial adaptation by wavelet shrinkage. _Biometrika_ , 81(3):425–455, 1994.
  * (21) M. Elad. Why simple shrinkage is still relevant for redundant representations? _IEEE Trans. Inform. Theory_ , 52(12):5559–5569, 2006.
  * (22) M. Elad. _Sparse and redundant representations: From theory to applications in signal and image processing._ Springer, New York, 2010.
  * (23) Y. Eldar and G. Kutyniok. _Compressed Sensing: Theory and Applications_. Cambridge University Press, Cambridge, UK, 2012.
  * (24) S. Foucart. Hard thresholding pursuit: An algorithm for compressive sensing. _SIAM J. Numer. Anal._ , 49(6):2543–2563, 2011.
  * (25) S. Foucart and H. Rauhut. _A mathematical introduction to compressive sensing._ Springer, New York, 2013.
  * (26) M. Grant and S. Boyd. CVX: matlab software for disciplined convex programming. Version 1.21, 2017.
  * (27) M. Gürbüzbalaban, A. Ozdaglar, and P. A. Parrilo. On the convergence rate of incremental aggregated gradient algorithms. _SIAM J. Optim._ , 27(2):1035–1048, 2017.
  * (28) N. Kuru, Ş. İ. Birbil, M. Gürbüzbalaban, and S. Yildirim. Differentially private accelerated optimization algorithms. _SIAM J. Optim._ , 32(2):795–821, 2022.
  * (29) A. Kyrillidis and V. Cevher. Matrix recipes for hard thresholding methods. _J. Math. Imaging Vis._ , 48:235–265, 2014.
  * (30) L. Lessard, B. Recht, and A. K. Packard. Analysis and design of optimization algorithms via integral quadratic constraints. _SIAM J. Optim._ , 26(1):57–95, 2016.
  * (31) S. Li, M. Amin, G. Zhao, and H. Sun. Radar imaging by sparse optimization incorporating MRF clustering prior. _IEEE Geosci. Remote Sens. Lett._ , 17(7):1139–1143, 2020.
  * (32) H. Li, H. Cheng, Z. Wang, and G. C. Wu. Distributed nesterov gradient and heavy-ball double accelerated asynchronous optimization. _IEEE Trans. Neural Netw. Learn. Syst._ , 32(12):5723–5737, 2021.
  * (33) N. Meng and Y. B. Zhao. Newton-step-based hard thresholding algorithms for sparse signal recovery. _IEEE Trans. Signal Process._ , 68:6594–6606, 2020.
  * (34) N. Meng and Y. B. Zhao. Newton-type optimal thresholding algorithms for sparse optimization problems. _J. Oper. Res. Soc. China_ , 10:447–469, 2022.
  * (35) N. Meng, Y. B. Zhao, M. Kočvara, and Z. F. Sun. Partial gradient optimal thresholding algorithms for a class of sparse optimization problems. _J. Global Optim._ , 84:393–413, 2022.
  * (36) H. Mohammadi, M. Razaviyayn, and M. R. Jovanović. Robustness of accelerated first-order algorithms for strongly convex optimization problems. _IEEE Trans. Automat. Contr._ , 66(6):2480–2495, 2021.
  * (37) D. Needell and J. A. Tropp. CoSaMP: Iterative signal recovery from incomplete and inaccurate samples. _Appl. Comput. Harmon. Anal._ , 26(3):301–321, 2009.
  * (38) S. Oymak, B. Recht, and M. Soltanolkotabi. Sharp time-data tradeoffs for linear inverse problems. _IEEE Trans. Inf. Theory_ , 64(6):4129–4158, 2018.
  * (39) B. T. Polyak. Some methods of speeding up the convergence of iteration methods. _USSR Comput. Math. Math. Phys._ , 4(5):1–17, 1964.
  * (40) P. Schniter, L. C. Potter, and J. Ziniel. Fast Bayesian matching pursuit. _2008 Inf. Theory Appl. Workshop_ , 326–333, 2008.
  * (41) T. Tirer and R. Giryes. Back-projection based fidelity term for ill-posed linear inverse problems. _IEEE Trans. Image Process._ , 29:6164–6179, 2020.
  * (42) J. A. Tropp and A. C. Gilbert. Signal recovery from random measurements via orthogonal matching pursuit. _IEEE Trans. Inform. Theory_ , 53(12):4655–4666, 2007.
  * (43) J. A. Tropp and S. J. Wright. Computational methods for sparse solution of linear inverse problems. _Proc. IEEE_ , 98(6):948–958, 2010.
  * (44) V. Ugrinovskii, I. R. Petersen, and I. Shames. Global convergence and asymptotic optimality of the heavy ball method for a class of nonconvex optimization problems. _IEEE Control Syst. Lett._ , 6:2449–2454, 2022.
  * (45) D. P. Wipf and B. D. Rao. Sparse Bayesian learning for basis selection. _IEEE Trans. Signal Process._ , 52(8):2153–2164, 2004.
  * (46) J. Xiang, Y. Dong, and Y. Yang. FISTA-net: Learning a fast iterative shrinkage thresholding network for inverse problems in imaging. _IEEE Trans. Med. Imag._ , 40(5):1329–1339, 2021.
  * (47) R. Xin and U. A. Khan. Distributed heavy-ball: A generalization and acceleration of first-order methods with gradient tracking. _IEEE Trans. Autom. Control_ , 65(6):2627–2633, 2020.
  * (48) Y. B. Zhao. Optimal $k$-thresholding algorithms for sparse optimization problems. _SIAM J. Optim._ , 30(1):31–55, 2020.
  * (49) Y. B. Zhao and Z. Q. Luo. Constructing new reweighted $\ell_{1}$-algorithms for the sparsest points of polyhedral sets. _Math. Oper. Res._ , 42(1):57–76, 2017.
  * (50) Y. B. Zhao and Z. Q. Luo. Improved RIP-based bounds for guaranteed performance of two compressed sensing algorithms. _Sci. China Math.,_ 2022\. https://doi.org/10.1007/s11425-021-1987-2.
  * (51) Y. B. Zhao and Z. Q. Luo. Analysis of optimal thresholding algorithms for compressed sensing. _Signal Process._ , 187:108148, 2021.
  * (52) Y. B. Zhao and Z. Q. Luo. Natural thresholding algorithms for signal recovery with sparsity. _IEEE Open J. Signal Process._ , 3: 417-431, 2022.