Claudio Estatico
email: estatico@dima.unige.it

Serge Gratton
email: serge.gratton@enseeiht.fr

Flavia Lenti
email: flavialenti@gmail.com

David Titley-Peloquin
email: david.titley-peloquin@mcgill.ca

A conjugate gradient like method for p-norm minimization in functional spaces

Keywords: Iterative regularization, Banach spaces, Conjugate gradient method

We develop an iterative algorithm to recover the minimum p-norm solution of the functional linear equation Ax = b, where A p : X -→ Y is a continuo r us linear operator between the two Banach spaces X = L , 1 < p < 2, and Y = L , r > 1, with x ∈ X and b ∈ Y. The algorithm is conceived within the same frame-work of the Landweber method for functional linear equations in Banach spaces proposed by Schöpfer, Louis and Schuster (Inverse Probl., 22:311-329, 2006). Indeed, the algorithm is based on using, at the n-th iteration, a linear combination of the steepest current "descent functional" A * J (b -Ax n ) and the previous descent functional, where J denotes a duality map of the Banach space Y. In this regard, the algorithm can be viewed as a generalization of the classical conjugate gradient method on the normal equations (CGNR) in Hilbert spaces. We demonstrate that the proposed iterative algorithm converges strongly to the minimum p-norm solution of the functional linear equation Ax = b and that it is also a regularization method, by applying the discrepancy principle as stopping rule.

According to the geometrical properties of L p spaces, numerical experiments show that the method is fast, robust in terms of both restoration accuracy and stability, promotes sparsity and reduces the over-smoothness in reconstructing edges and abrupt intensity changes.

Introduction

We discuss the problem of iteratively recovering a solution of the functional equation Ax = b ,

where A : X -→ Y is a linear and continuous operator between two functional spaces X and Y, with x ∈ X and b ∈ Y.

The problem we consider is ill-posed, since its solution is not unique or does not exist for arbitrary data or does not depend continuously on the data [START_REF] Bertero | Introduction to inverse problems in imaging[END_REF]. There is now an exhaustive and comprehensive literature about linear and nonlinear inverse ill-posed problems in Hilbert spaces. Over the last decade, there has been a growing interest in studying inverse ill-posed problems in Banach spaces, because many applicative problems are therein better modeled ( [START_REF] Bonesky | Minimization of Tikhonov functional in Banach spaces[END_REF], [START_REF] Daubachies | An iterative thresholding algorithm for linear inverse problems with sparsity constraint[END_REF], [START_REF] Resmerita | Regularization of ill-posed problems in Banach spaces: convergence rates[END_REF], [START_REF] Schöpfer | Nonlinear iterative methods for linear ill-posed problems in Banach spaces[END_REF], [START_REF] Schöpfer | Metric and Bregman projections onto affine subspaces and their computation via sequential subspace optimization methods[END_REF], [START_REF] Schöpfer | Fast regularizing sequential subspace optimization in Banach spaces[END_REF], [START_REF] Schöpfer | An iterative regularization method for the solution of the split feasibility problem in Banach spaces[END_REF], [START_REF] Schuster | Regularization methods in Banach spaces[END_REF], to name just a few references).

In [START_REF] Schöpfer | Nonlinear iterative methods for linear ill-posed problems in Banach spaces[END_REF], Schöpfer, Louis and Schuster computed a regularized solution in Banach spaces by means of the following generalized Landweber, i.e. gradient-type, iterative scheme

x n+1 = j X * " j X (x n ) -α n A * j Y (Ax n -b) " , (2)
where j X , j X * and j Y are duality mappings of the corresponding spaces X , X * , Y [START_REF] Schuster | Regularization methods in Banach spaces[END_REF], and α n > 0 is a appropriately chosen step size. The discrepancy principle is employed to obtain a suitable stopping index in the case of noisy data. Like in Hilbert spaces, the method turned out to have good regularizing properties but its convergence speed is generally very slow. The conjugate gradient (CG) algorithm for linear systems is known to enjoy better convergence properties than methods based solely on gradient descent, but its original formulation by Hestenes and Stiefel applies only in Hilbert spaces. In this respect, the algorithm we develop in this paper can be considered as an extension to Banach space settings of the conventional CG method on the normal equations for linear systems in Hilbert spaces (also known as CGNR method). Indeed, on the basis of the theoretical framework of [START_REF] Schöpfer | Nonlinear iterative methods for linear ill-posed problems in Banach spaces[END_REF], we suggest using a different "descent functional" than the steepest one A * j Y (b -Ax n ) at each iterative step, which is defined by the same kind of linear combination as in the CGNR method. More specifically, in [START_REF] Schöpfer | Nonlinear iterative methods for linear ill-posed problems in Banach spaces[END_REF] the spaces X is assumed to be smooth and uniformly convex and Y can even be an arbitrary Banach space. Although our convergence proof closely follows the arguments of [START_REF] Schöpfer | Nonlinear iterative methods for linear ill-posed problems in Banach spaces[END_REF], in this paper we restrict the analysis to the case of X = L p with 1 < p < 2 and Y = L r with r > 1. On the other hand, we argue that a generalization to other Banach spaces could be obtained without strong changes.

To give a numerical validation, we apply our algorithm to a standard image restoration problem. The exact solution presents discontinuities characterized by different dimensions and different intensities. Reconstructions obtained by using the proposed approach are compared with the ones obtained by CG in Hilbert spaces, and by both steepest descent and Landweber methods in Hilbert and Banach spaces. Since the proposed method is based on the minimization of p-norm (with 1 < p < 2), it is expected to overcome the typical over-smoothing drawback effects of regularization in Hilbert spaces and to enhance the sparsity of the reconstructed solution. Our first numerical evidences confirm such a positive behavior in both convergence speed and accuracy of the restoration.

Due to the intrinsic non-linearity of the duality maps, we mention that, differing from original CG method in Hilbert spaces, the short term recursion formula for computing the descent functionals (which have the same role of the descent directions in Hilbert spaces) does not guarantee now their full mutual conjugacy. Thus, the proposed algorithm does not converge within n iterations for n×n linear systems in Banach spaces. A concomitant proposal by Herzog and Wollner [START_REF] Herzog | A conjugate direction method for linear systems in Banach spaces[END_REF] for the linear system (1) in a reflexive Banach space with A : X -→ X * self-adjoint and positive, preserves the conjugacy of the descent directions, so that convergence in a finite number of steps holds as in the finite dimensional Hilbertian case. Another positive fact is that, differing form the proposed algorithm, the method in [START_REF] Herzog | A conjugate direction method for linear systems in Banach spaces[END_REF] does not require to solve any one-dimensional minimization problem for the computation of the optimal step size at each iteration, since a closed formula holds as in the Hilbertian case. On the other hand, to keep the mutual conjugacy of any new descent direction, such an algorithm requires the storage of all the previous descent directions, which are all used for the implementation of a modified Arnoldi conjugation procedure as the iterations go on. Although finite convergence is not assured, the CG-like method proposed in the present paper does not require any other storage than the last descent (functional) direction and can be used in a more general setting, since the operator A is here not required to be self-adjoint. Moreover, differing form [START_REF] Herzog | A conjugate direction method for linear systems in Banach spaces[END_REF], the iterations are here explicitly computed in the dual space X * , as generally done in [START_REF] Schöpfer | Nonlinear iterative methods for linear ill-posed problems in Banach spaces[END_REF], [START_REF] Schöpfer | Metric and Bregman projections onto affine subspaces and their computation via sequential subspace optimization methods[END_REF] and [START_REF] Bonesky | Minimization of Tikhonov functional in Banach spaces[END_REF].

The remainder of this paper is organized as follows. In Section 2, we give the necessary theoretical tools used in Section 3 and Section 4 to prove the strong convergence for noise-free and noisy data b, respectively. In Section 5, a practical implementation issue related to the step size computation is described and in Section 6 the numerical experiments are presented and discussed. Some conclusions are drawn in Section 7.

Preliminaries

The aim of this section is to briefly recall basic tools and classical notations usually used in the Banach space setting for the regularization of ill-posed problems. For details and proofs, we suggest the following monographs: [START_REF] Cioranescu | Geometry of Banach spaces, duality mappings and nonlinear problems[END_REF], [START_REF] Lindestrauss | Classical Banach spaces[END_REF], [START_REF] Schuster | Regularization methods in Banach spaces[END_REF]. Following the same notation as [START_REF] Schöpfer | Nonlinear iterative methods for linear ill-posed problems in Banach spaces[END_REF], throughout the paper, both X and Y are two real Banach spaces with dual spaces X * and Y * equipped with the corresponding operator norm . For x * ∈ X * and x ∈ X we denote by x * , x X * ×X and x, x * X ×X * the duality pairing defined as x * , x X * ×X = x, x * X ×X * = x * (x). In general, we omit subscripts indicating the space when this piece of information is implicitly clear. In particular, this will be done for any norm of vectors or oper-ators throughout the paper. The operator A : X -→ Y is continuous and linear, and A * : Y * -→ X * denotes its adjoint operator of A, that is, the operator such that Ax, y = x, A * y , ∀x ∈ X and y ∈ Y * . We have A = A * in each corresponding operator norm. For two real numbers a and b, we write a ∨ b = max {a, b} , a ∧ b = min {a, b} . Throughout the paper, for p, r ∈ (1, +∞), we usually denote by q = p * and s = r * their Hölder conjugates.

