 # Chebyshev Inertial Landweber Algorithm for Linear Inverse Problems

The Landweber algorithm defined on complex/real Hilbert spaces is a gradient descent algorithm for linear inverse problems. Our contribution is to present a novel method for accelerating convergence of the Landweber algorithm. In this paper, we first extend the theory of the Chebyshev inertial iteration to the Landweber algorithm on Hilbert spaces. An upper bound on the convergence rate clarifies the speed of global convergence of the proposed method. The Chebyshev inertial Landweber algorithm can be applied to wide class of signal recovery problems on a Hilbert space including deconvolution for continuous signals. The theoretical discussion developed in this paper naturally leads to a novel practical signal recovery algorithm. As a demonstration, a MIMO detection algorithm based on the projected Landweber algorithm is derived. The proposed MIMO detection algorithm achieves much smaller symbol error rate compared with the MMSE detector.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## I Introduction

A number of signal detection problems in wireless communications and signal processing can be classified into linear inverse problems. In a linear inverse problem, a source signal

( or ) is inferred from the noisy linear observation where and the additive noise .

One of the simplest approaches for the above task is to rely on the least square principle, i.e., one can minimize

which corresponds to the maximum likelihood estimation rule for estimating the source signal

under the Gaussian noise assumption. The Landweber algorithm   is defined by the fixed-point iteration

 x(k+1)=x(k)−ωH†(Hx(k)−y), k=0,1,2,… (1)

where if , otherwise . The notation indicates the Hermitian transpose of and the parameter is a real constant. Note that the Landweber algorithm can be defined not only on a finite dimensional Euclidean space but also on an infinite dimensional complex/real Hilbert space. The Landweber algorithm defined on a Hilbert space is especially important for linear inverse problems including convolutions of continuous signals.

The Landweber algorithm can be regarded as a gradient descent method for minimizing the objective function because is the gradient of the objective function. The Landweber algorithm has been widely employed as a signal reconstruction algorithm for image deconvolution 

, inverse problems regarding diffusion partial differential equations

, MIMO detectors , and so on. An advantage of the Landweber algorithm is that we can easily include a proximal or projection operation that utilizes the prior knowledge on the source signal after the gradient descent step (1), which is often called the projected Landweber algorithm .

One evident drawback of the Landweber algorithm is that the convergence speed is often too slow and we need to exploit an appropriate acceleration method. Recently, Takabe and Wadayama  found that a step-size sequence determined from Chebyshev polynomials can accelerate the convergence of gradient descent algorithms. Wadayama and Takabe  generalized the central idea of  for general fixed-point iterations. The method is called the Chebyshev inertial iteration. It would be natural to apply the Chebyshev inertial iteration for improving the Landweber algorithm in terms of the convergence speed.

The successive over relaxation (SOR) is a common method for accelerating a fixed-point iteration for solving linear equation such as Jacobi method. The SOR iteration corresponding to the above equation is given by where . A fixed SOR factors is commonly used in practice. The Chebyshev inertial iteration  is a natural generalization of the SOR method, which is a method employing defined based on the Chebyshev polynomials. The fundamental properties including the convergence rate are analyzed in [8, 9].

The goal of this paper is twofold. The first goal is to extend the theory of the Chebyshev inertial iteration to the Landweber algorithm on Hilbert spaces. The arguments in [8, 9] are restricted to the case where the underlying space is a finite dimensional Euclidean space. By extending the argument to a Hilbert space, the essential idea of the Chebyshev inertial iteration becomes applicable to iterative algorithms for infinite dimensional linear inverse problems, and we can obtain accelerated convergence of these algorithms. The Chebyshev inertial Landweber algorithm presented in this paper can be applied to wide class of signal reconstruction algorithms on complex/real Hilbert spaces such as deconvolution for continuous signals.

The second goal is to present a novel MIMO detection algorithm based on the projected Landweber algorithm with the Chebyshev inertial iteration for demonstrating that the theoretical discussion developed in this paper naturally leads to a novel practical signal detection algorithm.

## Ii Preliminaries

In this section, several basic facts on functional analysis required for the subsequent argument are briefly reviewed. Precise definitions regarding functional analysis presented in this section can be found in .

### Ii-a Hilbert space

Let be a real number satisfying . The set of infinite sequences in satisfying is denoted by where or

. If the length of sequences are finite, i.e., a complex or real vector space, we can define a finite dimensional norm space

