remarkRemark hypothesisHypothesis claimClaim Convergence and Semi-convergence
of a Class P-BIMM. Mirzapour, A. Cegielski & T. Elfving

# Convergence and Semi-convergence of a Class of Constrained Block Iterative
Methods

Mahdi Mirzapour Department of Mathematics, Faculty of Sciences, Bu-Ali Sina
University, Hamedan, Iran<EMAIL_ADDRESS>Andrzej Cegielski
Institute of Mathematics, University of Zielona Góra, Zielona Góra, Poland()
<EMAIL_ADDRESS>Tommy Elfving Department of Mathematics,
Linköping University, SE-581 83 Linköping, Sweden<EMAIL_ADDRESS>
###### Abstract

In this paper, we analyze the convergence properties of projected non-
stationary block iterative methods (P-BIM) aiming to find a constrained
solution to large linear, usually both noisy and ill-conditioned, systems of
equations. We split the error of the $k$th iterate into noise error and
iteration error, and consider each error separately. The iteration error is
treated for a more general algorithm, also suited for solving split
feasibility problems in Hilbert space. The results for P-BIM come out as a
special case. The algorithmic step involves projecting onto closed convex
sets. When these sets are polyhedral, and of finite dimension, it is shown
that the algorithm converges linearly. We further derive an upper bound for
the noise error of P-BIM. Based on this bound, we suggest a new strategy for
choosing relaxation parameters, which assist in speeding up the reconstruction
process and improving the quality of obtained images. The relaxation
parameters may depend on the noise. The performance of the suggested strategy
is shown by examples taken from the field of image reconstruction from
projections.

###### keywords:

Landweber type iteration, Block Iterative Method, Split Feasibility Problem,
Relaxation Parameters, Semi-convergence, Constraints, Tomographic Imaging

68Q25, 68R10, 68U05

##  1 Introduction

Among reconstruction methods in computerized tomography (CT) iterative methods
are particularly suited for their ability to incorporate constraints (e.g.
nonnegativity) on the sought solution. This is important in situations with
few and/or noisy data, as in limited data CT when exposition to a low dose of
X-rays is an issue. After discretization of the underlying integral transform
a large, sparse, unstructured and ill-posed (sensitive to noise) linear system
occurs:

(1) $Ax=b,$

where $A\in\mathbb{R}^{m\times n}$ and $b\in\mathbb{R}^{m}$ are known,
$x\in\mathbb{R}^{n}$ is the unknown image to be estimated. We will, in
connection with each algorithm/problem discuss possible consistency
assumptions on the rhs $b$.

###  1.1 Some background

The first reported use of iteration in CT seems to be the paper [29]. There,
in its basic form, the given method turned out to be identical to the Kaczmarz
method [36]. Then for each equation in (1) the new iterate is constructed such
that the previous equation is satisfied. This is done for each of the $m$
equations in turn, and then the process is repeated. For more on the history
see [28]. Somewhat later in [26] a fully simultaneous method, SIRT, was
proposed where the new iterate is constructed by adding to the previous one
the gradient of the least squares functional corresponding to (1). Later
block-iterative versions were developed and analyzed, where, instead of using
a single row, blocks of rows are used in each iteration. A very general form
of block-iteration is analyzed in [1]. For more on the development of block-
iteration see, e.g. the recent survey [20], and further in Subsection 1.4. To
cope with constraints it was, early on, suggested to project the iterate onto
the nonnegative orthant, either after each iterative step or after a whole
cycle, i.e. after completing a full sweep through the matrix $A$. The order in
which the rows/blocks are used is called control, and can have quite an impact
of the behavior of the method [33]. In cyclic control the rows/blocks are
picked up in their original order in the matrix. For some other controls, e.g.
almost cyclic control see [9]. A recent strategy is stochastic control which
has attracted much interest see, e.g. [44, 46, 55]. We will however not
consider this control in our paper.

###  1.2 Simultaneous Algorithms

Let $M$ and $N$ be given symmetric positive definite (SPD) matrices, and
denote by $Q^{1/2}$ the square root of an SPD matrix $Q$. Moreover
$\\{\lambda_{k}\\}_{k=0}^{\infty}$ is a sequence of positive relaxation
parameters, and $\tilde{A}=M^{1/2}AN^{1/2}$. The simultaneous iterative
reconstruction technique (SIRT) is a class of iterative and gradient-based
methods defined as

(2)
$x^{k+1}=x^{k}+\frac{\lambda_{k}}{\|\tilde{A}\|^{2}}NA^{T}M\left(b-Ax^{k}\right)~{}~{}~{}~{}k=0,1,2,\ldots,$

see e.g. [32]. Here $\|\tilde{A}\|=\|\tilde{A}\|_{2}$ (note that
$\|\tilde{A}\|^{2}=\|\tilde{A}^{T}\tilde{A}\|$). By premultiplying (2) by
$N^{-1/2}$ and putting $N^{-1/2}x^{k}=y^{k}$ one gets the equivalent iteration

$y^{k+1}=y^{k}+\frac{\lambda_{k}}{\|\tilde{A}\|^{2}}\tilde{A}^{T}(M^{1/2}b-\tilde{A}y^{k}),$

which as is well known converges for
$\lambda_{k}\in[\varepsilon,2-\varepsilon]$, where $\varepsilon$ is small
positive constant. For the stationary case, $\lambda_{k}=\lambda$ this is also
a necessary condition. It follows that the sequence $\\{x^{k}\\}$ converges
towards an $M-$weighted least squares solution of (1) provided
$\lambda_{k}\in[\varepsilon,2-\varepsilon]$. Hence no consistency assumption
is needed on $b$. For a more general result (for the nonstationary case) see
[52]. Both $M$ and $N$ can have an efficacious impact on the initial speed of
convergence. A computational analysis of SIRT applied to CT-problems appears
in [30].

Let diag$(B)$ be a diagonal matrix with $B_{j}$ being the $j$th diagonal
element, and let Id denote the identity matrix of appropriate dimension (or in
the case of infinite dimensional Hilbert spaces the identity operator).
Further $a^{i}$ is $i$th row of $A$, $||.||$ is the Euclidean 2-norm, and
$||.||_{W}$ is the weighted Euclidean norm. Here $W$ is a given SPD matrix.

Some well-known choices for $M$ and $N$ yielding fully simultaneous iterations
are listed below:

  * •

$N=M=$Id, lead to the Landweber method [5, 25].

  * •

$N=$Id and $M=D_{C}=\frac{1}{m}diag\left(\frac{1}{\|a^{i}\|^{2}}\right)$ lead
to the Cimmino method [19].

  * •

$N=$Id and $M=D_{W}=\frac{1}{m}diag\left(\frac{1}{\|a^{i}\|_{W}^{2}}\right)$
lead to CAV (Component Averaging Method) [16]. Here $W=diag(w_{j})$ where
$w_{j}$ equals the number of nonzero elements in the $j$th column of $A$.

  * •

$N=W^{-1}$ and $M=mD_{C}$ lead to DROP (Diagonally Relaxed Orthogonal
Projection) method [14].

  * •

$N=diag(\mbox{row sums})^{-1}$ and $M=diag(\mbox{column sums})^{-1}$ lead to
SART (Simultaneous Algebraic Reconstruction Technique) [37].

  * •

$N=$Id and $M=(2-\lambda)(D+\lambda L^{T})D(D+\lambda L)^{-1}$ where $L$ is
the left triangular part of $AA^{T}$, $D$ its diagonal and taking
$\lambda_{k}=\lambda$ fixed, lead to the symmetric Kaczmarz’s method [23].

###  1.3 Projected SIRT

Consider the linear system of equations (1) subject to $x\in C$ where $C$
denotes a closed convex set and $P_{C}$ the operator of metric projection onto
the set $C$. Inserting a projection onto a convex set after each iteration can
reduce the error and improve the quality of the reconstructed image, see, e.g.
[49, 51] and Figure 2 below. We will now assume that $N=\textrm{Id}$ (but
remark on the case $N\neq\textrm{Id}$ at the end of Section 2). Put

(3) $\hat{A}=M^{1/2}A,\ \hat{b}=M^{1/2}b.$

Then the P(Projected)-SIRT Algorithm is (note that
$\|\hat{A}\|^{2}=\|A^{T}MA\|$)

Algorithm 1 P-SIRT Algorithm

1:$x^{0}\in\mathbb{R}^{n}$ is an arbitrary starting point.

2:for $k\geq 0$ do
$x^{k+1}=U(x^{k})=P_{C}\left(x^{k}+\frac{\lambda_{k}}{\|\hat{A}\|^{2}}\hat{A}^{T}\left(\hat{b}-\hat{A}x^{k}\right)\right)$

3:end for

Let $\mathop{\mathrm{F}ix}\hat{U}$ denote the set of fixed points of an
operator $\hat{U}\ :\ \mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$. The following
characterization (adapted to our problem formulation) is shown in [6,
Proposition 2.1] (also given in [9, Proposition 4.7.2]),

(4) $\mathop{\mathrm{F}ix}U=\arg\min_{x\in C}\|b-Ax\|^{2}_{M}.$

We have the following convergence result for Algorithm 1.

######  Theorem 1.

Let $\\{x^{k}\\}_{k=0}^{\infty}$ be the sequence generated by Algorithm 1, and
let for $k\geq 0,\ \lambda_{k}\in(\varepsilon,2-\varepsilon)$ for some small
$\varepsilon>0$. Then for an arbitrary $x^{0}\in\mathbb{R}^{n}$ the sequence
$\\{x^{k}\\}_{k=0}^{\infty}$ converges to a fixed point of $U$ provided such a
fixed point exists.

###### Proof.

For $\lambda_{k}=\lambda$ see [6, Theorem 2.1]. For the non-stationary case
see [40, 49].

Note that the same interval for the relaxation parameters is valid both for
SIRT and P-SIRT.

######  Remark 1.

When $ker(A)=\emptyset$ then $U$ has a unique fixed point for any $C$. When
$ker(A)$ is not empty, e.g. when $A$ is underdetermined, then $U$ has a fixed
point when $C$ is compact, e.g. $C$ equals the unit hypercube (as in our
experiments). One may also note that with $C$ the rangespace of $A^{T}$ the
unique fixed point of $U$ is the pseudo-inverse solution.

###  1.4 Block-iterative Algorithms

There are two issues related to the practical use of SIRT-methods on large
data problems. One is the need for computer memory, and the other the often
slow rate of initial rate of convergence. By partitioning the matrix into row-
blocks and iterating sequentially over the blocks one can improve on both. A
careful study of implementation of block-iterative methods on a modern
multicore platform is presented in [57]. We mention also the possibility of
performing the iterations in parallel over the blocks [57, 47]. However, here
we do not consider this case. The use of block-iteration has been quite
successful in CT-applications, where it is often referred to as Ordered Subset
Iteration [8, 34]. A thorough convergence analysis of block-iteration used in
CT is [35]. A Landweber-Kaczmarz type block-iteration is proposed and analyzed
in [31]. Block-iterative versions of the examples listed in Section 1.2 are
also presented in the corresponding references.