Duality mapping

The key point in the generalization of the regularization method in Banach spaces is the duality mapping [START_REF] Asplund | Positivity of duality mappings[END_REF], [START_REF] Kien | The normalized duality mapping and two related characteristic properties of a uniformly convex Banach space[END_REF], [START_REF] Milicic | On moduli of expansion of the duality mapping of smooth Banach spaces[END_REF]. A duality map is an appropriate function which associates an element of a Banach space X with an element of its dual X * , and it is useful when the Banach space X is not isometrically isomorphous to its dual X * . Formally, we have the following definition. Definition 1 (Duality mapping) The (set-valued) mapping J X p : X -→ 2 X * with p ≥ 1 defined by

J X p (x) = n x * ∈ X * : x * , x = x x * , x * = x p-1 o , is called duality map of X with gauge function t -→ t p-1 .
In general, by j X p (x) ∈ X we will denote a single-valued selection of the subset J X p (x) ⊂ X . If X = l r or X = L r , with 1 < r < +∞, the duality map is a single-valued function which will be denoted as j r p (•). If r = p, the apex will be usually omitted.

Example 1 Let us consider x ∈ l r (R n ). For every p ∈ (1, +∞), the duality map j r p is given by j

r p (x) = x p-r |x| r-1 sign(x) , (3)
where sign(•) denotes the sign function and the product has to be considered as component-wise. In particular, for the Hilbert space l 2 , according to the Riesz representation Theorem, j 2 is the identity operator, i.e. j 2 (x) = x.

The meaning of any duality mapping is naturally related to the sub-gradient of the Banach norm. We first recall the following basic definition. Definition 2 (Subgradient of convex functional)

Let f : X -→ R ∪ {+∞} . Then, x * ∈ X * is a subgradient of f at x ∈ X if f (y) ≥ f (x) + x * , y -x , ∀y ∈ X .
The set ∂f (x) ⊂ X * of all subgradients of f at x is called subdifferential of f at x. The following important theorem give us a simple and heuristic way to understand the meaning of any duality map.

Theorem 1 (Asplund Theorem) [START_REF] Asplund | Positivity of duality mappings[END_REF] Let X be a Banach space and p > 1. Then

J X p = ∂ " 1 p • p « . (4)
The Asplund Theorem is a key step for an intuitive "geometrical" interpretation of any Banach space and its duality maps. In our context, it will be useful in the following sections for the computation of the subdifferential of the residual functional 1 p Axb p by the chaining rule.

Geometry of Banach spaces and Bregman distance

Regarding the regularization theory in Banach spaces, the geometrical properties of the spaces, like convexity and smoothness, play a crucial role [START_REF] Cioranescu | Geometry of Banach spaces, duality mappings and nonlinear problems[END_REF]. These properties, which can be viewed as an extension to Banach space of heuristic basic properties of Hilbert (or simply Euclidean) spaces, are strictly related to the duality maps and will be useful to prove the convergence of the proposed algorithm.

The "degree" of convexity or smoothness of a Banach space is given by its modulus of convexity and its modulus of smoothness, and we refer to Section 2.1 of [START_REF] Schöpfer | Nonlinear iterative methods for linear ill-posed problems in Banach spaces[END_REF] for their formal definitions and for the subsequent characterization of p-convex, p-smooth, uniformly convex and uniformly smooth Banach spaces.

In [START_REF] Burger | Convergence rates of convex variational regularization[END_REF], it has been argued that the Bregman distance is the correct measure for measuring the "quality" of the regularized solution in the case of convex regularization in Banach spaces. Moreover, due to the geometrical properties of Banach spaces, it is often more appropriate to exploit the Bregman distance between two vectors x and y instead of more conventional norm-distances like xy p or j X p (x)j X p (y) p to prove convergence of the algorithms [START_REF] Butnariu | Bregman distances, totally convex functions and a method for solving operator equations in Banach spaces[END_REF], [START_REF] Grasmair | Generalized Bregman distances and convergences rate for non-convex regularization methods[END_REF], [START_REF] Schöpfer | Nonlinear iterative methods for linear ill-posed problems in Banach spaces[END_REF]. Generally, the Bregman distance associated to a convex functional is defined as the difference between the functional and its linear approximation around x as follows [START_REF] Bregman | The relaxation method of finding the common points of convex sets and its application to the solution of problems in convex programming[END_REF].

Definition 3 (Bregman distance)

Let f : X -→ R be a convex and contiously-differentiable functional on a Banach space X . Then the Bregman distance ∆(x, y) of f between x ∈ X and y ∈ X is defined as

∆(x, y) = f (x) -f (y) -f ′ (y)(x -y),
where f ′ (y) is the first order Taylor expansion of f at y.

The role of a Bregman distance is similar to the role of any metric. However, any Bregman distance in general does not satisfy the triangle inequality nor symmetry. As shown by the following special case first example, the Bregman distance is a generalization of the square norm distance in Hilbertian contexts, when the basic square norm distance functional is considered.

Example 2 If X is an Hilbert space, then ∆ 2 (x, y)

= 1 2 x-y 2 , for f 2 (x) = 1 2 x 2 .
Example 3 In a general Banach space X , by considering the convex functional f p (x) = 1 p x p , thanks to the Asplund Theorem (4) the associated Bregman distance ∆ p is

∆ p (x, y) = 1 p x p - 1 p y p -j p (y), x -y , (5)
for any x, y ∈ X , where j p is a single-valued selection of J X p .

In general ∆ p (x, y) ≥ 0, and ∆ p (x, y) = 0 if and only if j p (y) ∈ J X p (x) [START_REF] Schuster | Regularization methods in Banach spaces[END_REF]. From that, for all x, z ∈ X we have

1 p z p -1 p x p -j p (x), z -x ≥ 0 , which, by setting y = -(z -x), yields 1 p x -y p - 1 p x p + j p (x), y ≥ 0. ( 6
)
In general, the Bregman distance in the primal space X and the Bregman distance in its dual space X * are strictly linked, since