. The set of measurable functions defined on the closed interval , is denoted by .

Let be a complex or real vector space. For any , an inner product is assumed to be given. The norm (inner product norm) associated with the inner product is defined by If the norm space is complete, the space is said to be Hilbert space. The norm spaces and are Hilbert spaces. In the case of , the inner products are defined by On the other hand, the inner product for is defined by

Let be a bounded linear operator on . The operator norm of is defined by where is an element of . For any and an operator , we have the sub-multiplicative inequality

 ∥Tx∥≤∥T∥∥x∥. (2)

### Ii-B Compact operator and spectral mapping theorem

Let be a linear operator on . A complex number

is an eigenvalue of

if and only if a nonzero vector in satisfies where

is said to be an eigenvector corresponding to

. In the following, the set of all eigenvalues of is denoted by . The spectral radius of is given as

 ρ(T):=max{|λ|:λ∈σ(T)}. (3)

Suppose that satisfies For any , let which is a linear operator on and the operator is a compact operator . If the kernel function satisfies then the operator is a compact self-adjoint operator. In the case of the finite dimensional space , a linear operator defined by a Hermitian matrix corresponds to a compact self-adjoint operator.

The next theorem plays an important role in our analysis described below.

###### Theorem 1 (Spectral mapping theorem )

Let be a compact operator defined on . If a complex valued function is analytic around , then the set of eigenvalues of is given by

Let be a compact self-adjoint operator on . It is known that the set of eigenvalues of a compact self-adjoint operator is a countable set of real numbers . Let . A polynomial defined on is analytic over whole . As a consequence of the spectral mapping theorem, the set of eigenvalues of is given as For a compact self-adjoint operator , we also have

 ∥f(T)∥=ρ(f(T)). (4)

### Ii-C Landweber iteration

Assume that is an element in a Hilbert space. The aim of the Landweber iteration is to recover the original signal from a measurement or noisy measurement where is a linear operator. The Landweber iteration is a gradient descent algorithm in a Hilbert space to minimize the error functional The Landweber iteration is defined by the fixed-point iteration

 x(k+1)=x(k)−ωT∗(Tx(k)−y). (5)

## Iii Chebyshev Inertial Iteration

### Iii-a Fixed-point iteration

Suppose that we have a fixed-point iteration defined on a Hilbert space :

 x(k+1)=Ax(k)+b,k=0,1,2,…, (6)

where . The operator is a compact self-adjoint operator on .

Assume also that is a contraction mapping which has a fixed point satisfying where the existence of the fixed point is guaranteed by Banach fixed-point theorem . The inertial iteration  corresponding to the original iteration (6) is defined by

 x(k+1)=x(k)+ωk(Ax(k)+b−x(k)), (7)

where are called the inertial factors. It should be remarked that the fixed point of the inertial iteration exactly coincides with the fixed point of the original iteration (6).

### Iii-B Analysis on spectral radius

From the inertial iteration (7), we have the equivalent update equation

 x(k+1)=(I−ωkB)x(k)+ωkb, (8)

where . From the assumption that is a compact self-adjoint operator, we can show is also compact self-adjoint. Since the fixed point satisfies

 x∗=(I−ωkB)x∗+ωkb, (9)

for any nonnegative and subtracting (9) from (8), we immediately have a recursive formula on the residual errors:

 x(k+1)−x∗=(I−ωkB)(x(k)−x∗). (10)

In the following discussion, we will assume the following inertial factors satisfying the periodical condition:

 ωℓT+j=ωj,ℓ=0,1,2,…,j=0,1,…,T−1, (11)

where a positive integer represents the period.

Recursively applying (10) multiple times, we can get

 x(T)−x∗=(T−1∏k=0(I−ωkB))(x(0)−x∗). (12)

The norm of the left-hand side can be upper bounded by

 ∥x(T)−x∗∥ = ∥∥ ∥∥(T−1∏k=0(I−ωkB))(x(0)−x∗)∥∥ ∥∥ (13) ≤ ∥∥ ∥∥T−1∏k=0(I−ωkB)∥∥ ∥∥∥x(0)−x∗∥ (14) = ρ(T−1∏k=0(I−ωkB))∥x(0)−x∗∥. (15)

The above inequality is based on the sub-multiplicative inequality (2) and the last equality is due to (4).

### Iii-C Chebyshev inertial factors