Let the matrix $A$ be partitioned into $p$ blocks of equations, which may
contain common equations but each equation should appear at least in one of
the blocks. Choose $A_{t}\textrm{ and }b^{t}$ as the $t$th row-block of $A$
and $b,$ respectively. Let $\\{M_{t}\\}_{t=1}^{p}$ be a set of given SPD
weight matrices. Further $\\{C_{t}\\}_{t=1}^{p}$ is a family of closed convex
sets. Then (1) subject to $x\in C$ is equivalent to the following block
problems:

(5) $\text{find }x\in C_{t}\text{ satisfying }A_{t}x=b^{t}\text{,
}t=1,2,...,p\text{.}$

Defining $C=\bigcap_{t=1}^{p}C_{t}$ problem (5) can be also written as find
$x\in C$ satisfying $Ax=b$. In our study, we prefer, however, form (5),
because in algorithms for solving (5) we need only a sequential access to the
data. Recently, a more general form of this problem has been defined in [53,
Eq. (1.3)], see also [61, Algorithm 1].

Define cyclic control of the data by

(6) $i_{k}=[k]:=k(\mathop{\mathrm{m}od}p)+1.$

Put

(7) $\hat{A}_{[k]}=M_{[k]}^{1/2}A_{[k]},\ \hat{b}^{[k]}=M_{[k]}^{1/2}b^{[k]}.$

With these notations the Projected Block-Iterative Method (P-BIM) is

Algorithm 2 Projected Block Iterative Method (P-BIM)

1:$x^{0}\in\mathbb{R}^{n}$ is an arbitrary starting point

2:for $k\geq 0$ do
$x^{k+1}=U_{k}(x^{k})=P_{C_{[k]}}\left(x^{k}+\frac{\lambda_{k}}{\|\hat{A}_{[k]}\|^{2}}\hat{A}_{[k]}^{T}\left(\hat{b}^{[k]}-\hat{A}_{[k]}x^{k}\right)\right)$

3:end for

Note that the $k-$dependence is in the relaxation parameters, and the cyclic
access of the data $\\{C_{t},A_{t},M_{t},b^{t}\\},\ t=1,2,\dots,p$. A special
case of the above algorithm has been studied in [48, Eq. (7)], where
$C_{t}=\mathbb{R}^{n}$, $t=1,2,\cdots,p$. A parallel gradient based projection
algorithm has been studied in [18, Algorithm 4] and [53]. The difference with
respect to our algorithm is that they use the full gradient, whereas in our
approach only a part of the gradient (of the underlying least squares
functional) is used, which gives rise to the block structure. In the
optimization literature such methods are sometimes called incremental methods.
The method $(1.4)$ and $(1.5)$ in [45] with $f=\|Ax-b\|^{2}_{M}$, and assuming
the relaxation parameter is constant during a whole cycle is identical to
P-BIM. In [45] several relaxation strategies are investigated. The one most
similar to our choice is using a constant parameter $\lambda_{k}=\lambda$.
Here however, the convergence result is rather weak see [45, Proposition 2.1].

We next characterize the fixed points of the operator $U_{k}$ in P-BIM.
Similarly as in (4) we get

(8) $\mathop{\mathrm{F}ix}U_{k}=\arg\min_{x\in
C_{[k]}}\|b^{[k]}-A_{[k]}x\|^{2}_{M_{[k]}}.$

Put

$\hat{U}_{t}(x)=P_{C_{t}}\left(x+\frac{\lambda}{\|\hat{A}_{t}\|^{2}}\hat{A}_{t}^{T}\left(\hat{b}^{t}-\hat{A}_{t}x\right)\right),~{}t=1,2,\dots,p.$

Then

$\mathop{\mathrm{F}ix}\hat{U}_{t}=\arg\min_{x\in
C_{t}}\|b^{t}-A_{t}x\|^{2}_{M_{t}},\ t=1,2,\dots,p.$

It holds

$\bigcap_{k=0}^{\infty}\mathop{\mathrm{F}ix}U_{k}=\bigcap_{t=1}^{p}\mathop{\mathrm{F}ix}\hat{U}_{t}.$

The following result will be shown in Section 2 in Hilbert space settings.

######  Theorem 2.

Let $\varepsilon\in(0,1)$. Assume that
$\lambda_{k}\in[\varepsilon,2-\varepsilon]$. If
$\bigcap_{t=1}^{p}\mathop{\mathrm{F}ix}\hat{U}_{t}\neq\emptyset$, then the
sequence $\\{x^{k}\\}_{k=0}^{\infty}$ produced by the Algorithm P-BIM
converges to $x^{*}\in\bigcap_{t=1}^{p}\mathop{\mathrm{F}ix}\hat{U}_{t}$.
Further if $\\{C_{t}\\}$ is polyhedral the iteration converges linearly.

###  1.5 Noise-error

A common situation is when the right-hand side $b$ is corrupted by noise, i.e.
$b=\bar{b}+\delta b$, where $\bar{b}$ denotes the exact right-hand side. Let
$\bar{x}^{k}$ be the iteration vector using the exact right-hand side
$\bar{b}$ and let $x^{*}$ be a fixed point of the unperturbed iteration. The
total error can be decomposed into two terms

(9) $x^{k}-x^{*}=(x^{k}-\bar{x}^{k})+(\bar{x}^{k}-x^{*}).$

The first term is called the _noise error_ , and the second the _iteration
error_. During the first iterations of a convergent method the iteration error
dominates, and hence the total error decreases – but after a while the noise
error starts to grow resulting in so-called semi-convergence, as will be
demonstrated in Section 5. This phenomenon originally defined in [43, Section
IV.1] is further studied in [43, 25, 5, 33]. It has been observed
experimentally that the minimum error, i.e. $\min_{k\geq 0}\|x^{k}-x^{*}\|$,
is almost independent of both relaxation parameters and the weight matrix.
However the rate of decrease of the iteration error and the rate of increase
of the noise error are affected by both. The noise-error measures the
deviation due to errors in the data. Here we will only consider errors in the
rhs $b$. Possible errors in $A$ are usually caused by modelling errors, and
should then be dealt with using model refinements. We will study how the
noise-error evolves as a function of the iteration index $k$. Rather trivially
it is proportional to $k$. However we will derive a bound proportional to
$\sqrt{k}$, see (41). We next describe some earlier contributions. In [24], an
upper bound for the noise error of SIRT with a fixed relaxation parameter was
derived. Based on minimizing this upper bound, two strategies for choosing
relaxation parameters were suggested. The obtained sequences of relaxation
parameters are both non-negative and non-ascending. Later in [21] these
strategies were analyzed and used in P-SIRT. An upper bound for the noise
error of the block version of SIRT was derived in [48], and the above two
strategies for picking relaxation parameters together with a third choice were
studied. Recently, the stationary case of SIRT and P-SIRT were studied in
[22], [38] when the M-matrix is non-symmetric. This case includes the
classical Kaczmarz method [36].

###  1.6 Organization

In Section 2 we define and study a sequential block-iterative iteration (of
which P-BIM is a special case) adapted for solving split feasibility problems
in Hilbert space. The proof of Theorem 2 then follows from a more general
result in Section 2. We further study the noise error of P-BIM in Section 3,
and derive a new upper bound. In particular we extend the noise error analysis
of block-iteration to projected block-iteration, and also to the case when the
relaxation parameters are allowed to depend on the noise. In Section 4 we
propose a new rule for picking relaxation parameters, and present its
properties. We demonstrate the performance of P-BIM using this and other
relaxation strategies with examples taken from tomographic imaging in Section
5.

##  2 Hilbert space analysis

In this section we will first formulate a quite general feasibility problem.
Then we continue to study algorithms, and their convergence properties for the
problem. It will be shown that P-BIM comes out as a special case.

Let $\mathcal{H}$ be a real Hilbert and
$C_{i}\subseteq\mathcal{H},i=1,2,...,m,$ be closed convex subsets. The convex
feasibility problem (CFP) is to find $x\in\bigcap_{i=1}^{m}C_{i}$ if such a
point exists. CFP is a classical problem and was studied within the last
decades by many researchers, see, e.g. [3, 9, 4] and the references therein,
where several methods for solving the CFP and their convergence properties
were presented. First we recall some notions regarding operators in a Hilbert
space. Let $X\subseteq\mathcal{H}$ be a nonempty convex set and let
$T:X\rightarrow\mathcal{H}$ be an operator with a fixed point. Then $T$ is
called $\alpha$-strongly quasi-nonexpansive ($\alpha-$SQNE), where $\alpha\geq
0$, if

$\|Tx-z\|^{2}\leq\|x-z\|^{2}-\alpha\|Tx-x\|^{2},\ \forall x\in X,\text{ and
}z\in\mathop{\mathrm{F}ix}T,$

[9, Definition 2.1.38]. The operator $T$ is called a cutter [9, Definition
2.1.30] if for all $x\in X$ and for all $z\in\mathop{\mathrm{F}ix}T$ it holds

$\langle x-Tx,z-Tx\rangle\leq 0\text{.}$

By [27, Section 11], a firmly nonexpansive operator having a fixed point is a
cutter, see e.g., [9, Section 2.2] or [4, Section 4.1] for a definition of a
firmly nonexpansive operator.

Further, $T$ is called weakly regular if $\mathop{\mathrm{I}d}-T$ is demi-
closed at $0$, i.e., $x_{k}\rightharpoonup x$ and $x^{k}-Tx_{k}\rightarrow 0$
yields $x\in\mathop{\mathrm{F}ix}T$. Here $\rightharpoonup$ is a symbol for
weak convergence. $T$ is called linearly regular if there is a constant
$\delta>0$ such that for all $x\in X$ it holds

$d(x,\mathop{\mathrm{F}ix}T)\leq\delta\|Tx-x\|\text{,}$

cf. [11, Definition 3.1], where $d(x,C):=\inf_{y\in C}\|x-y\|$ denotes the
distance of $x$ to a subset $C$. In a similar way one can define weak and
linear regularity of sequences of operators, see [11, Definition 4.1] for
details.

###  2.1 Split feasibility problem

Let $\mathcal{H},\mathcal{H}_{0}$ be real Hilbert spaces,
$C\subseteq\mathcal{H}$ and $Q\subseteq\mathcal{H}_{0}$ be closed convex
subsets and $A:\mathcal{H}\rightarrow\mathcal{H}_{0}$ be a nonzero bounded
linear operator. The split feasibility problem (SFP) is to

$\text{find }x\in C\text{ such that }Ax\in Q\text{.}$

A more general form of the SFP is to

(10) $\text{minimize}~{}\|P_{Q}(Ax)-Ax\|^{2}~{}\text{subject to}~{}x\in C,$