∆ p (x, y) = ∆ q `jp (x), j p (y)
´.

This duality, together with basic relationships about strong convergence in norm and convergence in Bregman distance (see [START_REF] Schuster | Regularization methods in Banach spaces[END_REF], and [START_REF] Schöpfer | Nonlinear iterative methods for linear ill-posed problems in Banach spaces[END_REF] Theorem 2.12) and with the Xu-Roach inequality (see [START_REF] Xu | Characteristic inequalities of uniformly convex and uniformly smooth Banach spaces[END_REF], and [START_REF] Schöpfer | Nonlinear iterative methods for linear ill-posed problems in Banach spaces[END_REF] Theorem 2.8), will have an important role in the convergence analysis of our algorithm.

CG and its convergence: the noise-free case

We restrict to the case X = L p with 1 < p < 2 and Y = L r , with r > 1. By implicitly considering the classical isomorphisms, in the following L q and L s denote the dual spaces of L p and L r , respectively. Concerning our functional equation ( 1), we are interested in finding the minimum p-norm solution of (1), denoted hereinafter as x, with exact data b ∈ R(A). In this respect, we first recall the following existence and characterization result.

Theorem 2 [START_REF] Schöpfer | Nonlinear iterative methods for linear ill-posed problems in Banach spaces[END_REF] If b ∈ R(A), then there exists the minimum p-norm solution x ∈ L p of (1) and

j p (x) ∈ R(A * ). Moreover, if x ∈ L p fulfils j p (x) ∈ R(A * ) and x -x ∈ N (A), then x = x.
To recover the minimum p-norm solution x, we propose the following algorithm based on the well known conjugate gradient method in Hilbert spaces.

Algorithm 1 Conjugate Gradient based method in Banach spaces

Choose a constant C ∈ (0, 1), and let

γ = Cr 2 r -1 + Cr (7)
be a relaxation parameter and d an arbitrary constant satisfying

0 < d ≤ " 1 - 2 r -1 + r r γ « 1 A . Set n = 0, x * 0 = 0, p * 0 = A * j r (b), R 0 = b , and
α 0 = arg min α∈ h 0, q p-1 A p R p-r 0 i Aj q (x * 0 + αp * 0 ) -b r . ( 8
) While R n > 0 do: Update n ← n + 1. Compute 8 < : x * n = x * n-1 + α n-1 p * n-1 , x n = j q (x * n ) , p * n = -A * j r (Ax n -b) + β n p * n-1 , (9)
where

α n = arg min α∈[0,T n ] Aj q (x * n + αp * n ) -b r (10)
and

β n = γ R r n R r n-1 , (11)
with

R n = Ax n -b , T n = min ( R 2-r n (V n -d A Q n ) G q 2 q-2 x * n q-2 (q -1) A 2 Q 2 n , x * n A R r-1 n Q n ) , V n = 1 - 2 r -1 r γQ n-1 , Q n = 1 -γ n+1 1 -γ ,
and G q the constant value of the Banach space defined in [START_REF] Schöpfer | Nonlinear iterative methods for linear ill-posed problems in Banach spaces[END_REF].

End while

Before providing the convergence proof, which follows the idea developed in [START_REF] Schöpfer | Nonlinear iterative methods for linear ill-posed problems in Banach spaces[END_REF], some comments are useful. First, it is interesting to notice that the functional p * n is a linear combination of the current steepest descent functional -A * j r (Ax nb) and the previous descent functional p * n-1 . This is a key point in the definition of the CG method in Hilbert spaces, where the descent functionals are just descent directions. On the other hand, there are two important facts to mention as main differences between (the proposed) CG in Banach and (the conventional) CG in Hilbert spaces. The first is that in Banach spaces the optimal step size α n cannot be directly computed by means of an explicit formula, since the corresponding one-dimensional minimization problem [START_REF] Grasmair | Generalized Bregman distances and convergences rate for non-convex regularization methods[END_REF] is no more quadratic. The second is related to the weight β n of [START_REF] Gratton | An active-set trust-region method for derivative-free nonlinear bound-constrained optimization[END_REF]. Now β n is different from the Fletcher Reeves formula for CG in Hilbert spaces. Indeed β n now concerns the ratio of the norms of the two last residuals (instead of the norms of the last two gradients) and it requires the relaxation factor γ < 1/2 defined in [START_REF] Cioranescu | Geometry of Banach spaces, duality mappings and nonlinear problems[END_REF], which does not appear in the original CG method. Hence, recalling the optimality condition on the step size [START_REF] Grasmair | Generalized Bregman distances and convergences rate for non-convex regularization methods[END_REF], differing from the original CG method, β n is always lower than 1/2. Although such a constraint is necessary only for our theoretical convergence proof, our numerical tests show that it helps the method to be more stable.

We prove the following main convergence result.

Theorem 3 Let X = L p with 1 < p < 2 and Y = L r , with r > 1.
The sequence of the iterations (x n ) n of Algorithm 1 defined in [START_REF] Brianzi | Preconditioned iterative regularization in Banach spaces[END_REF], either stops at or converges strongly to the minimum p-norm solution x of (1).

Proof First of all, by simple induction on the definition of x * n and p * n , we have that any

x * n = j p (x n ) ∈ R(A * ), since x * 0 = 0 ∈ R(A * ) and p * 0 = A * j r (b) ∈ R(A * ). Then, if the algorithm at step n gives R n = 0, by Theorem 2 we have that x n = x. Indeed, by 0 = R n = Ax n -b = A(x n -x) , we know that x n -x ∈ N (A), so that x n = x.
Otherwise we have that R n > 0 for all n > 0. In this case, as similarly done in [START_REF] Schöpfer | Nonlinear iterative methods for linear ill-posed problems in Banach spaces[END_REF] for the Landweber method in Banach space, the proof of convergence will be structured in the following four steps:

1. the sequence of the Bregman distances (∆ p (x, x n )) n obeys a recursive inequality which implies its convergence; 2. the sequence of the residuals (R n ) n is such that lim n-→+∞ R n = 0 ; 3. the sequence of the iterates (x n ) n has a Cauchy subsequence; 4. the sequence of the iterates (x n ) n converges strongly to x .

Before going on, we itemize some notations used in the proof.

g * n = -A * j r (Ax n -b), ∆ n = ∆ p (x, x n ), c k,n = n Y l=k+1 β l = γ n-k R r n R r k , τ n = α n A R r-1 n Q n x * n , W = G q 2 q-2 (q -1), S n = W τ n x * n q-1 R n Q n V n A .
1.

First step: the sequence of Bregman distances (∆ n ) n is nondecreasing and convergent. By (5), we can write

∆ n+1 = 1 p x p - 1 p x n+1 p + x * n+1 , x -x n+1 = 1 q x n+1 p + 1 p x p -x * n+1 , x = 1 q x * n+1 q + 1 p x p -x * n+1 , x = 1 q x * n + α n p * n q + 1 p x p -x * n , x -α n p * n , x .
For n = 0, since

x 0 = 0, then R 0 = b > 0, p * 0 = A * j r (b) and ∆ 0 = 1 p x p . Hence ∆ 1 = 1 q α 0 p * 0 q + 1 p x p -α 0 p * 0 , x = 1 q α q 0 A * j r (b) q + ∆ 0 -α 0 j r (b), Ax = 1 q α q 0 A * j r (b) q + ∆ 0 -α 0 j r (b), b ≤ 1 q α q 0 A q R (r-1)q 0 + ∆ 0 -α 0 R r 0
To obtain a reduction of the Bregman distance at first iteration, we consider the upper-bound for α 0 > 0 defined in (8), since

1 q A q α q 0 R (r-1)q 0 -α 0 R r 0 < 0 ⇐⇒ α 0 < q p-1 A p R p-r 0 =⇒ ∆ 1 < ∆ 0 . (12)
Now we consider n > 0. Thanks to the Xu-Roach inequality ([20], Theorem 2.8),

∆ n+1 = 1 q x * n + α n p * n q + 1 p x p -x * n , x -α n p * n , x ≤ 1 q h x * n q + qG q Z 1 0 ( x * n + tα n p * n ∨ x * n ) q t ̺ q " t α n p * n x * n + tα n p * n ∨ x * n « dt -q x n , -α n p * n i + 1 p x p -x * n , x -α n p * n , x ,
that is

∆ n+1 ≤ ∆ n + G q Z 1 0 ( x * n + tα n p * n ∨ x * n ) q t ̺ q " t α n p * n x * n + tα n p * n ∨ x * n « dt + α n p * n , x n -x . ( 13
)
We now study the last two addenda of the inequality [START_REF] Herzog | A conjugate direction method for linear systems in Banach spaces[END_REF].

A) Let us consider the term

p * n , x n -x . (14)
From

p * n = g * n + β n p * n-1 = g * n + n-1 X k=0 2 4 n Y l=k+1 β l 3 5 g * k = g * n + n-1 X k=0 c k,n g * k = -A * " j r (Ax n -b) + n-1 X k=0 c k,n j r (Ax k -b) # ,
thanks to (6) we have

p * n , x n -x = -j r (Ax n -b) + n-1 X k=0 c k,n j r (Ax k -b), Ax n -b = -R r n - n-1 X k=0 c k,n j r (Ax k -b), Ax n -b ≤ -R r n - n-1 X k=0 c k,n » 1 r Ax k -b r - 1 r Ax k -b -Ax n + b r - ≤ -R r n - n-1 X k=0 c k,n » 1 r R r k - 1 r " Ax k -b + Ax n -b " r - = -R r n + n-1 X k=0 c k,n » 1 r (R k + R n ) r - 1 r R r k - ≤ -R r n + 2 r -1 r n-1 X k=0 c k,n R r k , (15)
where the last inequality holds by construction of the optimal step size [START_REF] Grasmair | Generalized Bregman distances and convergences rate for non-convex regularization methods[END_REF] which guarantees that the residuals R k decrease.

B) Let us consider the term

G q Z 1 0 ( x * n + tα n p * n ∨ x * n ) q t ̺ q " t α n p * n x * n + tα n p * n ∨ x * n « dt . (16)
Since 0 ≤ t ≤ 1 and

x * n ≤ x * n + tα n p * n ∨ x * n ≤ x * n + α n p * n , we have p * n = -A * " j r (Ax n -b) + n-1 X k=0 c k,n j r (Ax k -b) # , so that p * n = ˛˛- A * " j r (Ax n -b) + n-1 X k=0 c k,n j r (Ax k -b) #˛≤ A ˛˛j r (Ax n -b) + n-1 X k=0 c k,n j r (Ax k -b) ˛≤ A " j r (Ax n -b) + n-1 X k=0 c k,n j r (Ax k -b) # = A " R r-1 n + n-1 X k=0 γ n-k R r n R r k R r-1 k # ≤ A R r-1 n " 1 + n-1 X k=0 γ n-k R n R k # ≤ A R r-1 n " 1 + n-1 X k=0 γ n-k # = A R r-1 n " 1 -γ n+1 1 -γ # = A R r-1 n Q n , (17)
that is,

p * n ≤ A R r-1 n Q n .
This last inequality, together with condition on T n of (10), gives

α n p * n ≤ α n A R r-1 n Q n ≤ x * n . ( 18
)
By means of the latter, we have that

Z 1 0 ( x * n + tα n p * n ∨ x * n ) q t ̺ q " t α n p * n x * n + tα n p * n ∨ x * n « dt ≤ Z 1 0 ( x * n + α n p * n ) q t ̺ q " t α n p * n x * n « dt ≤ Z 1 0 ( x * n + α n p * n ) q t ̺ q tα n A R r-1 n Q n x * n ! dt ≤ ` x * n + α n p * n ´q Z 1 0 1 t ̺ q (tτ n ) dt = ` x * n + α n p * n ´q Z τ n 0 1 t 1 ̺ q (t 1 ) dt 1 ≤ ` x * n + α n p * n ´q Z 1 0 1 t ̺ q (tτ n ) dt = ` x * n + α n p * n ´q Z τ n 0 1 t 1 ̺ q (t 1 ) dt 1 ≤ ` x * n + α n p * n ´q Z τ n 0 1 t 1 q -1 2 t 2 1 dt 1 = ` x * n + α n p * n ´q q -1 2 2 τ 2 n ≤ (2 x * n ) q (q -1) 2 2 τ 2 n = 2 q-2 x * n q-2 (q -1)α 2 n A 2 R 2(r-1) n Q 2 n . (19)
Now, reconsidering (13), we can write

∆ n+1 ≤ ∆ n + G q Z 1 0 ( x * n + tα n p * n ∨ x * n ) q t ̺ q " t α n p * n x * n + tα n p * n ∨ x * n « dt + α n p * n , x n -x ≤ ∆ n + G q h 2 q-2 x * n q-2 (q -1)α 2 n A 2 R 2(r-1) n Q 2 n i + α n " -R r n + 2 r -1 r n-1 X k=0 c k,n R r k # .
We have that (∆ n ) n is a non-increasing sequence, that is,

∆ n+1 ≤ ∆ n , if α n α n h G q 2 q-2 x * n q-2 (q -1) A 2 R 2(r-1) n Q 2 n i -R r n + 2 r -1 r n-1 X k=0 c k,n R r k ! ≤ 0. ( 20
) The latter holds if R r n - 2 r -1 r n-1 X k=0 c k,n R r k ≥ 0 ( 21
)
and

0 ≤ α n ≤ min ( R r n -2 r -1 r P n-1 k=0 c k,n R r k G q 2 q-2 x * n q-2 (q -1) A 2 R 2(r-1) n Q 2 n , x * n A R r-1 n Q n ) , (22)
∀n = 1, 2, 3, . . ., where condition [START_REF] Milicic | On moduli of expansion of the duality mapping of smooth Banach spaces[END_REF] has been considered. First we show that the relaxation parameter γ < 1 2 of ( 7) allows ( 21) to be always true. Indeed we have

n-1 X k=0 c k,n R r k = n-1 X k=0 γ n-k R r n R r k R r k = n-1 X k=0 γ n-k R r n = γ -γ n+1 1 -γ R r n ≤ γ 1 -γ R r n , so that R r n - 2 r -1 r n-1 X k=0 c k,n R r k ≥ " 1 - 2 r -1 r γ 1 -γ « R r n = (1 -C)R r n > 0 . ( 23
)
By means of these conditions on γ and α n , we can finally write

0 ≤ ∆ n+1 ≤ ∆ n + W x * n q τ 2 n - x * n τ n A Q n V n R n ≤ ∆ n , (24)
where τ n , W and V n have been previously properly defined. According to the aim of the first step of this proof, this last inequality shows that the sequence (∆ n ) n is both non-increasing and bounded by below, so that it is convergent.

Remark 1 From ( 15) and ( 23), we obtain thatp * n , x nx > 0, which shows that p * n is a descent functional at the current point x n , as expected by any gradienttype iteration scheme.

2.

Second step: the sequence of residuals (R n ) n vanishes, that is, lim n-→+∞ R n = 0. On the ground of [START_REF] Schöpfer | Nonlinear iterative methods for linear ill-posed problems in Banach spaces[END_REF] Theorem 2.12 (b), since by first step we know that (∆ n ) n is bounded, then the sequence (x n ) n is bounded. In addition, since x * n = j p (x n ) = x n p-1 and R n = Ax nb , then both the sequences (x * n ) n and (R n ) n are bounded too.

We rewrite [START_REF] Scheinberg | Self-correcting geometry in model based algorithms for derivative-free unconstrained optimization[END_REF] as follows

∆ n+1 ≤ ∆ n - 1 -S n A V n Q n τ n x * n R n ≤ ∆ n that is ∆ n -∆ n+1 ≥ 1 -S n A V n Q n τ n x * n R n ≥ 0 .
For all n ∈ N we have

0 ≤ n X k=0 1 -S k A V k Q k τ k x * k R k ≤ n X k=0 (∆ k -∆ k+1 ) = ∆ 0 -∆ n+1 ≤ ∆ 0 which yields that +∞ X k=0 1 -S k A V k Q k τ k x * k R k < +∞, (25)
that is, the numerical series is convergent. Since S k ≥ 0 by construction, in order to first ensure that the limit inferior of (R n ) n is zero we impose the following condition

0 ≤ S k ≤ 1 - d A Q k V k so that 1 -S k A V k Q k ≥ d > 0, ( 26
)
where d is a positive constant independent from n. This guarantees that the product of the first two factors of any element of the series (25) stays uniformly bounded away from zero. The latter inequality leads to the following constructive condition on the constant d

1 - d A Q k V k ≥ 0 which is satisfied when d is chosen such that 0 < d ≤ " 1 - 2 r -1 + r r γ « 1 A (27) since 0 < 1 - 2 r -1 + r r γ ≤ V k Q k . ( 28
)
On this ground, we can write

1 - d A Q k V k ≥ S k = W τ k x * k q-1 R k Q k V k A ,
which is equivalent to

τ k = α k A q R r-1 k Q k x * k ≤ (V k -d A Q k ) R k W x * k q-1 Q k A , that is, α k ≤ R 2-r k (V k -d A Q k ) W x * k q-2 Q 2 k A 2 , (29)
as requested by the the upper bound T n of [START_REF] Grasmair | Generalized Bregman distances and convergences rate for non-convex regularization methods[END_REF]. From ( 25), we finally have

d +∞ X k=0 τ k x * k R k ≤ n X k=0 1 -S k A V k Q k τ k x * k R k < +∞. ( 30
)
Suppose now that lim inf n-→+∞ R n > 0 . Then there exists n 0 ∈ N and ε > 0 such that R n ≥ ε ∀n ≥ n 0 . In this case, from (30), we would have

εd ∞ X k=n 0 τ k x * k ≤ d +∞ X k=n 0 τ k x * k R k < +∞.
Since the series is convergent, then we should have true at least one of the following two conditions:

1. (x k ) k is a null sequence. 1 p x p = ∆ p (0, x) = lim k-→+∞ ∆ p (x k , x) < 1 p x p p
But this is a contradiction, so that (x n ) n cannot be vanishing.

the numerical sequence (τ

k ) n is a null sequence. Since τ k = α k A R r-1 k Q k x * k , if τ k -→ 0 then α k -→ 0 because R k ≥ ε, and x *
k ≤ c. But the step size α k , which solves the one-variable minimization problem [START_REF] Grasmair | Generalized Bregman distances and convergences rate for non-convex regularization methods[END_REF], cannot be a vanishing sequence when the residual does not vanish, that is, if R k ≤ ε definitively, as supposed. Indeed, from ( 14) and ( 23), we have that

p * k , x k -x ≤ -(1 -C)R r k ≤ -(1 -C)ε r < 0 , which means that any functional p *
k is a descent functional with negative slope uniformly bounded away from zero. Moreover, since the duality map j q is continuous, then the map φ k : R -→ R defined as

α -→ φ k (α) = Aj q (x * k + αp * k ) -b r is continuous with respect to α, with φ k (0) = R r
k ≥ ε r and φ ′ k (0) ≤ c < 0, with c = c(ε) a fixed constant. On this grounds, by sign permanence Theorem, we have that there exists a δ > 0 such that φ n (α) is decreasing in [0, δ], so that α k ≥ δ > 0 for any k. So, we have necessarily that lim inf n-→+∞ R n = 0 . On the other hand, by construction of the step size based on the one dimensional minimization problem [START_REF] Grasmair | Generalized Bregman distances and convergences rate for non-convex regularization methods[END_REF], we recall that 0

< R n+1 ≤ R n , so that lim n-→+∞ R n = 0. ( 31
)
3. Third step: the sequence (x n ) n has a Cauchy subsequence. By boundedness of (x n ) n and (j p (x n )) n we can find a subsequence (x n k ) k such that (a) the sequence of the norm ( x n k ) k is convergent, (b) the sequence (j p (x n k )) k is weakly convergent, The sequence (R n k ) k is also null by (31). We want to show that (x n k ) k is a Cauchy (sub)sequence. By virtue of [START_REF] Schöpfer | Nonlinear iterative methods for linear ill-posed problems in Banach spaces[END_REF] Theorem 2.12 (e), we have for all k, l ∈ N, with l > k,

∆ p (x n k , x n l ) = 1 q ( x n k p -x n l p ) + j p (x n l ) -j p (x n k ), x n l .
The first addendum satisfies [START_REF] Schöpfer | Nonlinear iterative methods for linear ill-posed problems in Banach spaces[END_REF] Theorem 2.12 (d) (ii).

( x n k p -x n l p ) -→ 0 as k -→ +∞ for
Regarding the second term, we have

j p (x n l ) -j p (x n k ), x n l = j p (x n l ) -j p (x n k ), x + j p (x n l ) -j p (x n k ), x n l -x .
(32) Here j p (x n l )j p (x n k ), x -→ 0 as k -→ +∞ for (b), and

| j p (x n l ) -j p (x n k ), x n l -x | = ˛nl -1 X n=n k j p (x n+1 ) -j p (x n ), x n l -x = ˛nl -1 X n=n k α n p * n , x n l -x ≤ ˛nl -1 X n=n k α n ' ' ' ' ' ' j r (Ax n -b) + n-1 X j=0 c j,n j r (Ax j -b) ' ' ' ' ' ' Ax n l -b ≤ n l -1 X n=n k α n R r-1 n Q n R n = n l -1 X n=n k τ n x * n R n A → 0 as k → +∞ by (30).
Hence we can write

j p (x n l ) -j p (x n k ), x n l -→ 0,
which prove that the (sub)sequence (x n k ) k is a Cauchy sequence. Finally, since it is Cauchy, it converges strongly to a point x ∈ X .

4.

Forth step: the sequence (x n ) n converges strongly to the minimum p-norm solution x We know that the subsequence (x n k ) n k converges to an element x ∈ X . We prove now that x = x and lim n-→+∞

x n -x = 0. From Ax n k -b = A (x n k -x) = R n k -→ 0 as k -→ +∞ ,
then, by continuity of A, we have that (xx) ∈ N (A) and j p (x) ∈ R(A * ). Hence, by recalling Theorem 2, x = x. It remains to prove that not only the subsequent (x n k ) k converges to x, but the full sequence (x n ) n converges strongly to x. By continuity of the Bregman distance, we can state that lim k-→+∞

∆ n k = lim k-→+∞ ∆ p (x, x n k ) = ∆ p (x, lim k-→+∞ x n k ) = ∆ p (x, x) = 0.
Hence the sequence (∆ n ) n is convergent (by first point of the proof), and has a subsequence converging to zero. Theorem 2.12 (d) of [START_REF] Schöpfer | Nonlinear iterative methods for linear ill-posed problems in Banach spaces[END_REF], allows us to state that x n -→ x strongly, that is,

x nx -→ 0 as n -→ +∞ , which concludes the proof.

CG and its regularization properties: the noisy case

In the previous section, we proved that the CG method, in the case X = L p with 1 < p < 2 and Y = L r , converges to the minimum p-norm solution of (1) for exact data b ∈ R(A). In the same setting, now we consider the case of noisy data b δ ∈ Y with a known noise level δ > 0, that is,

b -b δ ≤ δ.
According to the well known semi-convergence behavior of any iterative regularization algorithm, in the case of noisy data an early stop of the iterations prevents noise amplification in the restoration process. Hence, our algorithm for noisy data is based on the discrepancy principle stopping rule, that is, the iterations are stopped as soon as Ax nb δ ≤ τ δ , where τ > 1 is a fixed constant value. This way, the restoration is accepted as soon as its residual is equal or less than τ times the magnitude of the noise of the (noisy) data b δ (indeed, searching for a solution x n such that the corresponding residual is smaller than the noise level δ is unreasonable). After introducing the algorithm, in the subsequent theorem we will prove that it belongs to the class of iterative regularization methods, by showing that its iterations either stop at or converge strongly to the minimum p-norm solution x of ( 1) as the noise level δ goes to zero.

Algorithm 2 Regularization by Conjugate Gradient method in Banach spaces

Choose two constants C ∈ (0, 1) and τ ∈ (1, +∞) such that C + τ -1 ∈ (0, 1), and let

γ = Cr (2 + τ -1 ) r -1 + Cr (33)
be a relaxation parameter and d an arbitrary constant satisfying

0 < d ≤ (1 -γ)(1 -τ -1 ) - (2 + τ -1 ) r -1 r γ ! 1 A . Set n = 0, x * 0 = 0, p * 0 = A * j r (b δ ), R 0 = b δ , and
α 0 = arg min α∈ » 0, q p-1 (1-τ -1 ) p-1 R p-r 0 A p -Aj q (x * 0 + αp * 0 ) -b δ r , ( 34
) While R n > τ δ do: Update n ← n + 1. Compute < : x * n = x * n-1 + α n-1 p * n-1 , x n = j q (x * n ) , p * n = -A * j r (Ax n -b δ ) + β n p * n-1 , (35)
where

α n = arg min α∈[0,T n ] Aj q (x * n + αp * n ) -b δ r (36)
and

β n = γ R r n R r n-1 , (37)
with

R n = Ax n -b δ , T n = min ( R 2-r n (V n -d A Q n ) G q 2 q-2 x * n q-2 (q -1) A 2 Q 2 n , x * n A R r-1 n Q n ) , V n = 1 -τ -1 - (2 + τ -1 ) r -1 r γQ n-1 , Q n = 1 -γ n+1 1 -γ ,
and G q the constant value of the Banach space defined in [START_REF] Schöpfer | Nonlinear iterative methods for linear ill-posed problems in Banach spaces[END_REF]. End while Theorem 4 Let X = L p with 1 < p < 2 and Y = L r , with r > 1. The sequence of the iterations (x n ) n of Algorithm 2 is a regularization method.

Proof According to [START_REF] Schöpfer | Nonlinear iterative methods for linear ill-posed problems in Banach spaces[END_REF], to show that Algorithm 2 is a regularization method we consider a sequence of "less and less" noisy data (b h ) h∈N such that bb h ≤ δ h , with lim h-→+∞ δ h = 0, and 0 < δ h+1 < δ h . We then will prove that for any data δ h there exists an iteration x n such that Ax nb h < τ δ h , and that the sequence of iterates (x n ) n either stops at or converges strongly to the minimum p-norm solution x of (1) as the noise level δ h goes to zero. To this aim, we consider the following adaptive setting of the noisy data (b h ) h as the iterations go on.

REPEAT, for n = 0, 1, 2, . . ., the following steps:

IF for all h > h n-1 (where h -1 = 0) the discrepancy principle is satisfied, that is, Ax n -b h < τ δ h , (38)
THEN STOP the iterations; ELSE choose the smallest integer value

h n > h n-1 such that Ax n -b h n ≥ τ δ h n , and consider Algorithm 2 with noisy data b δ = b h n , that is 8 < : x * n = x * n-1 + α n-1 p * n-1 x n = j q (x * n ) p * n = -A * j r (Ax n -b h n ) + β n p * n-1 (39)
We notice that the discrepancy principle related to δ h is satisfied for any h n < h ≤ h n-1 , and the algorithm find a regularized solution for the associated data b h with h n < h ≤ h n-1 . On the other hand, if the inequality (38) holds for all h > h n-1 at a certain iteration n, then x n = x as in the case of exact data (indeed the sequence of noise levels is such that lim h-→+∞ δ h = 0), that is, the iterations stop at the minimum p-norm solution of (1).

So we have to consider only the case where inequality (38) is never fulfilled. In such a case R n = Ax nb h n > 0 for all n. The proof follows the same steps of the noiseless case of Theorem 3, so we report the basic different arguments.

First we show that, since R n = Ax nb h n > τ δ h n , the sequence of Bregman distances (∆ n ) n is non-decreasing and convergent.

As already shown at the beginning of the first step of the proof of Theorem 3, we have

∆ n+1 = 1 q x * n + α n p * n q + 1 p x p -x * n , x -α n p * n , x .
For n = 0, it becomes

∆ 1 = 1 q α q 0 A * j r (b h 0 ) q + 1 p x p -α 0 j r (b h 0 ), Ax = 1 q α q 0 A * j r (b h 0 ) q + ∆ 0 -α 0 j r (b h 0 ), b h 0 + α 0 j r (b h 0 ), b h 0 -b ≤ 1 q α q 0 A q R (r-1)q 0 + ∆ 0 -α 0 R r 0 + α 0 R r-1 0 δ h 0 .
Similarly to [START_REF] Hansen | Rank-deficient and discrete ill-posed problem: numerical aspects of linear inversions[END_REF], this gives the upper-bound for α 0 > 0, in order to obtain a reduction of the Bregman distance at the first iteration. Indeed, since R 0 > τ δ h 0 we have that, if

1 q α q 0 A q R (r-1)q 0 -α 0 R r 0 + α 0 R r 0 τ -1 < 0 ⇐⇒ α 0 < q p-1 `1 -τ -1 ´p-1 R p-r 0 A p then =⇒ ∆ 1 < ∆ 0 .
For n > 0, we have again ( 13)