As we discussed, the operator is a compact self-adjoint operator and it has countable real eigenvalues, which are denoted by . Hereafter, the minimum and the maximum eigenvalues of are denoted by and , respectively.

In the following discussion, we will use the inertial factors defined by

 ω∗k:=[λ++λ−cos(2k+12Tπ)]−1, k=0,1,…,T−1, (16)

where These inertail factors are called the Chebyshev inertial factors. A Chebyshev inertial factor is the inverse of a root of a Chebyshev polynomial [8, 9].

Let . If the Chebyshev inertial factors are employed in the inertial iteration, we can use Lemma 2 proved in  to bound as

 |βT(λ)|≤sech(Tcosh−1(λ+λ−)). (17)

for . Combining the spectral mapping theorem and the above inequality, we immediately obtain

 ρ(βT(B))=ρ(T−1∏k=0(I−ωkB))≤% sech(Tcosh−1(λ+λ−)). (18)

The above argument can be summarized in the following theorem.

###### Theorem 2

Let be a positive integer. For , the Landweber iteration with the Chebyshev inertial factors satisfies the residual error bound:

 ∥x(ℓT)−x∗∥≤U(T)ℓ∥x(0)−x∗∥, (19)

where .

This theorem implies that the error norm is tightly upper bounded if the iteration index is a multiple of .

Combining the Landweber iteration with the Chebyshev inertial iteration, we have the Chebyshev inertial Landweber algorithm:

 x(k+1) = x(k)+ω∗k(x(k)−ωT∗(Tx(k)−y)−x(k)) (20) = x(k)−ω∗kωT∗(Tx(k)−y).

The iteration is almost the same as the original Landweber iteration except for the use of the Chebyshev inertial factors . It is common to set the inverse of the maximum eigenvalue of to the fixed factor .

## Iv Experiments

### Iv-a Deconvolution by Landweber iteration

In this subsection, deconvolution by the Landweber iteration over the Hilbert space over is studied. Suppose that are measurable functions on . The convolution of and are defined as

 y(x):=∫∞−∞f(u)g(x−u)du, (21)

which is a compact operator on . We thus can write where is a compact operator defined by (21) and . In a context of an inverse problem regarding a linear system involving a convlution, the function can be considered as a source signal and represents a point spread function (PSF) or impulse response. In the context of a partial differential equation, represents the Green’s function or the integral kernel corresponding to a given partial differential equation.

We further assume that is an even function satisfying . This implies that the operator is a compact self-adjoint operator. In the following, we try to recover the source signal from the blurred signal by using the Landweber iteration on given by where the initial condition is Note that holds due to the assumption on .

In the following experiments shown below, the closed range from to are discretized into 16384 bins, and convolution integration (21

) is approximated by cyclic convolution with 16384-points FFT and frequency domain products. Figure

1 presents the source signal , the convolutional kernel , and the convolved signal . Fig. 1: Plots of the source signal f, convolutional kernel g, and the convolved signal y. The definitions of f and g are given by f(x):=12exp(−x2)+exp(−(x−2)2)−exp(−(x−3)2)+12exp(−(x−4)2),g(x):=exp(−x2). Fig. 2: Deconvolution process: (upper) original Landweber iteration, (lower) Chebyshev inertial Landweber iteration. The coefficient ω is set to 0.3 in the both cases. For Chebyshev inertial Landweber algorithm, lmin=0.1 and lmax=0.9 are used. The period T is set to 8. The index k represents the number of iteration.

Figure 2 presents a deconvolution process by the original Landweber algorithm (upper) and the Chebyshev inertial Landweber algorithm (lower). The tentative results for every 30 iterations are shown. In both cases, we can observe that gradually approaches to the source signal . This is because is a contractive mapping under the setting . Comparing two figures, we can see that the Chebyshev inertial Landweber algorithm shows much faster convergence to the source signal . This observation is also supported by the error curves depicted in Fig.3. The Chebyshev inertial Landweber algorithm shows faster convergence compared with the original Landweber iteration. For example, the error is achieved with with the original Landweber iteration but the Chebyshev iteration requires only iterations. Fig. 3: Error norms ||s(k)−f|| as functions of number of iterations in deconvolution processes.

### Iv-B Finite-dimensional least square problem

In this subsection, we will discuss a least square problem closely related to MIMO detection problems. We here assume a finite dimensional Hilbert space on .

Let

be a complex matrix whose elements follow the normal complex Gaussian distribution