or, equivalently, to find a fixed point of the operator
$T:\mathcal{H}\rightarrow\mathcal{H}$ defined by
$T(x)=P_{C}(x+\frac{1}{\|A\|^{2}}A^{\ast}(P_{Q}(Ax)-Ax)$. The SFP was
introduced by Censor and Elfving [13] in the finite dimensional settings.
Byrne [6] proved the weak convergence of the so called $CQ$-method,
$x^{k+1}=T_{\lambda}(x^{k})$, to $\mathop{\mathrm{a}rgmin}_{x\in
C}\|P_{Q}(Ax)-Ax\|$ which is a solution of the SFP. Here $\lambda\in(0,2)$ is
a relaxation parameter and $\|A\|$ denotes the spectral norm of $A$. In the
last $20$ years the SFP and its variants was intensively studied by many
researchers, see, e.g. [41, 60, 59, 10, 12, 54] to mention some of them.

###  2.2 Block SFP

Let $\mathcal{H},\mathcal{H}_{i}$, $i\in I:=\\{1,2,\dots,p\\}$ be Hilbert
spaces. Consider a family of problems:

(11) $\text{find }x\in C_{i}\text{ such that }A_{i}x\in Q_{i}\text{,}$

$i\in I$, where $C_{i}\subseteq\mathcal{H},Q_{i}\subseteq\mathcal{H}_{i}$ are
nonempty closed convex subsets, and
$A_{i}:\mathcal{H}\rightarrow\mathcal{H}_{i}$ are nonzero bounded linear
operators, $i\in I$. We call this family of problems, the block SFP. Denote by

(12) $F_{i}:=C_{i}\cap A_{i}^{-1}(Q_{i})\text{,}$

the solution set of problem (11), $i\in I$. Define

$\mathcal{H}_{0}:=\mathcal{H}_{1}\times\mathcal{H}_{2}\times\dots\times\mathcal{H}_{p},$

and a bounded linear operator $A:\mathcal{H}\rightarrow\mathcal{H}_{0}$ by

(13) $Ax=(A_{1}x,A_{2}x,\dots,A_{p}x).$

Put

(14) $Q=Q_{1}\times Q_{2}\times,\dots\times Q_{p},\ C=C_{1}\cap
C_{2}\cap\dots\cap C_{p}\text{,}$

and suppose that there is a common solution of all problems (11), $i\in I$,
i.e.,

(15) $\Omega:=\bigcap_{i=1}^{p}F_{i}=C\cap A^{-1}(Q)\neq\emptyset.$

The family of problems (11), $i\in I$, can be written as the following split
feasibility problem [6, 9, 15]

(16) $\text{find }x\in C\text{ such that }Ax\in Q\text{.}$

This problem formulation includes problem (5) by setting $Q_{i}=\\{b^{i}\\}$,
$i\in I$. By picking $Q=\\{b\\}$ (and using the M-norm) problem (10) becomes
minimization of $\|Ax-b\|_{M}$ subject to $x\in C$, i.e. problem (4).
Similarly when $p>1$ we retrieve the corresponding block-problem as expressed
in Theorem 2. In fact one can reformulate the block SFP into an optimization
problem (as for the case $p=1$, and problem (10)) see [15] for details.
Several extensions of (16) have been proposed, [17, 7]. In, e.g. [17] the sets
$C,Q$ were extended to fixed points sets of certain operators, $S$ and $T$.
With $S=P_{C}$, $T=P_{Q}$ the original formulation is retrieved. An
interesting extension is compressed sensing, where a sparse solution of a
consistent linear system is sought. Here one uses $l_{1}$-regularization which
yields the $Q$-lasso problem $\min_{x\in
C}\|\left(I-P_{Q}\right)Ax\|^{2}_{2}+\gamma\|x\|_{1}$. With $C=\mathbb{R}^{n}$
and $Q=\\{b\\}$ the original lasso problem [58] is recovered. For some further
work along these lines see [42].

In order to prove Theorem 2, we present below some general convergence results
which follow from [11]. Let $S_{i}:\mathcal{H}_{i}\rightarrow\mathcal{H}_{i}$
be a cutter with $\mathop{\mathrm{F}ix}S_{i}=Q_{i}$, $i\in I$. Define an
operator $T_{i}:\mathcal{H}\rightarrow\mathcal{H}$ by

(17)
$T_{i}(x)=\mathcal{L}\\{S_{i}\\}(x):=x+\frac{1}{\|A_{i}\|^{2}}A_{i}^{\ast}(S_{i}(A_{i}x)-A_{i}x)\text{,}$

i.e., $T_{i}$ is the Landweber operator corresponding to $S_{i}$, $i\in I$
(cf. [12, Definition 4.1]). Then the $\lambda$-relaxation of $T_{i}$ is
$T_{i,\lambda}:=(1-\lambda)\mathop{\mathrm{I}d}+\lambda T_{i}$. Hence

(18)
$T_{i,\lambda}(x)=\mathcal{L}\\{S_{i,\lambda}\\}(x):=x+\frac{\lambda}{\|A_{i}\|^{2}}A_{i}^{\ast}(S_{i}(A_{i}x)-A_{i}x)\text{,}$

where $\lambda\in(0,2)$, $i\in I$. By [59, Lemma 3.1] and by (15), $T_{i}$ is
a cutter and
$\mathop{\mathrm{F}ix}T_{i,\lambda}=\mathop{\mathrm{F}ix}T_{i}=A_{i}^{-1}(Q_{i})\neq\emptyset$.
Let $U_{i}$ be a cutter with $\mathop{\mathrm{F}ix}U_{i}=C_{i}$, $i\in I$, and
$U_{i,\mu}$ be the $\mu$-relaxation of $U_{i}.$ Define

(19) $V_{i}^{\lambda,\mu}:=U_{i,\mu}T_{i,\lambda}\text{,}$

$i\in I$, where $\lambda,\mu\in(0,2)$. Let
$\\{i_{k}\\}_{k=0}^{\infty}\subseteq\\{1,2,...,p\\}$ be a control sequence and
$x^{k}$ be generated by the iteration

(20)
$x^{k+1}=U_{i_{k},\mu_{k}}(x^{k}+\frac{\lambda_{k}}{\|A_{i_{k}}\|^{2}}A_{i_{k}}^{\ast}((S_{i_{k}}(A_{i_{k}}x^{k})-A_{i_{k}}x^{k})))$

or shorter

(21) $x^{k+1}=V_{k}x^{k}\text{,}$

where $V_{k}=V_{i_{k}}^{\lambda_{k},\mu_{k}}$,
$\lambda_{k},\mu_{k}\in[\varepsilon,2-\varepsilon]$, $k\geq 0$, for some
$\varepsilon\in(0,1)$ and $x^{0}\in\mathcal{H}$ is arbitrary. Recall that a
sequence $\\{i_{k}\\}_{k=0}^{\infty}\subseteq I=\\{1,2,...,p\\}$ is called an
almost cyclic control if there is a constant $s\geq p$ such that for any
$k\geq 0$ the set $\\{i_{k+1,}i_{k+2},...,i_{k+s}\\}\subseteq I$ (cf. [9,
Definition 5.6.10]). Theorems 3 and 5 presented below apply results of [11].

######  Theorem 3.

Let $S_{i}$ and $U_{i}$ be weakly regular cutters, $i\in I$,
$\lambda_{k},\mu_{k}\in[\varepsilon,2-\varepsilon]$, $k\geq 0$, for some
$\varepsilon\in(0,1)$ and $\\{i_{k}\\}_{k=0}^{\infty}$ be an almost cyclic
control. Then the sequence $x^{k}$ generated by iteration (21) converges
weakly to some $x^{\ast}\in\Omega$.

###### Proof.

Let $\lambda,\mu\in(0,2)$. Because $T_{i}$ is a cutter, its
$\lambda$-relaxation $T_{i,\lambda}$ defined by (18) is $\alpha$-SQNE, $i\in
I$, where $\alpha=(2-\lambda)/\lambda$ (see, e.g., [9, Theorem 2.1.39]) and
$\mathop{\mathrm{F}ix}T_{i,\lambda}=A_{i}^{-1}(Q_{i})\neq\emptyset$, because
$A_{i}^{-1}(Q_{i})\neq\emptyset$, $i=1,2,...,p$. Similarly, $U_{i,\mu}$ is
$\beta$-SQNE, where $\beta=(2-\mu)/\mu$. By $C_{i}\cap
A_{i}^{-1}(Q_{i})\neq\emptyset$ (see (15)), the operator $V_{i}^{\lambda,\mu}$
is $\rho$-SQNE, where $\rho=(\alpha^{-1}+\beta^{-1})^{-1}$ and
$\mathop{\mathrm{F}ix}V_{i}^{\lambda,\mu}=C_{i}\cap A_{i}^{-1}(Q_{i})$ [9,
Corollary 2.1.47 and Theorem 2.1.26(ii)], $i\in I$. Thus, one can easily check
that $V_{i_{k}}^{\lambda_{k},\mu_{k}}$ is $\rho_{k}$-SQNE, where
$\rho_{k}\geq\varepsilon/(4-2\varepsilon)>0$. Let $n_{-1}^{i}=-1$ and

(22) $n_{k}^{i}=\min\\{n>n_{k-1}^{i}:i_{n}=i\\}\text{,}$

$i\in I$, $k\geq 0$. Because $\\{i_{k}\\}_{k=0}^{\infty}$ is almost cyclic,
$n_{k}^{i}$ is well defined and there is $s\geq p$ such that
$n_{k+1}^{i}-n_{k}^{i}\leq s$ for all $k\geq 0$, $i=1,2,...,p$. The sequences
$\\{T_{i,\lambda_{n_{k}^{i}}}\\}_{k=0}^{\infty}$ and
$\\{U_{i,\mu_{n_{k}^{i}}}\\}_{k=0}^{\infty}$ are weakly regular [11,
Proposition 4.7]. By
$V_{n_{k}^{i}}=V_{i}^{\lambda_{n_{k}^{i}},\mu_{n_{k}^{i}}}=U_{i,\mu_{n_{k}^{i}}}T_{i,\lambda_{n_{k}^{i}}}$,
the sequence $\\{V_{n_{k}^{i}}\\}_{k=0}^{\infty}$ is weakly regular [11,
Corollary 5.5]. Now [11, Theorem 6.2(i)] yields the weak convergence of
$x^{k}$ to $x^{\ast}\in\Omega$.

Setting $S_{i}=P_{Q_{i}}$, $x\in\mathcal{H}$, $U_{i}=P_{C_{i}}$, $\mu_{k}=1$
and $\\{i_{k}\\}_{k=0}^{\infty}$ a cyclic control, i.e.,
$i_{k}=[k]:=k(\mathop{\mathrm{m}od}p)+1$, we obtain the following.

######  Corollary 4.

Let $\lambda_{k}\in[\varepsilon,2-\varepsilon]$, $k\geq 0$, for some
$\varepsilon\in(0,1)$. Then the sequence $x^{k}$ generated by iteration

(23)
$x^{k+1}=P_{C_{[k]}}(x^{k}+\frac{\lambda_{k}}{\|A_{[k]}\|^{2}}A_{[k]}^{\ast}((P_{Q_{[k]}}(A_{[k]}x^{k})-A_{[k]}x^{k})).$

converges weakly to some $x^{\ast}\in\Omega$.

###### Proof.

The metric projection, as a nonexpansive operator is weakly regular [50, Lemma
2]. Thus, the Corollary follows directly from Theorem 3.

Now taking $Q_{[k]}=\\{\hat{b}^{[k]}\\}$ and $A_{[k]}=\hat{A}_{[k]}$ in (23)
and assuming $\mathcal{H}_{i}$ being of finite dimension we retrieve P-BIM and
this proves Theorem 2.

Next suppose that $C_{i}$ is a polyhedral set (the intersection of a finite
family of half spaces), $U_{i}=P_{C_{i}}$, $\mathcal{H}_{i}$ is finite
dimensional, $Q_{i}=\\{b^{i}\\}$ and $S_{i}=P_{Q_{i}}$, $i\in I$. Then
$P_{Q_{i}}(y)=b^{i}$, $y\in\mathcal{H}_{i}$, the range of $A_{i}$ is closed,
$A_{i}^{-1}(Q_{i})$ is polyhedral, $i\in I$, and iteration (20) can be written
as

(24)
$x^{k+1}=P_{C_{i_{k}},\mu_{k}}(x^{k}+\frac{\lambda_{k}}{\|A_{i_{k}}\|^{2}}A_{i_{k}}^{\ast}(b_{i_{k}}-A_{i_{k}}x^{k}))\text{.}$

Denote $b=(b^{1},b^{2},...b^{p})\in\mathcal{H}_{0}$. Then
$A^{-1}(Q)=\\{x\in\mathcal{H}:Ax=b\\}$ and the solution set has the form
$\Omega=\\{x\in C:Ax=b\\}$.

######  Theorem 5.

Let $C_{i}$, $i\in I$, be polyhedral,
$\lambda_{k},\mu_{k}\in[\varepsilon,2-\varepsilon]$, $k\geq 0$, for some
$\varepsilon\in(0,1)$ and $\\{i_{k}\\}_{k=0}^{\infty}\subseteq I$ be almost
cyclic. Then the sequence $x^{k}$ generated by iteration (24) converges
linearly to some $x^{\ast}\in\Omega$.

###### Proof.

Similarly as in the proof of Theorem 3, $V_{i_{k}}^{\lambda_{k},\mu_{k}}$ is
$\rho_{k}$-SQNE, where $\rho_{k}\geq\varepsilon/(4-2\varepsilon)>0$. Moreover,
the sequences $\\{T_{i,\lambda_{n_{k}^{i}}}\\}_{k=0}^{\infty}$ and
$\\{U_{i,\mu_{n_{k}^{i}}}\\}_{k=0}^{\infty}$ are linearly regular [11,
Proposition 4.7]. Because, $C_{i}\ $and $A_{i}^{-1}(Q_{i})$ are polyhedral
sets, the family $\\{C_{i},A_{i}^{-1}(Q_{i})\\}$ is linearly regular, $i\in I$
[3, Corollary 5.26]. Thus, the sequences $\\{V_{n_{k}^{i}}\\}_{k=0}^{\infty}$,
$i\in I$, are linearly regular [11, Corollary 5.5(iii)] (for a definition of a
linearly regular sequence of operators, see [11, Definition 4.1]). Because
$F_{i}:=C_{i}\cap A_{i}^{-1}(Q_{i})$, $i\in I$, are polyhedral sets, the
family $\\{F_{i}:i\in I\\}$ is linearly regular [3, Corollary 5.26]. Now [11,
Theorem 6.2(iii)] yields the linear convergence of $x^{k}$ to
$x^{\ast}\in\Omega$.

We end this section by considering the case $N\neq\textrm{Id}$. For ease of
notation we only consider the simultaneous case. Let
$A:\mathcal{H}_{1}\rightarrow\mathcal{H}_{2}$ be a bounded linear operator and
$b\in\mathcal{H}_{2}$. Define
$T_{A,b}:\mathcal{H}_{1}\rightarrow\mathcal{H}_{1}$ by

(25) $T_{A,b}(x)=x+\frac{1}{\|A\|^{2}}A^{\ast}(b-Ax)\text{,}$

$x\in\mathcal{H}_{1}$. The operator $T_{A,b}$ is firmly nonexpansive ([9,
Definition 2.2.1]), and
$\mathop{\mathrm{F}ix}T_{A,b}=\\{x\in\mathcal{H}_{1}:A^{\ast}(b-Ax)=0\\}$,
see, e.g., [9, Lemma 4.6.2 and Theorem 4.6.3]. In particular, if
$A^{-1}(\\{b\\})\neq\emptyset$ then
$\mathop{\mathrm{F}ix}T_{A,b}=A^{-1}(\\{b\\})$ [59, Lemma 3.1].

Let $N:\mathcal{H}_{1}\rightarrow\mathcal{H}_{1}$ and
$M:\mathcal{H}_{2}\rightarrow\mathcal{H}_{2}$ be strongly positive symmetric
operators and $N^{\frac{1}{2}}$ and $M^{\frac{1}{2}}$ their square roots. It
is well known that $N^{\frac{1}{2}}$ and $M^{\frac{1}{2}}$ are also strongly
positive definite symmetric operators [62, page 618]. Let
$\tilde{A}:=M^{\frac{1}{2}}AN^{\frac{1}{2}}$ and $\hat{b}:=M^{\frac{1}{2}}b$.
Define $\hat{T}:=T_{\tilde{A},\hat{b}\,}$, i.e.,

$\hat{T}(x)=x+\frac{1}{\|\tilde{A}\|^{2}}\tilde{A}^{\ast}(\hat{b}-\tilde{A}x)=x+\frac{1}{\|M^{\frac{1}{2}}AN^{\frac{1}{2}}\|^{2}}N^{\frac{1}{2}}A^{\ast}M(b-AN^{\frac{1}{2}}x)\text{.}$

Clearly, $\hat{T}$ is firmly nonexpansive and

$\mathop{\mathrm{F}ix}\hat{T}=\\{x\in\mathcal{H}_{1}:\hat{A}^{\ast}(\hat{b}-\hat{A}x)=0\\}=\\{x\in\mathcal{H}_{1}:N^{\frac{1}{2}}A^{\ast}M(b-AN^{\frac{1}{2}}x)=0\\}\text{.}$

Let $\hat{T}_{\lambda}$ denote the $\lambda$-relaxation of $\hat{T}$, where
$\lambda\in[0,2]$, and recall the definition of an averaged operator ([9,
Definition 2.2.16]).

######  Proposition 6.

The operator $V:\mathcal{H}_{1}\rightarrow\mathcal{H}_{1}$ defined by
$V=P_{C}\hat{T}_{\lambda}$, where $\lambda\in(0,2)$, is
$\frac{2}{4-\lambda}$-averaged.

###### Proof.

$V$ is $\frac{4}{4-\lambda}$-relaxed firmly nonexpansive being a composition
of a firmly nonexpansive operator $P_{C}$ and a $\lambda$-relaxed firmly
nonexpansive operator $\hat{T}_{\lambda}$ [9, Theorem 2.2.37]. Consequently,
$V$ is $\frac{2}{4-\lambda}$-averaged [9, Corollary 2.2.17].

Consider now the iterates generated by $x^{k+1}=Vx^{k}.$ It follows that
$x^{k}\rightharpoonup x^{*}\in\mathop{\mathrm{F}ix}V$ provided
$\mathop{\mathrm{F}ix}V\neq\emptyset.$

##  3 P-BIM and its noise error

In this section we will study how errors in the unperturbed right-hand-side
$\bar{b}$ affect the iterates in P-BIM. To this end let

(26) $b=\bar{b}+\delta b.$

The noise will apart from the iterates also possibly affect the relaxation
parameters (depending on their definition). To cope with possible rank-
deficiency we will, similarly as in [21, 48] add a regularization term
$\alpha\|x\|^{2}$. Let

(27) $u^{k}(x,b^{[k]})=A_{[k]}^{T}M_{[k]}\left(b^{[k]}-A_{[k]}x\right)-\alpha
x,$

where the last term is due to regularization. Next we introduce the notations

(28) $\theta_{k}=\frac{\lambda_{k}}{\|\hat{A}_{[k]}\|^{2}+\alpha},\
\bar{\theta}_{k}=\frac{\bar{\lambda}_{k}}{\|\hat{A}_{[k]}\|^{2}+\alpha}.$

We will in the sequel refer to $\theta_{k}$ and $\bar{\theta}_{k}$ as the
noisy and noise-free relaxation parameter, respectively. Note that these are
just a scaled version of the original relaxation parameters $\lambda_{k}$ and
$\bar{\lambda}_{k}$. For noise-free data $\lambda_{k}=\bar{\lambda}_{k}$. Now
the noise-free P-BIM becomes

(29)
$\bar{x}^{k+1}=P_{C_{[k]}}\left(\bar{x}^{k}+\bar{\theta}_{k}u^{k}(\bar{x}^{k},\bar{b}^{[k]})\right).$

The noisy version of P-BIM is

(30) $x^{k+1}=P_{C_{[k]}}\left(x^{k}+\theta_{k}u^{k}(x^{k},b^{[k]})\right).$

Next define

(31) $\hat{\delta}=\max_{1\leq t\leq p}\|A_{t}^{T}M_{t}\delta b^{t}\|,$ (32)
$\gamma_{k}=\max_{0\leq j\leq k}|\bar{\theta}_{j}-\theta_{j}|,\
\underline{\eta}_{k}=\min_{0\leq j\leq k}\theta_{j},\
\overline{\eta}_{k}=\max_{0\leq j\leq k}\theta_{j}.$

Let $\underline{\sigma}_{t}$ be the smallest nonzero singular value of
$M_{t}^{1/2}A_{t}$, and define

(33) $\underline{\sigma}=\min_{1\leq t\leq
p}\underline{\sigma}_{t},\quad\overline{\sigma}=\max_{1\leq t\leq
p}\|M_{t}^{1/2}A_{t}\|,$

and

(34) $\hat{u}_{k}=\max_{0\leq s\leq k}\|u^{s}(\bar{x}^{s},\bar{b}^{[s]})\|.$

######  Theorem 7.

Assume that $\alpha=\underline{\sigma}^{2}$. Then the noise-error in P-BIM is
bounded above by

(35)
$e^{k}_{N}:=\|x^{k}-\bar{x}^{k}\|\leq\frac{1}{\underline{\sigma}}\left(\frac{\gamma_{k-1}\hat{u}_{k-1}+\overline{\eta}_{k-1}\hat{\delta}}{\underline{\eta}_{k-1}}\right)\Psi^{k}(\underline{\sigma},\underline{\eta}_{k-1}),$

where $\Psi^{k}(x,y)$ is defined by

(36) $\Psi^{k}(x,y)\equiv\frac{1-(1-yx^{2})^{k}}{x}.$

###### Proof.

Since $P_{C}$ is nonexpansive we get

$\displaystyle e^{k}_{N}$
$\displaystyle\leq\|x^{k-1}+\theta_{k-1}u^{k-1}(x^{k-1},b^{[k]})-\left(\bar{x}^{k-1}+\bar{\theta}_{k-1}u^{k-1}(\bar{x}^{k-1},\bar{b}^{[k]})\right)\|$
$\displaystyle=\|(x^{k-1}-\bar{x}^{k-1})+\theta_{k-1}(u^{k-1}(x^{k-1},b^{[k]})-u^{k-1}(\bar{x}^{k-1},\bar{b}^{[k]}))$
$\displaystyle+(\theta_{k-1}-\bar{\theta}_{k-1})u^{k-1}(\bar{x}^{k-1},\bar{b}^{[k]})\|.$

Now

$\displaystyle u^{k-1}(x^{k-1},b^{[k]})-u^{k-1}(\bar{x}^{k-1},\bar{b}^{[k]})$
$\displaystyle=A_{[k-1]}^{T}M_{[k-1]}\left(b^{[k-1]}-A_{[k-1]}x^{k-1}\right)$
$\displaystyle\quad-
A_{[k-1]}^{T}M_{[k-1]}\left(\bar{b}^{[k-1]}-A_{[k-1]}\bar{x}^{k-1}\right)$
$\displaystyle\quad+\alpha(\bar{x}^{k-1}-x^{k-1})$
$\displaystyle=A_{[k-1]}^{T}M_{[k-1]}\delta b^{[k-1]}$
$\displaystyle\quad+(\alpha\textrm{Id}+A_{[k-1]}^{T}M_{[k-1]}A_{[k-1]})(\bar{x}^{k-1}-x^{k-1}).$

Hence

$\displaystyle e^{k}_{N}$ $\displaystyle\leq$
$\displaystyle\|\left((1-\theta_{k-1}\alpha)\textrm{Id}-\theta_{k-1}A_{[k-1]}^{T}M_{[k-1]}A_{[k-1]}\right)(x^{k-1}-\bar{x}^{k-1})$
$\displaystyle+$ $\displaystyle\theta_{k-1}A_{[k-1]}^{T}M_{[k-1]}\delta
b^{[k-1]}+(\theta_{k-1}-\bar{\theta}_{k-1})u^{k-1}(\bar{x}^{k-1},\bar{b}^{[k-1]})\|.$

Here

(37)
$Q_{k}=\left((1-\theta_{k}\alpha)\textrm{Id}-\theta_{k}A_{[k]}^{T}M_{[k]}A_{[k]}\right),\quad
q_{k}=\|Q_{k}\|.$

Then it follows

$e^{k}_{N}\leq
q_{k-1}e_{N}^{k-1}+|\theta_{k-1}-\bar{\theta}_{k-1}|.\|u^{k-1}(\bar{x}^{k-1},\bar{b}^{[k-1]})\|+\theta_{k-1}\hat{\delta}.$

Assuming $e^{0}_{N}=x^{0}-\bar{x}^{0}=0$ it follows by induction

$\displaystyle e^{k}_{N}$
$\displaystyle\leq\sum^{k-2}_{s=0}|\theta_{s}-\bar{\theta}_{s}|.\|u^{s}(\bar{x}^{s},\bar{b}^{[s]})\|\prod_{j=s+1}^{k-1}q_{j}$
$\displaystyle+\hat{\delta}\sum^{k-2}_{s=0}\theta_{s}\prod_{j=s+1}^{k-1}q_{j}+|\theta_{k-1}-\bar{\theta}_{k-1}|.\|u^{k-1}(\bar{x}^{k-1},\bar{b}^{[k-1]})\|+\theta_{k-1}\hat{\delta}.$

With (34), and putting $\hat{q}_{k}=\max_{0\leq j\leq k}q_{j}$ we therefore
get

(38) $\displaystyle e^{k}_{N}$
$\displaystyle\leq(\gamma_{k-2}\hat{u}_{k-2}+\bar{\eta}_{k-2}\hat{\delta})\sum^{k-2}_{s=0}\prod_{j=s+1}^{k-1}\hat{q}_{k-1}+(\gamma_{k-1}\hat{u}_{k-1}+\bar{\eta}_{k-1}\hat{\delta}).$

Now

(39)
$\sum^{k-2}_{s=0}\prod_{j=s+1}^{k-1}\hat{q}_{k-1}=\left(\hat{q}_{k-1}^{k-1}+\hat{q}_{k-1}^{k-2}+\ldots+\hat{q}_{k-1}^{2}+\hat{q}_{k-1}^{1}\right).$

It follows from the properties of geometric progression that,

$e^{k}_{N}\leq\left(\gamma_{k-1}\hat{u}_{k-1}+\overline{\eta}_{k-1}\delta\right)\frac{1-\hat{q}_{k-1}^{k}}{1-\hat{q}_{k-1}}.$

We will now use the following result, derived in [48, Lemma 1], for block-
iteration, and [21, Lemma 3.9] for simultaneous iteration,

(40) $\|Q_{k}\|=1-\theta_{k}\underline{\sigma}^{2},\textrm{ assuming
}\alpha=\underline{\sigma}^{2}.$

Using (40) it holds $\hat{q}_{k-1}=\max_{0\leq j\leq
k-1}(1-\theta_{j}\underline{\sigma}^{2})=1-\underline{\eta}_{k-1}\underline{\sigma}^{2}.$
It follows

$e^{k}_{N}\leq\left(\gamma_{k-1}\hat{u}_{k-1}+\overline{\eta}_{k-1}\hat{\delta}\right)\frac{1-(1-\underline{\eta}_{k-1}\underline{\sigma}^{2})^{k}}{\underline{\eta}_{k-1}\underline{\sigma}^{2}}.$

Therefore the result follows using (36).

######  Remark 2.

In the proof of Theorem 7, we assumed cyclic control. However, provided the
same control sequence $\\{i_{k}\\}$ (where $1\leq i_{k}\leq p$) is used in
both the noise-free (29), and the noisy (30) iteration it is obvious that
Theorem 7 also holds. Examples of other controls are almost cyclic control as
defined above and stochastic control [39, Definition 5.1]. Stochastic control
has proven to be quite successful and we refer to [55, 44, 46] for some recent
developments.

######  Remark 3.

When the regularization parameter $\alpha$ is chosen positive the convergence
of P-BIM follows from the contraction mapping theorem as follows. By (40)
$Q_{k}$ is a $\rho_{k}$-contraction with
$\rho_{k}=1-\theta_{k}\underline{\sigma}^{2}.$ Since $P_{C}$ is a nonexpansive
it follows that $U_{k}$ (defined in Algorithm 2) also is a contraction.
Because $\rho_{k}<\rho<1,$ it follows that $U_{k}$ is a $\rho-$contraction.
Thus, assuming feasibility, $x^{k}$ converges to
$x^{*}\in\cap_{k=1}^{\infty}\mathop{\mathrm{F}ix}U_{k}$.

######  Remark 4.

If $\lambda_{k}=\bar{\lambda}_{k}$, i.e. the relaxation parameters do not
depend on the noise, then $\gamma_{k}=0$, and
$\underline{\eta}_{k}=\min_{0\leq j\leq
k}\theta_{j},\overline{\eta}_{k}=\max_{0\leq j\leq k}\theta_{j}$ so that

(41)
$\|x^{k}-\bar{x}^{k}\|\leq\frac{1}{\underline{\sigma}}\left(\frac{\overline{\eta}_{k-1}\hat{\delta}}{\underline{\eta}_{k-1}}\right)\Psi^{k}(\underline{\sigma},\underline{\eta}_{k-1})$

This bound, assuming decreasing relaxation parameters $\theta_{k}$, coincides
with the bound given in [48, Theorem 3].

Next we remind the reader of the inequality
$(1-(1-x^{2})^{k})/x\leq\sqrt{k},x\in(0,1)$. It follows

$\Psi^{k}(\underline{\sigma},\underline{\eta}_{k-1})=\frac{1-(1-\underline{\eta}_{k-1}\underline{\sigma}^{2})^{k}}{\sqrt{\underline{\eta}_{k-1}}\underline{\sigma}}\sqrt{\underline{\eta}_{k-1}}\leq\sqrt{\underline{\eta}_{k-1}}\sqrt{k},$

provided $\underline{\sigma}\leq 1/\underline{\eta}_{k-1}$. Hence we arrive at
the bound (with $c$ a small constant)

(42)
$\|x^{k}-\bar{x}^{k}\|\leq\frac{c}{\underline{\sigma}}\sqrt{k},~{}~{}~{}~{}~{}~{}\underline{\sigma}\leq
1/\underline{\eta}_{k-1},$

cf [22, 38] for Kaczmarz’s method and [25] for SIRT. Due to the factor
$1/\underline{\sigma}$ this bound is unrealistically large, but we have not
been able to sharpen it. See also [20] for a more detailed discussion.

We next give an alternative upper bound for the noise error presented in
Theorem 7, when the relaxation parameters are decreasing, i.e.
$0<\theta_{k+1}\leq\theta_{k}$. At first, following [24, Propositions 2.3,
2.4], we consider

(43) $g_{k-1}(y)=(2k-1)y^{k-1}-(y^{k-2}+\cdots+y+1)$

which has a unique root $\zeta_{k}\in(0,1)$ and satisfies
$0<\zeta_{k}<\zeta_{k+1}<1$ and $\lim_{k\to\infty}\zeta_{k}=1$.

######  Theorem 8.

If the relaxation parameters are decreasing and satisfy

(44) $0<\theta_{k}\leq\frac{1}{\|M_{t}^{1/2}A_{t}\|^{2}},$

then we have

(45)
$\|x^{k}-\bar{x}^{k}\|\leq\displaystyle\frac{1}{\underline{\sigma}}\left(\frac{(\theta_{0}-\theta_{k-1})\hat{u}_{k-1}+\theta_{0}\hat{\delta}}{{\sqrt{\theta_{k-1}}}}\right)\frac{1-{\zeta}^{k}_{k}}{\sqrt{1-\zeta_{k}}}.$

###### Proof.

It holds $\underline{\eta}_{k-1}=\min\theta_{s}=\theta_{k-1}$,
$\overline{\eta}_{k-1}=\max\theta_{s}=\theta_{0}$ and
$\gamma_{k-1}\leq|\theta_{0}-\theta_{k-1}|$. Based on [24, Proposition 2.4],
we have

(46) $\displaystyle\max_{0<\sigma\leq
1/{\sqrt{\theta_{k-1}}}}\Psi^{k}(\sigma,\theta_{k})\leq\sqrt{\theta_{k-1}}\frac{1-\zeta^{k}_{k}}{\sqrt{1-\zeta_{k}}}$

where $\zeta_{k}\in(0,1)$ is the unique root of (43). By applying (46) in (35)
we obtain (45).

##  4 Relaxation parameters

In [48] three step-size rules for use in block iterative methods were studied
(the first two were originally proposed in [24], and in [21] (constrained
case) for simultaneous iteration). The main idea with these rules is to
control the propagated noise component of the error. Here we use the
formulation (29) for P-BIM (to conform with the notations in the preceding
Section and in [24], [21]).

(47)
$\Psi_{1}-\textrm{rule}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\theta_{k}=\left\\{\begin{array}[]{ll}\displaystyle\sqrt{2}\overline{\sigma}^{-2},&~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\textrm{for}~{}~{}k=0,1\vspace{0.1cm}\\\
\displaystyle
2\overline{\sigma}^{-2}(1-\zeta_{k}),&~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\textrm{for}~{}~{}k\geq
2,\end{array}\right.$ (48)
$\Psi_{2}-\textrm{rule}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\theta_{k}=\left\\{\begin{array}[]{ll}\displaystyle\sqrt{2}\overline{\sigma}^{-2},&~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\textrm{for}~{}~{}k=0,1\vspace{0.1cm}\\\
\displaystyle
2\overline{\sigma}^{-2}\frac{1-\zeta_{k}}{(1-\zeta_{k}^{k})^{2}},&~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\textrm{for}~{}~{}k\geq
2,\end{array}\right.$ (49)
$\Psi_{3}-\textrm{rule}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\theta_{k}=\left\\{\begin{array}[]{ll}\displaystyle\sqrt{2}\overline{\sigma}^{-2},&~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\textrm{for}~{}~{}k=0,1\vspace{0.1cm}\\\
\displaystyle
2\overline{\sigma}^{-2}\frac{(1-\zeta_{k}^{k})^{2}}{(1-\zeta_{k})^{1-r}},&~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\textrm{for}~{}~{}k\geq
2.\end{array}\right.$

Here $1<r\leq 2.$ All three rules are descending, i.e.
$0<\theta_{k+1}<\theta_{k},\ k\geq 2$. Also for these three rules the
relaxation parameters do not depend on the noise.

The following upper bounds for the noise error were derived in [48] (for the
simultaneous case see also [21]).

(50)
$\left\|x^{k}-\bar{x}^{k}\right\|\leq\left\\{\begin{array}[]{ll}\displaystyle\frac{\beta_{\delta
b}}{\underline{\sigma}}\left(1-\zeta_{k}^{k}\right)/\left(1-\zeta_{k}\right),&~{}~{}~{}~{}~{}~{}~{}\Psi_{1}-\textrm{rule}\vspace{0.1cm}\\\
\displaystyle\frac{\beta_{\delta
b}}{\underline{\sigma}}\left(1-\zeta_{k}^{k}\right)^{2}/\left(1-\zeta_{k}\right),&~{}~{}~{}~{}~{}~{}~{}\Psi_{2}-\textrm{rule}\vspace{0.1cm}\\\
\displaystyle\frac{\beta_{\delta
b}}{\underline{\sigma}}\left(1-\zeta_{k}\right)^{\frac{-r}{2}}.&~{}~{}~{}~{}~{}~{}~{}\Psi_{3}-\textrm{rule}\vspace{0.1cm}\end{array}\right.$

Here

(51) $\beta_{\delta b}=\max_{1\leq t\leq p}\|M_{[t]}^{\frac{1}{2}}\delta
b^{[t]}\|,\ \beta_{b}=\max_{1\leq t\leq p}\|M_{[t]}^{\frac{1}{2}}b^{[t]}\|.$

We next derive a new relaxation parameter rule which will be based on the
bound (45). We will then use the following heuristic
$\hat{u}_{k}\leq\hat{u}_{0}.$ This was fulfilled in all our numerical
experiments (but we have no formal proof). Also we use, for simplicity
$x^{0}=0$ (as in the experiments). Hence

$\hat{u}_{k}\approx\hat{u}_{0}\leq\max_{1\leq t\leq
p}\|A_{t}M_{t}b^{t}\|\leq\overline{\sigma}\beta_{b}.$

Further using (31) $\hat{\delta}\leq\overline{\sigma}\beta_{\delta b}.$
Summarizing we get

(52) $\hat{u}_{k}\approx\bar{\sigma}\beta_{b},\
\hat{\delta}\approx\bar{\sigma}\beta_{\delta b}.$

Again the basic idea is to control the noise error. Therefore consider the
function

(53) $\frac{\beta_{\delta
b}}{\underline{\sigma}}\left(1-\zeta_{k}^{k}\right)^{\frac{1}{2}}\left(1-\zeta_{k}\right)^{\frac{-r}{2}}$

We now enforce that the noise error bound (45) equals (53). It follows, also
using (52)

$\frac{\overline{\sigma}}{\underline{\sigma}}\left(\frac{(\theta_{0}-\theta_{k-1})\beta_{b}+\theta_{0}\beta_{\delta
b}}{{\sqrt{\theta_{k-1}}}}\right)\frac{1-{\zeta}^{k}_{k}}{\sqrt{1-\zeta_{k}}}=\frac{\beta_{\delta
b}}{\underline{\sigma}}\left(1-\zeta_{k}^{k}\right)^{\frac{1}{2}}\left(1-\zeta_{k}\right)^{\frac{-r}{2}}$

After simplification we get

$\overline{\sigma}\left(\frac{(\theta_{0}-\theta_{k-1})\beta_{b}+\theta_{0}\beta_{\delta
b}}{\sqrt{\theta_{k-1}}}\right)=\beta_{\delta
b}\frac{(1-\zeta_{k})^{\frac{1-r}{2}}}{\sqrt{1-{\zeta}^{k}_{k}}}.$

By solving for $\theta_{k-1}$, and using that
$\theta_{0}=\sqrt{2}/\overline{\sigma}^{2}$ one gets (after elementary
calculations)

$\theta_{k-1}=\frac{\mathcal{B}+\mathcal{Z}_{k,r}^{2}\beta_{\delta
b}^{2}-\mathcal{Z}_{k,r}\beta_{\delta
b}\sqrt{\mathcal{Z}_{k,r}^{2}\beta_{\delta
b}^{2}+2\mathcal{B}}}{2\overline{\sigma}^{2}\beta_{b}^{2}},$

where

$\mathcal{B}=2\sqrt{2}\beta_{b}\left(\beta_{b}+\beta_{\delta b}\right),\
\mathcal{Z}_{k,r}=\frac{(1-\zeta_{k})^{\frac{1-r}{2}}}{\sqrt{(1-\zeta_{k}^{k})}}.$

Summarizing we have the following rule

(54)
$\Gamma-\textrm{rule}~{}~{}~{}~{}~{}~{}~{}~{}\theta_{k}=\left\\{\begin{array}[]{ll}\displaystyle\frac{\sqrt{2}}{\overline{\sigma}^{2}},&~{}~{}~{}~{}~{}~{}~{}\textrm{for}~{}~{}k=0,1\vspace{0.1cm}\\\
\displaystyle\frac{\mathcal{B}+\mathcal{Z}_{k,r}^{2}\beta_{\delta
b}^{2}-\mathcal{Z}_{k,r}\beta_{\delta
b}\sqrt{\mathcal{Z}_{k,r}^{2}\beta_{\delta
b}^{2}+2\mathcal{B}}}{2\overline{\sigma}^{2}\beta_{b}^{2}},&~{}~{}~{}~{}~{}~{}~{}\textrm{for}~{}~{}k\geq
2,\end{array}\right.$

where $1<r\leq 2$. Note that the parameters depend on the noise. With no noise
present $\theta_{k}=\theta_{0},\ k\geq 1$. We mention here the recent paper
[2, Proposition 5] where also knowledge of $\|\delta b\|$ is utilized. By
construction the noise upper bound for the $\Gamma-$rule is

(55) $\|x^{k}-\bar{x}^{k}\|\leq\frac{\beta_{\delta
b}}{\underline{\sigma}}\left(1-\zeta_{k}^{k}\right)^{\frac{1}{2}}\left(1-\zeta_{k}\right)^{\frac{-r}{2}},~{}~{}~{}~{}~{}~{}~{}\Gamma-\textrm{rule}.$

Let $||x^{k}-\bar{x}^{k}||_{\Omega}$ denote the upper bound for the noise
error when $\Omega$-rule is used for choosing relaxation parameters. Then the
following inequalities obviusly hold,

(56)
$||x^{k}-\bar{x}^{k}||_{\Gamma}\leq||x^{k}-\bar{x}^{k}||_{\Psi_{3}}\textrm{
and
}||x^{k}-\bar{x}^{k}||_{\Psi_{2}}\leq||x^{k}-\bar{x}^{k}||_{\Psi_{1}}~{}~{}~{}~{}~{}\textrm{for}~{}~{}~{}k\geq
0.$

The bounds for the $\Psi_{3}$ and the $\Psi_{2}$ rules respectively may
intersect see [48, Figure 2]. We next show that the $\Gamma-$rule also is a
diminishing step-size strategy,

######  Lemma 9.

Let $r\in(1,2],\ k\geq 2$. The relaxation parameters in the $\Gamma$-rule are
descending, i.e. $0<\theta_{k+1}<\theta_{k}$.

###### Proof.

Since $\zeta_{k}\in(0,1)$, we have

$\frac{d\mathcal{Z}_{k,r}}{d\zeta_{k}}=\left(\frac{r-1}{(1-\zeta_{k})}+\frac{k\zeta_{k}^{k-1}}{\left(1-\zeta_{k}^{k}\right)}\right)\frac{(1-\zeta_{k})^{\frac{1-r}{2}}}{2\sqrt{1-\zeta_{k}^{k}}}>0.$

Showing that $\mathcal{Z}_{k,r}(\zeta_{k})$ is increasing. Next

$\frac{d\mathcal{\theta}_{k}}{d\mathcal{Z}_{k,r}}=\frac{1}{2\overline{\sigma}^{2}\beta_{b}\sqrt{\beta_{\delta
b}^{2}\mathcal{Z}_{k,r}^{2}+2\mathcal{B}}}\tau_{k,r},$

with

$\tau_{k,r}=\beta_{\delta b}\mathcal{Z}_{k,r}\sqrt{\beta_{\delta
b}^{2}\mathcal{Z}_{k,r}^{2}+2\mathcal{B}}-\left(\mathcal{B}+\beta_{\delta
b}^{2}\mathcal{Z}_{k,r}^{2}\right).$

We now show that $\tau_{k,r}<0$. To this end assume the contrary, i.e.

$\beta_{\delta b}\mathcal{Z}_{k,r}\sqrt{\beta_{\delta
b}^{2}\mathcal{Z}_{k,r}^{2}+2\mathcal{B}}\geq\left(\mathcal{B}+\beta_{\delta
b}^{2}\mathcal{Z}_{k,r}^{2}\right),$

After squaring and simplifying one gets ${\mathcal{B}}^{2}\leq 0,$ a
contradiction. Hence $\frac{d\mathcal{\theta}_{k}}{d\mathcal{Z}_{k,r}}<0.$ As
noted previously $\zeta_{k}(k)$ is increasing. Therefore the result follows by
the chain rule.

In Section 5, in addition to the above rules also the ’$\theta-$opt’ rule is
used. This means finding that constant value of $\theta$ which give rise to
the smallest relative error within a fixed number of cycles ($cmax$). Here a
cycle denotes one pass through all $p$ row-blocks. This value of $\theta$ is
found by searching over the interval $0<\theta_{opt}<2/\overline{\sigma}^{2}$.
This strategy requires knowledge of the exact solution, so for real data one
would first need to train the algorithm using simulated data, see [32, Section
6]. In our tests we took $cmax=100$ and used the exact phantom for training.
This is clearly unrealistic but we think it is of interest to include this
rule for comparison with our other rules.

##  5 Numerical Results

We will report on tests using examples from the field of image reconstruction
from projections. To create the projection matrix $A$ and the right hand side
$b$ the paralleltomo function in the MATLAB package AIR tools [32] is used. We
take the Shepp-Logan phantom as the original image $x^{*}$ discretized into
$365\times 365$ square pixels. In the first test problem we use $88$ views
uniformly distributed over $180$ degrees, and $516$ projections per view.
Since zero rows do not contribute to the reconstruction, after identifying and
removing these, the resulting matrix $A$ has dimension $40796\times 133225$
(case-one). Although in our application iterative methods usually are more
competitive the more underdetermined the system is this might not be the case
in other applications. Therefore we also consider taking more rays (264) per
projection leading to a matrix of dimension $122388\times 133225$ (case-two).
Apart from using noise-free data we added independent Gaussian noise of mean
$0$ and relative noise-level ($\|\delta b\|/\|b\|$) $2\%$ and $5\%$
respectively. In the experiments we set $x^{0}=0$ and $C_{[t]}=[0,1]^{n}$ for
$t=1,\ldots,p$. We further partition the matrix $A$ and the right hand side
$b$ into $8$ and $22$ blocks, respectively. The largest and smallest (nonzero)
singular values of each block is estimated using the power method [56]. We
further used Cimmino’s $M$-Matrix. The error metric is Relative Error defined
by

(57) $\mbox{Relative Error}=\frac{\|x^{k}-x^{*}\|}{\|x^{*}\|}.$

All codes are written in MATLAB(R2015a) and conducted on a PC with a Intel
Core i7-7700K CPU @ 4.2 GHz and 16 GB RAM.

When using the $\Gamma-$rule we need to estimate $\|M^{\frac{1}{2}}\delta
b\|.$ Therefore we randomly generated a new vector $\overline{\delta b}$
analogous with original $\delta b$. To show that the estimate
$\|M^{\frac{1}{2}}\overline{\delta b}\|$ does not strongly influence the final
results, for each case, three random $\overline{\delta b}$ with different
noise levels are generated as follows. We first generated $e=randn(size(b))$,
then put $\overline{\delta b}=g\|b\|e/\|e\|$ where $g$ denotes guessed noise
level. We take $g=0.01,0.02$ and $0.03$ for $2\%$ noise and $g=0.03,0.05$ and
$0.07$ for $5\%$ noise. Table 1 shows exact (i.e. $\|M^{\frac{1}{2}}\delta
b\|$) and estimated values (i.e. $\|M^{\frac{1}{2}}\overline{\delta b}\|$). It
will be seen from Figures below that the error-curves are quite robust versus
the estimates of the noise.

Table 1: Exact and estimated values for $\|M^{\frac{1}{2}}\delta b\|$. |  |  | Estimated value
---|---|---|---
Problem/ noise | # of blocks | Exact value | g=1% | g=2% | g=3%
case-one/ 2% | 8 | 14.04 | 5.07 | 10.14 | 15.22
22 | 19.71 | 8.70 | 17.40 | 26.10
case-two/ 2% | 8 | 8.75 | 4.87 | 9.56 | 14.33
22 | 15.82 | 7.59 | 15.17 | 22.77
|  |  | g=3% | g=5% | g=7%
case-one/ 5% | 8 | 28.36 | 19.29 | 32.14 | 45.01
| 22 | 45.91 | 31.55 | 52.59 | 73.64
case-two/ 5% | 8 | 15.10 | 8.95 | 14.91 | 20.88
22 | 21.57 | 13.48 | 22.47 | 31.46

Following [48], we take $r=1.5$ in the $\Psi_{3}$ rule. As noted above
$r\in(1,2]$ in the $\Gamma-$rule.

Figure 1: The effect of using different values for $"r"$ on $\Gamma$ rule.

In Figure 1 the effect of picking different $r-$values is displayed. As seen,
the choice of proper value of $r$ can have efficacious impact on the rate of
convergence. The value of $r$ should also be chosen in such a way that the
property (44) is satisfied, at least after a few iterations. To insure this,
we take $r=1.5$ and $1.75$ for $2\%$ and $5\%$ noise, respectively. Then (44)
will be satisfied, for $k\geq 2$. We further remark that our computational
experience shows that the unregularized problem $(\alpha=0)$ gives results
that are indistinguishable from those of the regularized problem with
$\alpha=\underline{\sigma}^{2}$. In the sequel we study relative error curves
during $cmax=100$ cycles, but will also extend to $cmax=500$ cycles to study
semiconvergence. Figure 2 illustrates the importance of incorporating
constraints, especially for highly underdetermined systems, during the
iterations.

Figure 2: The effect of using constraints (solid curves) and not using
constraints (dotted curves) in $\Gamma$ and $\Psi_{3}$ rules.

We next discuss Figures 3 and 4 where we display error-curves, and relaxation
parameters for both case-one and case-two using 2$\%$ noise-level. Since the
relaxation parameters showed a similar behavior for the same noise-level we
only display these for case-one. We observe that the relative error is much
smaller for the $\Gamma-$rule than for the best $\Psi-$rule ($\Psi_{3}$). The
reason for this behavior could be that the corresponding relaxation parameters
are bigger using the $\Gamma-$rule than using the $\Psi-$rule.

Figure 3: Case-one with $2\%$ noise. Relaxation parameter behavior (first row)
and relative error history (second row) for $8$ blocks (left column) and $22$
blocks (right column).

Figure 4: Case-two with $2\%$ noise. Relative error history for $8$ blocks
(left column) and $22$ blocks (right column).

In Figures 5 and 6 we show error curves and relaxation parameters (for case-
two) using 5$\%$ noise. It is seen that now there is smaller difference in
relative error between the $\Gamma$-rule and the $\Psi_{3}$-rule. We see from
the figures that the relaxation parameters using the two rules are quite close
in value. We have also listed the minimal error and corresponding cycle number
in Tables 2 and 3.

Figure 5: Case-one with $5\%$ noise. Relative error history for $8$ blocks
(left column) and $22$ blocks (right column).

Figure 6: Case-Two with $5\%$ noise. Relaxation parameter behavior (first row) and relative error history (second row) for $8$ blocks (left column) and $22$ blocks (right column). Table 2: Minimum relative error and corresponding cycle number , when $2\%$ noise applied. Test Problem | # of blocks | Strategy
---|---|---
$\theta_{opt}$ | $\Psi_{3}$ | $\Gamma$
$g=0.01$ | $g=0.02$ | $g=0.03$
Case One | 8 | (66, 0.1531) | (100, 0.2914) | (100, 0.1543) | (100, 0.1622) | (100, 0.1706)
22 | (29, 0.1538) | (100, 0.2295) | (100, 0.1530) | (100, 0.1567) | (100, 0.1613)
Case Two | 8 | (40, 0.1221) | (100, 0.2715) | (100, 0.1265) | (100, 0.1449) | (100, 0.1597)
22 | (15, 0.1219) | (100, 0.2128) | (64, 0.1217) | (100, 0.1237) | (100, 0.1300)
Table 3: Minimum relative error and corresponding cycle number , when $5\%$ noise applied. Test Problem | # of blocks | Strategy
---|---|---
$\theta_{opt}$ | $\Psi_{3}$ | $\Gamma$
$g=0.03$ | $g=0.05$ | $g=0.07$
Case One | 8 | (12, 0.2383) | (100, 0.2914) | (100, 0.2439) | (100, 0.2666) | (100, 0.2866)
22 | (5, 0.2392) | (100, 0.2557) | (97, 0.2398) | (100, 0.2495) | (100, 0.2639)
Case Two | 8 | (15, 0.1947) | (100, 0.2769) | (100, 0.2606) | (100, 0.2356) | (100, 0.2408)
22 | (6, 0.1948) | (100, 0.2559) | (100, 0.1952) | (100, 0.2200) | (100, 0.2313)

We note from Figures 3, 4, 5 and 5 that semi-convergence has hardly begun for
the $\Gamma$-rule. To see more clearly this effect we extend in Figure 7 the
error curves (for two cases) to $cmax=500$. In this Figure we have marked the
point where the error has its minimum, i.e. the cycle number where semi-
convergence starts. We see that the $\theta-$opt rule has the fastest
convergence (note again that we used the exact phantom for training). For this
rule the semi-convergence behavior is quite pronounced, and hence it requires
a reliable stopping criterion (the choice of stopping criterion is an
interesting issue but is not addressed in this paper). In contrast the $\Psi-$
and $\Gamma-$ rules show little effect of semi-convergence, and therefore it
is less critical where the iterations are stopped.

Figure 7: The symbol $\bullet$ indicates the point that semi-convergence
starts.

The behavior of the relative noise error
$\|x^{k}-\bar{x}^{k}\|/\|\bar{x}^{k}\|$, and the relative iteration error
$\|\bar{x}^{k}-x^{*}\|/\|x^{*}\|$ for case-one is shown in Figure 8. Figure 9
shows the phantom and reconstructions using the $\Psi_{3}$ and the
$\Gamma$-rule respectively for case-one. To better judge the quality of the
two reconstructions we display also the corresponding difference images. These
are defined as the difference between the phantom and the respective
reconstruction. More artifacts and noise can be seen in the left image
($\Psi_{3}$) than in the right image ($\Gamma$).

Figure 8: Relative noise error $\|x^{k}-\bar{x}^{k}\|/\|x^{\ast}\|$ (dashed),
and the relative iteration error $\|\bar{x}^{k}-x^{\ast}\|/\|x^{\ast}\|$.

Figure 9: (a) The original image (left), the reconstructed image using
$\Psi_{3}$ strategy (middle) and the reconstructed image using the $\Gamma$
strategy, and (b) the difference images for the $\Psi_{3}$ strategy (left) and
$\Gamma$ strategy (right).

##  6 Conclusion

We define a sequential block-iterative iteration (21) for solving split
feasibility problems in Hilbert space. In this respect it compliments the
simultaneous block iteration given in [15] (defined in the finite dimensional
case). The basic operators involved are weakly regular cutters which includes
e.g. metric projections. A complete convergence analysis is provided. When the
projecting sets are polyhedral it is also shown that the iterates converge
linearly. The projected Block Iterative method (P-BIM) Algorithm 2 is a
special case of our general method. We also consider the noise error of P-BIM,
and derive a new upper bound. This bound generalizes earlier bounds given in
[21, 24, 22, 48]. In particular we extend the noise error analysis of block-
iteration to projected block-iteration, and also to the case when the
relaxation parameters are allowed to depend on the noise. Based on the new
bound a new rule for picking relaxation parameters is derived. We demonstrate
the performance of P-BIM using this and other relaxation parameter rules on
examples taken from tomographic imaging.

## Acknowledgement

We thank two anonymous referees for useful suggestions, and professor Per
Christian Hansen and professor Touraj Nikazad for their comments on an earlier
version.

## References

  * [1] R. Aharoni and Y. Censor, Block-iterative projection methods for parallel computation of solutions to convex feasibility problems, Linear Algebra and Its Applications, 120 (1989), pp. 165–175.
  * [2] Z.-Z. Bai, A. Buccini, K. Hayami, L. Reichel, J.-F. Yin, and N. Zheng, Modulus-based method for constrained Tikhonov regularization, Journal of Computational and Applied Mathematics, 319 (2017), pp. 1–13.
  * [3] H. Bauschke and J. Borwein, On projection algorithms for solving convex feasibility problems, SIAM review, 38 (1996), pp. 367–426.
  * [4] H. H. Bauschke and P. L. Combettes, Convex analysis and monotone operator theory in Hilbert spaces, vol. 408, CMS Books in Mathematics, Springer, Cham, 2017.
  * [5] M. Bertero and P. Boccacci, Introduction to Inverse Problems in Imaging, CRC press, 1998.
  * [6] C. Byrne, Iterative oblique projection onto convex sets and the split feasibility problem, Inverse problems, 18 (2002), pp. 441–453.
  * [7] C. Byrne, Y. Censor, A. Gibali, and S. Reich, The split common null point problem, J. Nonlinear Convex Anal, 13 (2012), pp. 759–775.
  * [8] C. L. Byrne, Block-iterative methods for image reconstruction from projections, IEEE Transactions on Image Processing, 5 (1996), pp. 792–794.
  * [9] A. Cegielski, Iterative Methods for Fixed Point Problems in Hilbert Spaces, vol. 2057 of Lecture Notes in Mathematics, Springer, Heidelberg, 2012.
  * [10] A. Cegielski, General method for solving the split common fixed point problem, Journal of Optimization Theory and Applications, 165 (2015), pp. 385–404.
  * [11] A. Cegielski, S. Reich, and R. Zalas, Regular sequences of quasi-nonexpansive operators and their applications, SIAM Journal on Optimization, 28 (2018), pp. 1508–1532.
  * [12] A. Cegielski, S. Reich, and R. Zalas, Weak, strong and linear convergence of the CQ-method via the regularity of Landweber operators, Optimization, 69 (2020), pp. 605–636.
  * [13] Y. Censor and T. Elfving, A multiprojection algorithm using Bregman projections in a product space, Numerical Algorithms, 8 (1994), pp. 221–239.
  * [14] Y. Censor, T. Elfving, G. T. Herman, and T. Nikazad, On diagonally relaxed orthogonal projection methods, SIAM Journal on Scientific Computing, 30 (2008), pp. 473–504.
  * [15] Y. Censor, T. Elfving, N. Kopf, and T. Bortfeld, The multiple-sets split feasibility problem and its applications for inverse problems, Inverse Problems, 21 (2005), p. 2071.
  * [16] Y. Censor, D. Gordon, and R. Gordon, Component averaging: An efficient iterative parallel algorithm for large and sparse unstructured problems, Parallel computing, 27 (2001), pp. 777–808.
  * [17] Y. Censor and A. Segal, The split common fixed point problem for directed operators, J. Convex Anal, 16 (2009), pp. 587–600.
  * [18] Y. Censor and A. Zaslavski, String-averaging projected subgradient methods for constrained minimization, Optimization Methods and Software, 29 (2014), pp. 658–670.
  * [19] G. Cimmino, Cacolo approssimato per le soluzioni dei systemi di equazioni lineari, La Ricerca Scientifica, II 9 (1938), pp. 326–333.
  * [20] T. Elfving, Row and column based iterations, Applied Analysis and Optimization, 2 (2018), pp. 219–236.
  * [21] T. Elfving, P. C. Hansen, and T. Nikazad, Semiconvergence and relaxation parameters for projected SIRT algorithms, SIAM Journal on Scientific Computing, 34 (2012), pp. A2000–A2017.
  * [22] T. Elfving, P. C. Hansen, and T. Nikazad, Semi-convergence properties of Kaczmarz’s method, Inverse Problems, 30 (2014), p. 055007.
  * [23] T. Elfving and T. Nikazad, Properties of a class of block-iterative methods, Inverse Problems, 25 (2009), p. 115011.
  * [24] T. Elfving, T. Nikazad, and P. C. Hansen, Semi-convergence and relaxation parameters for a class of SIRT algorithms, Electronic Transactions on Numerical Analysis, 37 (2010), pp. 321–336.
  * [25] H. W. Engl, M. Hanke, and A. Neubauer, Regularization of Inverse Problems, vol. 375, Springer Science & Business Media, 1996.
  * [26] P. Gilbert, Iterative methods for the three-dimensional reconstruction of an object from projections, Journal of theoretical biology, 36 (1972), pp. 105–117.
  * [27] K. Goebel and R. Simeon, Uniform convexity, hyperbolic geometry, and nonexpansive mappings, Marcel Dekker, New York and Basel, 1984.
  * [28] R. Gordon, Stop breast cancer now! imagining imaging pathways toward search, destroy, cure, and watchful waiting of premetastasis breast cancer, in Breast Cancer, Springer, 2010, pp. 167–203.
  * [29] R. Gordon, R. Bender, and G. T. Herman, Algebraic reconstruction techniques (art) for three-dimensional electron microscopy and x-ray photography, Journal of theoretical Biology, 29 (1970), pp. 471–481.
  * [30] J. Gregor and T. Benson, Computational analysis and improvement of SIRT, Medical Imaging, IEEE Transactions on, 27 (2008), pp. 918–924.
  * [31] M. Haltmeier, Convergence analysis of a block iterative version of the loping Landweber–Kaczmarz iteration, Nonlinear Analysis: Theory, Methods & Applications, 71 (2009), pp. e2912–e2919.
  * [32] P. C. Hansen and J. S. Jørgensen, AIR Tools II: algebraic iterative reconstruction methods, improved implementation, Numerical Algorithms, 79 (2018), pp. 107–137.
  * [33] G. T. Herman, Fundamentals of Computerized Tomography: Image Reconstruction from Projections, Springer Science & Business Media, 2009\.
  * [34] H. M. Hudson and R. S. Larkin, Accelerated image reconstruction using ordered subsets of projection data, IEEE transactions on medical imaging, 13 (1994), pp. 601–609.
  * [35] M. Jiang and G. Wang, Convergence studies on iterative algorithms for image reconstruction, Medical Imaging, IEEE Transactions on, 22 (2003), pp. 569–579.
  * [36] S. Kaczmarz, Angenäherte Auflösung von Systemen linearer Gleichungen, Bulletin International de l’Academie Polonaise des Sciences et des Lettres, 35 (1937), pp. 355–357.
  * [37] A. C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging, SIAM, 2001.
  * [38] S. Kindermann and A. Leitao, Convergence rates for Kaczmarz-type regularization methods., Inverse Problems & Imaging, 8 (2014).
  * [39] V. I. Kolobov, S. Reich, and R. Zalas, Finitely convergent deterministic and stochastic iterative methods for solving convex feasibility problems, Mathematical Programming, (2021), pp. 1–21, https://doi.org/10.1007/s10107–021–01628–z.
  * [40] A. Moudafi, The split common fixed-point problem for demi-contractive mappings, Inverse problems, 26 (2010), p. 05507.
  * [41] A. Moudafi, A note on the split common fixed-point problem for quasi-nonexpansive operators, Nonlinear Analysis: Theory, Methods & Applications, 74 (2011), pp. 4083–4087.
  * [42] A. Moudafi and A. Gibali, $\ell_{1}$-$\ell_{2}$ regularization of split feasibility problems, Numerical Algorithms, 78 (2018), pp. 739–757.
  * [43] F. Natterer, The Mathematics of Computerized Tomography, vol. 32, SIAM, 1986.
  * [44] I. Necoara, Faster randomized block Kaczmarz algorithms, SIAM Journal on Matrix Analysis and Applications, 40 (2019), pp. 1425–1452.
  * [45] A. Nedic and D. P. Bertsekas, Incremental subgradient methods for nondifferentiable optimization, SIAM Journal on Optimization, 12 (2001), pp. 109–138.
  * [46] D. Needell and J. A. Tropp, Paved with good intentions: analysis of a randomized block Kaczmarz method, Linear Algebra and its Applications, 441 (2014), pp. 199–221.
  * [47] T. Nikazad, M. Abbasi, and T. Elfving, Error minimizing relaxation strategies in Landweber and Kaczmarz type iterations, Journal of Inverse and Ill-posed Problems, 25 (2017), pp. 35–56.
  * [48] T. Nikazad and M. Karimpour, Controlling noise error in block iterative methods, Numerical Algorithms, 73 (2016), pp. 907–925.
  * [49] T. Nikazad and M. Mirzapour, Projected non-stationary simultaneous iterative methods, Int. J. Nonlinear Anal. Appl, 7 (2016), pp. 243–251.
  * [50] Z. Opial, Weak convergence of the sequence of successive approximations for nonexpansive mappings, Bull. Amer. Math. Soc, 73 (1967), pp. 591–597.
  * [51] M. Piana and M. Bertero, Projected Landweber method and preconditioning, Inverse Problems, 13 (1997), p. 441.
  * [52] G. Qu, C. Wang, and M. Jiang, Necessary and sufficient convergence conditions for algebraic image reconstruction algorithms, IEEE Transactions on Image Processing, 18 (2008), pp. 435–440.
  * [53] S. Reich, M. T. Truong, and T. N. H. Mai, The split feasibility problem with multiple output sets in Hilbert spaces, Optimization Letters, (2020), pp. 1–19.
  * [54] S. Reich and T. M. Tuyen, Projection algorithms for solving the split feasibility problem with multiple output sets, Journal of Optimization Theory and Applications, 190 (2021), pp. 861–878.
  * [55] P. Richtárik and M. Takác, Stochastic reformulations of linear systems: algorithms and convergence theory, SIAM Journal on Matrix Analysis and Applications, 41 (2020), pp. 487–524.
  * [56] Y. Saad, Iterative Methods for Sparse Linear Systems, SIAM, 2003.
  * [57] H. H. B. Sørensen and P. C. Hansen, Multicore performance of block algebraic iterative reconstruction methods, SIAM Journal on Scientific Computing, 36 (2014), pp. C524–C546.
  * [58] R. Tibshirani, Regression shrinkage and selection via the lasso, Journal of the Royal Statistical Society: Series B (Methodological), 58 (1996), pp. 267–288.
  * [59] F. Wang and H.-K. Xu, Cyclic algorithms for split feasibility problems in Hilbert spaces, Nonlinear Analysis: Theory, Methods & Applications, 74 (2011), pp. 4105–4111.
  * [60] H.-K. Xu, Iterative methods for the split feasibility problem in infinite-dimensional Hilbert spaces, Inverse Problems, 26 (2010), p. 105018.
  * [61] M. Yukawa, K. Slavakis, and I. Yamada, Multi-domain adaptive filtering by feasibility splitting, in 2010 IEEE International Conference on Acoustics, Speech and Signal Processing, IEEE, 2010, pp. 3814–3817.
  * [62] E. Zeidler, Nonlinear Functional Analysis and its Applications, IIB: Nonlinear Monotone Operators, Springer, New York, 1990.