∆ n+1 ≤ ∆ n + G q Z 1 0 ( x * n + tα n p * n ∨ x * n ) q t ̺ q " t α n p * n x * n + tα n p * n ∨ x * n q « dt + α n p * n , x n -x .
We study the two latter addenda of this inequality.

A) Following the analogous computation of (15), since R k > τ δ h k we have

p n , x n -x = -A * " j r (Ax n -b h n ) + n-1 X k=0 c k,n j r (Ax k -b h k ) " , x n -x = -j r (Ax n -b h n ), Ax n -b h n + j r (Ax n -b h n ), b -b h n - n-1 X k=0 c k,n j r (Ax k -b h k ) , Ax n -b ≤ -R r n + R r-1 n δ h n - n-1 X k=0 c k,n » 1 r Ax k -b h k r - 1 r Ax k -b h k -Ax n + b r - ≤ -R r n + τ -1 R r n - n-1 X k=0 c k,n » 1 r R r k - 1 r ` Ax k -b h k + Ax n -b h n + δ h n ´r- ≤ (τ -1 -1)R r n + n-1 X k=0 c k,n » 1 r `Rk + (1 + τ -1 )R n ´r - 1 r R r k - ≤ (τ -1 -1)R r n + (2 + τ -1 ) r -1 r n-1 X k=0 c k,n R r k . (40)
B) The other term (which is equivalent to ( 16)) is estimated as in [START_REF] Lindestrauss | Classical Banach spaces[END_REF], where instead of j r (Ax nb) and j r (Ax kb), we respectively consider j r (Ax nb h n ) and j r (Ax kb h k ), obtaining again the same bound [START_REF] Milicic | On moduli of expansion of the duality mapping of smooth Banach spaces[END_REF] and the same inequality [START_REF] Resmerita | Regularization of ill-posed problems in Banach spaces: convergence rates[END_REF]. Now reconsidering (13), we have