. Our task is to recover a source signal from a noisy linear measurement where is a complex Gaussian noise vector. The problem can be seen as a MIMO detection problem where the numbers of transmit and receive antennas are . The least square problem is equivalent to the maximum likelihood estimation rule for the above setting.

In this case, the Landweber iteration is exactly same as a gradient descent process

 s(k+1)=s(k)−ωHH(Hs(k)−y) (22)

for minimizing the objective function. It is known that the asymptotic convergence rate is optimal if holds. In this case, the matrix becomes a Hermitian matrix and thus is a compact self-adjoint operator on . This means that we can apply the Chebyshev inertial iteration to the Landweber iteration.

The details of the experiment are summarized as follows. The dimension is set to . Each element of is chosen from 8-PSK constellation uniformly at random. The optimal factor is employed in the Landweber iteration. Each element of the noise vector follows where . The minimum and maximum eigenvalues of are used as and for determining the Chebyshev inertial factors.

Figure 4 summarizes the comparison on the squared error norms between the original Landweber iteration and the accelerated iterations. From Fig. 4 (left), we can immediately recognize the zigzag-shaped error curves of the Chebyshev inertial Landweber algorithm. This is because the error norm is exponentially upper bounded only when the iteration index is multiple of as shown in Theorem 2. It mean that the lower envelope of the zigzag curves can be seen as the guaranteed performance of the Chebyshev inertial Landweber algorithm. Figure 4 (right) indicates the squared errors up to . The proposed schemes (Chebyshev ) shows much faster convergence compared with the original Landweber iteration. Fig. 4: Squared error norms ∥s(k)−x∥2 as a function of number of iterations: (left) number of iteration is up to 50, (right) number of iteration is up to 30000. The squared error norms are averaged for 100-trials.

The speed of convergence of the Landweber iteration is dominated by the spectral radius On the other hand, indicates an upper bound of the normalized asymptotic convergence rate depending on the period . Figure 5 (left) presents both and as functions of . As becomes larger, the value of becomes strictly smaller than . This means that the Chebyshev inertial iteration certainly accelerates the asymptotic convergence rate. Figure 5 (right) shows approximate error norms derived from and which are given by . Although these values does not include the initial errors , these curves qualitatively explains the behavior of the actual error curves in Fig. 4 (right). These results support the usefulness of the theoretical analysis discussed in the previous section. Fig. 5: (left) Normalized convergence rate ρ(A) and U(T), (right) approximate error norms derived from ρ(A) and U(T): the approximated error norms are ρ(A)k,U(2)k,U(8)k, respectively.

### Iv-C Projected Landweber-based MIMO detection

It is natural to consider the projected Landweber algorithm  for making use of the prior information of the source signal to improve the reconstruction performance. We here examine the detection performance of a simple projected Landweber-based MIMO detection algorithm.

The soft projection operator used here is defined by

 η(r):=∑p∈Spexp(−|r−p|2/α2)∑p∈Sexp(−|r−p|2/α2), (23)

where is a signal constellation . The iteration of the projected Landweber algorithm is fairly simple; the soft projection operator is just element-wisely applied to the output in (22) and then, the output is passed to the inertial iteration process. The details of the experiments are as follows. The channel model is exactly the same as that used in the previous subsection ( MIMO channel with 8-PSK input). The minimum and maximum eigenvalues and are determined according to the Marchenko-Pastur law. Figure 6 presents the symbol error rate (SER) of the proposed scheme. The projected Landweber detectors achieves one or two order of magnitude smaller SER than those of the MMSE detector. The accelerated schemes (Chebyshev ) provide steeper error curves than that of the projected Landweber detector without Chebyshev inertial iteration (denoted by “Landweber”). This observation supports the potential of the Chebyshev inertial iteration. Fig. 6: Symbol error rate of the projected Landweber algorithm with the Chebyshev inertial iteration. The channel model is 32×32 MIMO channel with 8-PSK input. The number of iterations for all the algorithms (except for MMSE) is set to 100. As baselines, the performance of the MMSE detector defined by ^x:=(HHH+σ2I)−1HHy is included. The SNR snr (in dB) and σ is related as σ=√10−snr/10. In the projected Landweber algorithms, α2=0.5 is used when iteration index k is less than 20; for k≥20, α2=0.25 is used.

## Acknowledgement

This work was supported by JSPS Grant-in-Aid for Scientific Research (B) Grant Number 19H02138.