∆ n+1 ≤ ∆ n + G q h 2 q-2 x * n q-2 (q -1)α 2 n R 2(r-1) n Q 2 n i + α n " (τ -1 -1)R r n + (2 + τ -1 ) r -1 r n-1 X k=0 c k,n R r k # , so that ∆ n+1 ≤ ∆ n if G q h 2 q-2 x * n (q-2) α 2 n R 2(p-1) n Q 2 n i + α n " (τ -1 -1)R r n + (2 + τ -1 ) r -1 r n-1 X k=0 c k,n R r k # ≤ 0 . (41)
The relaxation parameter (33) , as shown by the same computations of ( 23), now leads to

(1 -τ -1 )R r n - (2 + τ -1 ) r -1 r n-1 X k=0 c k,n R r k ≥ 1 -(C + τ -1 ) > 0 , (42)
since C + τ -1 ∈ (0, 1) by definition, so that (41) holds if

0 ≤ α n ≤ min 8 < : (1 -τ -1 )R r n -(2+τ -1 ) r -1 r P n-1 k=0 c k,n R r k G q 2 q-2 x * n q-2 (q -1) A 2 R 2(r-1) n Q 2 n , x * n A R r-1 n Q n 9 = ; , (43)
∀n = 1, 2, 3, . . ., where condition [START_REF] Milicic | On moduli of expansion of the duality mapping of smooth Banach spaces[END_REF] has been considered. Thanks to the same results of ( 24), the sequence (∆ n ) n is convergent. The second step, as well as for the noiseless case of Theorem 3, is to show that lim n-→+∞ R n = 0. Following exactly the same formal computation, the series [START_REF] Schuster | Regularization methods in Banach spaces[END_REF] still converges. Similarly to ( 27) and (28), the constant d of ( 26) is now chosen such that

0 < d ≤ (1 -γ)(1 -τ -1 ) - (2 + τ -1 ) r -1 r γ ! 1 A , since now 0 < (1 -γ)(1 -τ -1 ) - (2 + τ -1 ) r -1 r γ ≤ V k Q k .
This leads to the same bound (29) on the step size α n and the inequalities (30). On this ground, the reductio ad absurdum leading to (31) can be done similarly, by noticing that, thanks to (40) and (42),

p * k , x k -x ≤ -(1 -(C + τ -1 ))R r k ≤ -(1 -(C + τ -1 ))ε r < 0 ,
that is, p * k would still be a descent functional with negative slope uniformly bounded away from zero.

The third step of the proof is to show that the sequence (x n ) n has a Cauchy subsequence. The proof is basically the same, where the convergence of the second addendum of (32) is given by

| j p (x n l ) -j p (x n k ), x n l -x | = ˛nl -1 X n=n k j p (x n+1 ) -j p (x n ), x n l -x = ˛nl -1 X n=n k α n p * n , x n l -x = ˛nl -1 X n=n k α n j r (Ax n -b h n ) + n-1 X j=0 c j,n j r (Ax j -b h j ), Ax n l -b ≤ ˛nl -1 X n=n k α n j r (Ax n -b h n ) + n-1 X j=0 c j,n j r (Ax j -b h j ), Ax n l -b h n l + ˛nl -1 X n=n k α n j r (Ax n -b h n ) + n-1 X j=0 c j,n j r (Ax j -b h j ), b h n l -b ≤ n l -1 X n=n k α n R r-1 n Q n (R n l + δ h n l ) ≤ n l -1 X n=n k α n R r-1 n Q n (R n + δ h n ) ≤ (1 + τ -1 ) n l -1 X n=n k τ n x * n R n A → 0 as k → +∞ .
This ground, the proof can be completed as well as in the step four of the noiseless case.

As already sketched, the classical CG method in Hilbert space gives a simple and closed form expression for the step size α n of [START_REF] Grasmair | Generalized Bregman distances and convergences rate for non-convex regularization methods[END_REF], which is the same optimal choice of the steepest descent method [START_REF] Bertero | Introduction to inverse problems in imaging[END_REF]. Indeed, in that Hilbertian case, the corresponding minimization problem arg min

α A (x n + αp n ) -b 2
is a simple one-dimensional quadratic and differentiable minimization problem, whose first derivative is linear. In our more general Banach settings, problem (10)

 arg min α Aj q (x * n + αp * n ) -b r = arg min α Φ(α) α ∈ [0, T ]
does not allow a closed form explicit solution, according to Galois' results on the solution of a polynomial equation. Thus, the description of the practical implementation of the proposed algorithm leaves a practical question open.

In order to find the minimum of Φ(α), we can evaluate the first derivative of Φ(α)

Φ ′ (α) = j r `Aj q `x * k + αp * k ´-b ´ * • `(q -1) A `jq-1 (|x * k + αp * k |) • p * k ´´(44)
where • is the component-wise product, and we find the critical point by solving the equation Φ ′ (α) = 0. If p = 2 or r = 2, (44) is no longer a linear equation, so we cannot compute its solution explicitly. Iterative methods can be used, such as the simplest secant method

α k+1 = α k -Φ ′ (α k ) α k -α k-1 Φ ′ (α k ) -Φ ′ (α k-1 )
, k = 1, 2, ... as first basic example. However, this approach becomes very expensive in computational load in real problems with large-scale matrices.

To avoid this high computational load, we use the recent and efficient derivativefree algorithm for bound-constrained optimization (BCDFO) developed in [START_REF] Gratton | An active-set trust-region method for derivative-free nonlinear bound-constrained optimization[END_REF], which is based on the iterative trust-region method. The idea behind a trust region method is very simple: the k + 1 iteration α k+1 = α k + s k is an approximate solution of the quadratic subproblem

m k (α k + s k ) = Φ(α k ) + g T k s k + 1 2 s T k H k s, (45)
inside the trust region

B ∞ (α k , ∆ k ) = {α ∈ R| α -α k ∞ ≤ ∆ k } .
where g k is an approximation of the gradient of Φ at point α k , H k is an approximation of its Hessian and ∆ k is a positive scalar representing the size of the trust region.

The point α k + s k is acceptable and ∆ k is increased if

̺ k = Φ(α k ) -Φ(α k + s k ) m k (α k ) -m k (α k + s k ) > η ,
where η > 0 is a suitable constant. If ̺ k < η 1 , the size of the trust region ∆ k is decreased and s k = 0. The key fact is that the function m k is evaluated not using the true derivatives but interpolating known function's values at a given set Y k . The geometry of this space has to cover the space to guarantee the convergence from an arbitrary starting point [START_REF] Scheinberg | Self-correcting geometry in model based algorithms for derivative-free unconstrained optimization[END_REF]. In our case, the BCDFO is a particular implementation where bounds on the variables are handled by an active set strategy (for details and convergence proofs, see [START_REF] Gratton | An active-set trust-region method for derivative-free nonlinear bound-constrained optimization[END_REF]). Although efficient, the algorithm BCDFO to solve the one dimensional minimization problem (10) cannot be considered as optimal, so that the appropriate computation of the step size still remains an open problem.

Numerical results

In this section, first numerical experiments illustrate the performance of the proposed conjugate gradient technique in Banach spaces. We choose a basic image restoration problem [START_REF] Bertero | Introduction to inverse problems in imaging[END_REF] and we compare the technique with the Landweber method in Banach spaces (2). Consider the operator equation of the first kind

g(t 1 , t 2 ) = Z Z G 1 (y 1 , t 1 )f (y 1 , y 2 )G 2 (y 2 , t 2 )dy 1 dy 2 , (46)
where sions and different intensities. Since the solution is almost-sparse, we look for the minimum p-norm solution with p ∈ (1, 2), in order to promote sparsity in the restored solution. For more information about the role of the sparsity in inverse problems, we refer to [START_REF] Schuster | Regularization methods in Banach spaces[END_REF], [START_REF] Daubachies | An iterative thresholding algorithm for linear inverse problems with sparsity constraint[END_REF], [START_REF] Lenti | Analysis of reconstructions obtained solving l p penalized minimization problems[END_REF]. To simulate noisy data b δ , we add white Gaussian noise and measure the noise level bb δ = δ. When the noise is Gaussian, it is well known that the least square solution is the maximum likelihood estimate of x. For this reason, Axb δ 2 2 is used. In our numerical tests, the parameters of the proposed CG algorithm are selected as follows: γ = Cr (2+τ -1 ) r -1+Cr , with C = 0.9 and τ = 10.

G 1 (y 1 , t 1 ) = 1 √ 2πσ 1 e - (y 1 -t 1 ) 2 2σ 2 1 and G 2 (y 2 , t 2 ) = 1 √ 2πσ 2 e - (y 2 -t 2 ) 2 2σ 2
d = d ≤ " (1 -γ)(1 -τ -1 ) -(2+τ -1 ) r -1 r γ " 1 A -10 -6
, G q is evaluated using the estimate given in [START_REF] Schöpfer | Nonlinear iterative methods for linear ill-posed problems in Banach spaces[END_REF]. Fig. 2 shows the reconstructed solution by using the conjugate gradient method in l 2 (a), and by using the proposed approach in l 1.8 (b), in l 1.5 (c) and in l 1.2 (d). As stopping criterion, the discrepancy principle is used [START_REF] Hansen | Rank-deficient and discrete ill-posed problem: numerical aspects of linear inversions[END_REF]. By direct visual inspection, it is clear that smaller p promotes sparsity and reconstruction of abrupt discontinuities. The reconstructed solutions in Banach spaces are less over-smoothed than the typical of Hilbert reconstructions, and few artifacts and oscillations are present, both in the spot images and in the background. where X ref and X rec are the reference and the reconstructed solutions of Fig. 2, respectively. The RMSEs of the proposed CG and the Landweber methods are very similar, and the corresponding reconstructions look almost identical (so that, to save space, we do not show the reconstructions of the Landweber method).

Anyway, as we can see in Fig. 3, the proposed CG method considerably outperforms the Landweber one both in terms of iterations number and computational time, for each value of γ (7). We recall that γ = Cr (2+τ -1 ) r -1+Cr , where c ∈ [0, 1) , then, in this case, γ ∈ [0, 0.4) for r = 2. It must be pointed out that the computational cost of each iteration is completely different for Landweber and the proposed technique. In the latter, a 1-dimensional non-linear problem [START_REF] Grasmair | Generalized Bregman distances and convergences rate for non-convex regularization methods[END_REF] must be minimized at each iteration. The dependence on the γ and p parameters in terms of computational time is Fig. 3 Iteration number(a) and computational time (b) against the p parameter for Conjugate gradient and Landweber method. Note that the conjugate gradient method with γ = 0 is the steepest descent one.

analyzed. In fig. 4, the iteration numbers and computational time against c and p are shown. The computational time decreases for increasing p and γ values. Like in Hilbert space, the conjugate gradient method in l p outperforms the Landweber method and the steepest descent method. These results are confirmed by different experiments (not shown here to save space) and an application of the proposed technique to a real remote sensed data could be found in [START_REF] Lenti | Conjugate gradient method in Hilbert and Banach spaces to enhance the spatial resolution of radiometer data[END_REF]. On these grounds, we can summarize that the conjugate gradient method is really more efficient of the Landweber method especially for large scale matrices . Moreover, its convergence speed can be further improved by means of the applications of a dual preconditioner as in [START_REF] Brianzi | Preconditioned iterative regularization in Banach spaces[END_REF].

Conclusions

An iterative method, based on a generalization of the conjugate gradient method for the minimization in l p Banach spaces is proposed. We demonstrate that the conjugate gradient method converges strongly to the minimum p-norm solution and that together with the discrepancy principle as stopping rule, the proposed method is a regularization method. Numerical experiments, undertaken on an image restoration problem, show that the method is robust in terms of reconstruction accuracy and it is faster than the Landweber method and the steepest descent one.

2 .

 2 The linear system Ax = b obtained by the discretization of the operator equation (46) is considered. Fig.1(a) shows the exact solution x † for the given right-hand side b shown in Fig.1(b). The exact solution x † presents discontinuities characterized by different dimen-

Fig. 1

 1 Fig. 1 Exact solution x † of the linear system Ax = b (a), exact data b (b).

Fig. 2

 2 Fig.2Reconstructed solution using conjugate gradient in Hilbert space (a) and using the generalization of the conjugate gradient in l 1.8 (b), l 1.5 (c) and l 1.2 (d).

Fig. 4

 4 Fig.4Iteration number (a) and computational time (b) against the c and p parameter for Conjugate gradient method. Note that the conjugate gradient method with γ = 0 is the steepest descent one.

Table 1

 1 RMSE for CG and Landweber methods in l pTable1shows the Root Mean Square Error (RMSE):RM SE(X rec , X ref ) = X rec -X ref

	Method	p = 2	p = 1.8	p = 1.5	p = 1.2
	CG	0.4091	0.3568	0.3246	0.2454
	LW	0.3957	0.3552	0.3311	0.2452
				2 X ref 2	2 2	(47)

Acknowledgements The work of C. Estatico is partly supported by PRIN 2012 2012MTE38N and GNCS-INdAM. The work of F. Lenti is supported by the AVENUE project of the RTRA STAE foundation.


Wilfried Uhring
email: wilfried.uhring@unistra.fr

Laurent Millet

Bertrand Misischi

Fatah Rarbi

Fabrice Guellec

Daniel Dzahini

Octavian Maciu

Jean-Baptiste Kammerer

Gilles Sicard

A Scalable Architecture for Multi Millions Frames per Second CMOS Sensor With Digital Storage

Keywords: CMOS video sensor, ultra fast burst imaging, digital storage, CMOS image sensor, 3D stacking I

This paper describes a 3D Integrated Circuit (3D-IC) architecture of a burst image sensor (BIS) with embedded digitization and digital storage. This architecture also proposes a new technique to further increase both the frame rate and the stored image capacity at the cost of a spatial resolution reduction. A 2D monolithic demonstrator that takes into account the constraints of a future 3D-IC imager has been fabricated. Experimental results are presented showing that a frame rate from 5 up to 45 Mega frames per second can be achieved. This fully functional approach paves the way to the very first in-focalplane digital BIS.

INTRODUCTION

Conventional video systems extract data from the image sensor at each frame and store it into an external memory for post processing or display. The constrained number of Input/Output (I/O) pads of the sensors limits the data rate and consequently the pixel readout rate to about 25 Gpixel/s for the fastest currently available video camera [1]. This bottleneck results in the limitation of the frame rate to a few kfps (kilo frame per second) when the frame resolution is close to 1 million pixel or about 1 Mfps (Mega frames per second) when the frame resolution is reduced to about 1 kpixel. If a higher pixel rate is required, the I/O issue can be circumvented if the data is stored in the sensor at a very high data rate during the acquisition and readout afterward at a relaxed speed.

This approach is known as Burst Image Sensor (BIS) and consists to use an on chip memory to store the video stream. At this time, published BIS use analog CMOS capacitances [START_REF] Zlatanski | Sub-500 ps temporal resolution streak-mode optical sensor[END_REF] [START_REF] Kuroda | A 20Mfps Global Shutter CMOS Image Sensor with Improved Light Sensitivity and Power Consumption Performances[END_REF] or CCD channel memories [START_REF] Crooks | Kirana: a solid-state megapixel uCMOS image sensor for ultrahigh speed imaging[END_REF] [START_REF] Etoh | Toward 1Gfps: Evolution of ultra-high-speed image sensors -ISIS, BSI, multi-collection gates, and 3D-stacking[END_REF] and are designed in 2D monolithic integrated circuit (IC).

A BIS that stores the recorded data frames within its chip can benefit from a digital storage as the information density can be much higher in the digital domain compared to the analog domain. For instance, in a 28 nm digital process, the surface of respectively a static and a dynamic bit cell are 0.12 µm² and 0.035 µm². As a consequence, while analog BIS can hardly store more than 200 frames in in-focal-plane storage topologies, a digital BIS can easily record respectively more than 1000 or 4000 frames by using a SRAM or a DRAM based memory bank and by assuming a pixel pitch of 50x50 µm and 8 bits dynamic [START_REF] Bonnard | New 3Dintegrated burst image sensor architectures with in-situ A/D conversion[END_REF]. On the other hand, the surface occupied by the Analog to Digital Converter (ADC) within the pixel will reduce the available area to the memory in a single 2D IC sensor. A 3D stacked IC approach overcomes this limitation by allowing a distribution of electronic functions on dedicated tiers. Top tier photodiodes can benefit from imagery technology and back side illumination for a maximal sensitivity, the ADC gains area and partitioning design flexibility on middle tier, and in the bottom tier almost all the pixel area is available to store the frames thanks to advanced digital node technology.

The paper describes a 3D-IC architecture of a scalable design for burst image sensors with on the fly digitization and embedded digital storage applied to a 3D stack of 3 integrated circuits. In order to demonstrate the benefits of this approach, a 2D monolithic circuit integrating the different electronic parts of the sensor side by side has been designed in a 65 nm CMOS technology. The overall sensor architecture is presented in section II and its different parts are described in section III to VI. The section VII presents the demonstrator and experimental results are shown in section VIII.

II. OVERALL SENSOR ARCHITECTURE

The 3D design allows a better repartition of the electronic parts and the targeted cross-section of the sensor is presented on fig. 1. A first IC is dedicated to the analog pixel array. A second one is devoted to the ADC and the last one is dedicated to the embedded digital memory.

A straight forward implementation of this architecture is to store an ADC and its associated memory in each pixel. This topology is not necessarily optimal in term of frame rate and sequence depth, i.e. the number of storable images. An on-thefly analog to digital conversion at a targeted frame rate of several tens of Mfps requires a high speed ADC which results in a high silicon area and power consumption as explained in section IV. Likewise, segmenting the memory floorplan into a unit for each pixel is not optimal in term of memory capacity since space is lost for addressing blocks and boundary separations.

For these reasons, a cluster approach has been adopted for the BIS. This line has been recently used with a 2 layers stacked global-shutter image sensor [START_REF] Kondo | A 3D stacked global-shutter image sensor with pixellevel interconnection technology for high-speed image capturing[END_REF]. The detailed architecture of the sensor according to the different tiers is shown on fig. 2. A cluster of pixels is composed of P=10 frontend pixels connected to one ADC. A bloc of L clusters, i.e. LxP pixels, shares the same digital memory for an optimal floorplan configuration that allows embedding the maximal amount of memory.

Furthermore this cluster of pixels approach sharing the same memory unit brings memory addressing flexibility since this feature allows to digitize an adaptable number of pixels. In "full frame" mode an ADC operating at a sampling rate Fs and reading all the P pixels gives an achievable frame rate of Fs/P. But we implemented some modes to improve the frame rate at the cost of a spatial resolution reduction. By assuming a constant data conversion rate, the frame rate increases as the number of read pixels in the cluster is reduced and is given by (1)

RC RC Fs P Fs Fps P P P = ⋅ = (1)
Where P is the number of pixels in the cluster and PRC∈[1;P] is the number of read pixels within the cluster. Meanwhile, the movie length is also increased by a factor P/PRC. Increasing the frame rate is obviously very interesting in high speed imaging to detect faster events whereas expanding the movie length allows the observations of fast and complex event that last for a longer time. While using this feature, the ADC and memory efforts are kept constant whereas the pixel front-end operates at a higher speed. Thus the front-end bandwidth has to be high enough to operate at the highest frame rate or it has to be increased as the frame rate increases to keep the sensor sensitivity constant.

III. ANALOG FRONT-END DESIGN

In high speed imaging the analog front-end has to operate at a frame rate of up to several tens of MHz. As a consequence, the exposure time is much less than 1 µs and thus the integrated light energy is potentially very low. Therefore, the front-end must provide both a high bandwidth and a high conversion gain. A Capacitive Trans-Impedance Amplifier (CTIA) provides a very low impedance node, a fixed bias and a very high conversion gain as well as a high bandwidth. The main drawback of a CTIA is the area and the power consumption which both are not critical for this design as the pixel pitch is 50 µm and as a backside illumination is expected for the final sensor. For all this reasons, a CTIA with a transconductance stage based on a cascode differential pair with 2.5V power supply voltage has been implemented in the pixel. In the 2D demonstrator the photodetector is a 35×40 µm PSub-Nwell placed beside the CTIA and a buffered analog multiplexer that allows to select one of the pixel output to the shared ADC input. A sample and hold circuitry with a ping pong structure to store the even and odd frames makes possible a pipeline integration, i.e. the odd frames are digitized during the integration of the even frames and vice versa. The CTIA integrating capacitor is 27 fF and the biasing current is 20 μA. An open loop gain of 40 dB and a bandwidth of 30 MHz are achieved. The 10 pixels cluster layout is shown in Fig. 4.

IV. ANALOG TO DIGITAL CONVERTER

In the case of a digital BIS the main bottleneck is the ADC. As the cluster approach allows to use a single ADC shared for several pixels, it can occupy the corresponding silicon surface minus the space required by some Trough Silicon Via (TSV) mandatory to convey signals from the top tier.

Regarding the state of the art from the Murmann data base [START_REF] Murmann | ADC Performance Survey[END_REF] (Fig. 3), the successive approximation register (SAR) and pipeline architectures seem to be two good candidates for the targeted frame rates between 10 and 100 Mfps with an 8 bits resolution. Pipeline architecture is more subject to missing code error than SAR architecture and delta sigma are too slow for our application. Moreover, the SAR ADC architecture disposes Fig. 2. IC repartition of a 3D CMOS stacked BIS. The first tier is for the pixel frontend, the second tier shares one ADC for one pixel clusters and the third tier is for the deigital memory shared for several ADC clusters. Fig. 3. State of the art of 8-bit ADC from [START_REF] Murmann | ADC Performance Survey[END_REF]. At the targeted frame rate, SAR architecture offers the best signal to noise ratio. intrinsically of an interesting single bit stream output that limit the number of required TSV to transmit the converted data. Therefore, this architecture is preferred for this range of frame rate. The designed SAR ADC is a single-ended architecture based on a feedforward principle, and is described in detail in [START_REF] Rarbi | An 8-Bit, 100-MSPS Fully Dynamic SAR ADC for Ultra-High Speed Image Sensor[END_REF]. A binary capacitor array segmented by two reduces the total value of the capacitors seen by the input nodes during the sampling phase and also to reduce the total area of the chip. The core's size is about 210×100 μm 2 so that two ADCs fit the area of two clusters, i.e. 500×100 µm², and a space of more than 6000 µm² is available for the TSV required to the 3D stacked version of the final sensor. The post-layout noise simulation shows an effective number of bit (ENOB) of 7.5-b. The total power consumption of this SAR ADC is only 1.55 mW at 100 MS/s.

V. DIGITAL MEMORY

In the final 3D stacked sensor, the digital memory should be designed with an advanced CMOS process such as a 28 nm technology. In our monolithic demonstrator realized in a 65 nm process, a not optimized SRAM memory bank has been used. The memory receives the single bit stream provided by the ADCs and stores each pixel value into a 16384 bits large memory bank for each cluster. By assuming an 8 bits dynamic, the storable movie length is then 16384/8/10 = 204 frames long.

VI. SENSOR CONTROL

The sensor operation is controlled with 4 embedded state machines (SM) that allow to precisely control the fast different unit. These SMs are clocked by a 100 MHz to 1 GHz master clock provided by an external PLL and almost all the timing of these SM are configurable by writing to their specific register. The main SM waits for a "start" signal to initiate the record of a sequence. The integration SM controls the reset signal at each new frame. The write SM controls the ADC operation, in particular the sampling signal and the number of digit of the conversion in the range from 1 to 8 bits. Finally, the pixel read SM allows to choose the pixel to read into the cluster at each frame. Several read modes are available. The first mode, the so called "full frame" mode, consists to read all pixels. The second one, reads out only half of the cluster leading to a 2 × higher frame rate. The third one reads out only 2 pixels of the cluster with a 5 times higher frame rate. The last mode deals with only 1 pixel of the cluster to offer a 10 × higher frame rate. For those different modes biasing adjustments are provided to readout circuits.

VII. MONOLITHIC DEMONSTRATOR

The photography of the designed monolithic demonstrator is depicted in Fig. 5. It is composed of 12 clusters of 10 pixels with their 12 associated ADC. Some additional test structure are also included.

The memory block merged the 12 clusters data with a total of 12×16384 bits. Two clusters (pink) and their two associated ADCs (blue) are highlighted to show their shape factor and especially the area consumption. The configuration of the sensor and the extraction of the acquired data are passing through a SPI port.

VIII. EXPERIMENTAL RESULTS

The response of the sensor has been tested under a pulsed illumination. A 635 nm laser diode driven by a custom electronic emits a 15 pJ, 100 ps Full Width at Half Maximum (FWHM) long pulse of light that illuminates the whole sensor sensitive area [START_REF] Uhring | A low-cost high-repetition-rate picosecond laser diode pulse generator[END_REF]. The acquisition mode is "full frame" and the ADC dynamic is set to 8 bits with a 5 Mfps frame rate. The pulse of light is generated during the frame 4 exposure time. For clarity reason, only a single cluster is represented over the 7 first frames (Fig. 6). We clearly see that the light pulse is detected only on the frame 4. The pixel response non uniformity is meanly due to the speckles pattern of the laser diode. The pixel aperture against time is precisely measured thanks to the very short 100 ps FWHM light pulse swept across two frame periods. The result (Fig. 7) indicates that the aperture of the pixel is very sharp at the opening and decreases slowly at the end of the integration time. This is because when the laser pulses occurs at the end of the frame, the photodiode and the front-end do not have the time to setup completely before the sampling. When the pulse arises during the reset time of the frame i+1, a few crosstalk of about 1% is still present on the frame i, while the signal of frame i+1 slightly increases. This is due to an incomplete and not fully efficient reset of the CTIA. No influence is observed on the frame i+2 that benefits from a second reset period. A higher frame rate has been achieved by setting the master clock to 900 MHz, and only 2 pixels of the cluster are read by frame. The resulting frame rate is 900/10/2=45 Mfps. The measured pixel aperture is given in Fig. 8. Measured noise floor is 0.65 lsb rms. Our demonstrator specifications and readout characteristics are summarized in TABLE I.

IX. CONCLUSION

An 3D-IC architecture of a burst imaging CMOS sensor with a digital storage has been presented. A 2D monolithic demonstrator integrating the electronic of the future 3 different tiers has been designed, realized and characterized. The results show that the proposed architecture is fully functional and can operate at the frame rate of 5 Mfps. The used cluster approach allows to increase both the frame rate and the movie length at the cost of a spatial resolution reduction. The CTIA-based high speed pixel allows to reach a frame rate of 45 Mfps. This work paves the way to the very first 3D stacked BIS with digital storage.

Fig. 1 .

 1 Fig. 1. Cross section of the targeted 3D CMOS stacked

Fig. 5 .

 5 Fig. 5. Demonstrator die photography. Highlighted: two clusters of 10 pixel (pink) and their two associated ADCs (blue) and the shared memory bank for the block of 12 digital clusters.

Fig. 4 .Fig. 6 .

 46 Fig. 4. tier 1 (top) and tier 2 (bottom) floorplan partitionning (100µm×500µm double cluster)

Fig. 7 .Fig. 8 .

 78 Fig. 7. Aperture of a pixel at 5 Mfps and a reset time of 20 ns

TABLE I .

 I

		TECHNICAL SPECIFICATIONS
	Scalable Pixel Block specifications	pixel characteristic
	technology	7M1P 65nm	Pixel size	50 × 50 µm
	Pixel count	40	Fill factor	64%
	Area	0.1mm²	Readout characteristics
		5Mfps (full frame)	Read noise	2.5 mV rms
	Fame rate	45Mfps (2/10 pixel count)	Dyamic range	1 V
	Memory depth	204 full frame 8bit (1024 in 0.28nm) images	Conversion gain	5.9 µV/e-

ACKNOWLEDGMENT

We thank the Agence Nationale de la Recherche for financial support under grant ANR-14-CE26-0024-01