Some approximate Gauss – Newton-type methods for nonlinear ill-posed problems

This paper treats numerical methods for solving the nonlinear ill-posed equation F(x) = 0, where the operator F is a Fréchet differentiable operator from one Hilbert space into another Hilbert space. Two parametric approximate Gauss–Newtontype methods are developed, a local convergence theorem is proved under certain conditions on a test function and the required solution, and some computational aspects are discussed. The validity of the theoretical convergence rate estimates is illustrated by the numerical results of solving two sample problems, one in a finite-dimensional and the other in an infinite-dimensional Hilbert space.


INTRODUCTION
An important problem in scientific computing is that of solving a nonlinear equation where F is a differentiable operator from an abstract space into another.Further we shall suppose that F is Fréchet differentiable and it is acting from a Hilbert space H into another Hilbert space.Depending on the properties of the Fréchet derivative F , the problem under consideration may have a unique solution or infinitely many solutions or no solution in the classical sense and one has to seek for a generalized solution.
One way to determine a generalized solution is to seek it in the least squares sense by minimizing the functional Real-life problems are frequently ill-posed.For example, those arising in inverse problems are often illposed because they typically involve the estimation of certain quantities based on indirect measurements.Illposedness is an essential aspect of computation that is often seriously expressed in instability of the solution (the solution is unstable under data perturbations) and the usage of standard, non-regularized methods can produce uncontrollable and unacceptable error propagation.For solving ill-posed problems regularizing algorithms containing several variable parameters may be fruitful, because a proper selection of parameters involved sometimes improves the convergence properties, reduces the amount of computation, and provides a wider choice of initial guesses.Much work has been devoted to the study of the convergence of the iteratively regularized methods.Not long ago, many monographs, e.g.[1,7], and a number of papers [3,6,11,12,16,17], were devoted to this topic.In the works mentioned above, as a rule, problems with the inaccurate data are considered and the convergence of the Gauss-Newton method is proved under the assumption on source-like representation of the exact solution x * : x holding for an element v ∈ H.The writing F * (x) denotes the adjoint operator of F (x).In order to obtain simple estimates of the terms including approximations to the inverse operators, here we use the value p = 1.
It is well known that without extra assumptions, usually source-type assumptions or some assumptions concerning the smoothness of the solution, one cannot get a specific rate of convergence.In contrast to many other papers on ill-posed problems, we do not consider the question of choosing the level of regularization based on the known or estimated level of the noise in the data.In this paper the data are assumed to be exact and we consider the question of finding a solution of (1) in the least squares sense when F * (x)F (x) is not necessarily continuously invertible for all x in a sufficiently large neighbourhood of the exact solution x * .This paper is a continuation of papers [8,10,13,14] and it treats approximate Gauss-Newton-type methods for solving nonlinear ill-posed problems for which procedural and rounding errors are unavoidable.Frequently the use of finite-difference approximations to the derivatives gives rise to an inexact method.Approximations may also be regarded as a result of inevitable inaccuracy of computation or they are due to the strategy used for solving linear problems at each iteration, i.e. associated linear equations are solved intentionally approximately by taking finitely many steps of an iterative procedure or the inverse operator appearing in a method is approximated by a formula that is less time-consuming to compute.The purpose of this paper is to discuss basic concepts of the approach and therefore the emphasis is on the theoretical background rather than the practical solution.

METHODS AND CONVERGENCE THEOREM
A large number of real-life problems are solved by finding the solution of certain equation or by minimizing a functional under consideration.Here for finding stationary points of the functional (2) or, equivalently, for solving the nonlinear normal equation F * (x)F(x) = 0 a family of iterative methods based on the notion of regularization is developed.Recall that here F * (x) denotes the adjoint operator of F (x).
We shall consider the case where the range of F (x) is not necessarily closed and the product operator F * (x)F (x) may not have the bounded inverse operator for every x.The Gauss-Newton method 2 at every iteration point x k for the correction term h.In order to cope with the possibility of having a noninvertible or ill-conditioned F (x k ) at some iteration step it is desirable to develop algorithms based on the Tikhonov functional where α > 0 and u 0 is an element of H, the so-called test function.
Following the idea of iterative regularization [2,5,15,16], we assume that α in the Tikhonov functional is changed at each iteration step according to a sequence {α k } of properly chosen positive numbers.From the stationarity condition for the Tikhonov functional we obtain a one-parameter iterative Gauss-Newton method where M k = B(x k ) + α k I and I denotes the identity mapping.As B(x) is positive semidefinite for every x ∈ H, the operator M k is clearly invertible for every k.
In this paper we shall study a class of two-parameter regularized Gauss-Newton-type methods where 0 < ε k ≤ 2 and D k is an approximation to M −1 k satisfying the condition The motivation for introducing the approximate inverse D k of M k is the fact that often it is not possible or reasonable to apply M −1 k to an element x ∈ H exactly and the result is found by solving the equation M k y = x approximately by a numerical method.
To get convergence results for the method (3), let us suppose that the operator F is Fréchet differentiable with Lipschitz continuous derivative in the region S ⊂ H under consideration with and that equation (1) has a solution x * with the representation holding for an element v ∈ H and a test function u 0 .We shall show that for an appropriate choice of α k the sequence {x k } defined by (3) converges to the solution x * of (1) provided the parameters ε k satisfy the condition and that the parameters µ k and v are sufficiently small, i.e. and In order to make the statement and the proof of the convergence result shorter, let us introduce some additional constants: Now we are ready to state the convergence theorem.
Theorem 1.Let F be a Fréchet differentiable operator with Lipschitz continuous derivative and the relations (4)-( 8) hold.If the regularization parameter α k is chosen according to where the coefficient c β satisfies the inequality the sequence {x k } defined by (3) converges to the solution x * with If the test function u 0 is not equal to x * , parameter δ is well-defined since for β ≥ q we have x * , we can choose the parameters α k according to the conditions of the theorem only for β > q.
Proof.Denote h k = x k − x * and τ k = x k − x * .Note that according to (13), we have τ 0 ≤ δ α 0 .Assume that τ k ≤ δ α k for some k ≥ 0 and show that then also Since, according to assumption (4), we have we can write where Taking into account that F(x * ) = 0 and the relation (3), we have On account of the relation ( 5) we obtain Therefore we can estimate Note that from (4) it follows that the inequality is valid.Making use of ( 15) and the inequalities Using the assumption τ k ≤ δ α k and the notations ( 9)-( 11), we get Since δ is a solution of the equation Thus, according to the principle of mathematical induction, we have proved that The statement of the theorem is proved.

NUMERICAL EXPERIMENTS
In order to illustrate the accuracy of our theoretical convergence result, we consider two test problems.
Our first test problem is solving the algebraic system of equations where the function F is the Powell singular function (see [9]) The exact solution is x * = (0, 0, 0, 0) T and the Jacobian F , given by is singular at x * .Note that for such problems also unmodified Newton's method has a linear convergence rate (see [4]).For the initial approximation we took the assumption (5) requires us to choose u 0 in the form u 0 = (a, 10a, b, −b) for some real numbers a and b.As for this problem the coefficients K and L in (4) are not small (by direct computations it is possible to show that K can be taken to be 13 in the unit sphere and L is √ 160), the norm of the vector v in (5) should be quite small for the local convergence result to apply.We used a = 0.001 and b = 0.001 in our computations; this corresponds to v ≈ 1.7 × 10 −4 .
First we present the results obtained for method (3) with the exact inverse D k = M −1 k .We considered three different constant step sizes ε k (namely 0.7, 1.0, and 1.3) and also the method with ε k chosen randomly in the interval [0.7, 1.3].According to Remark 1, we can use c β = 690 for all values of β .Thus the regularization parameter was chosen in the form α k = 690 β k for β ∈ {0.9, 0.7, 0.5, 0.3}.The results corresponding to 30 iterations are presented in Table 1.Here τ k is the Euclidean distance between x * and x k .We can see the observed reductions in errors during the last four iterations, the average observed reduction coefficient based on the last 10 iterations, and the final error.The last line of the table corresponds to the standard Gauss-Newton method (ε k = 1, β = 0).According to our result, the theoretical overall reduction rate of the error is guaranteed to be at least β for β ≥ 0.684 if ε k = ε ∈ {0.7, 1.3} and for β ≥ 0.295 if ε k = ε = 1.We see that the numerical results confirm the theoretical estimates.It is not surprising that for the concrete problem we may get a higher convergence rate by choosing a value of β that is smaller than the ones allowed by the theorem but the last line for the standard Gauss-Newton method shows clearly that a too small value may cause slower convergence speeds.
In order to test our theoretical results in the case of approximate inversion of M k , we used the method with D 0 = M −1 0 and The computations were done for β ∈ {0.9, 0.8, 0.7, 0.6, 0.5} and it was verified that the accuracy of approximation of M −1 k was not smaller than µ 0 = 0.1.The corresponding results are presented in Table 2. From the convergence theorem it follows that in the case µ 0 = 0.1 the linear convergence rate β is guaranteed for β ≥ 0.43 if ε k = 1 and for β ≥ 0.86 if ε k = ε ∈ {0.7, 1.3}.We see that the statement of the theorem is confirmed.
For our second test problem we consider the equation with the exact solution We view it as an operator equation in the space L 2 (0, 0.5) of the square integrable functions.By denoting the former equation can be viewed as an operator equation of the form Clearly the Fréchet derivative is and its adjoint operator is By direct calculation we get that the composition F * (x)F (x) is the integral operator where the function κ x depends on x and is defined by Note that for any x, z ∈ L 2 (0, 0.5) the function [F * (x)F (x)z](t) is 0 at t = 0.5 and its (one-sided) derivative is 0 at t = 0, so the condition for some v ∈ L 2 (0, 0.5) requires us to choose u 0 so that it is continuously differentiable and satisfies u 0 (0.5) = 2 √ 3 , u 0 (0) = 0.It can be shown that any function u 0 that satisfies the given conditions and has a second derivative in L 2 (0, 0.5) is an acceptable choice.We use the function for our test function.The corresponding function v is and by direct computation we get v ≈ 0.49.It is relatively easy to check that the conditions (5) are satisfied with K = 1 √ 8 , L = 1 4 , so from Remark 1 we get that any coefficient c β ≥ 0.23 can be used for defining α k in Theorem 1.We are going to use c β = 1.
We start our iteration process by choosing x 0 (t) = 0, t ∈ [0, 0.5].In order to compute approximately we solve the integral equation numerically as follows.We fix the number of intervals n = n(k) and consider the uniform grid and determine a piecewise constant approximate solution with the values ξ i on the subintervals (t i−1 ,t i ], i = 1, . . ., n, to the integral equation by requiring with This procedure corresponds to using a fully discretized collocation method and can be viewed as using an approximation D k instead of M −1 k .In order to keep the approximation error under control as α k gets smaller, we have to increase the dimension n with k.We used n(1) = 2 and doubled the value of n at even values of k (so n(2) = 4, n(3) = 4, n(4) = 8 and so on).We used the method (3) with a constant step size ε k = ε for three different values of ε.The regularization parameter α k was chosen to be of the form α k = β k for five different values of β .The numerical results are presented in Table 3.
Here τ k = x k −x * L 2 .In total, 25 iteration steps were carried through.The observed reductions in errors during the last four iterations, the average observed reduction coefficient based on the last 10 iterations and the final error are presented in the table.According to Theorem 1, the theoretical asymptotic reduction rate of the error is guaranteed to be at least β for sufficiently large β < 1 and the numerical results confirm the validity of the estimate: for β ∈ {0.9, 0.8, 0.7} the observed reduction rate of the error is approaching the theoretical value.Our theorem does not cover the behaviour of the method if the parameter β is too small (less than q defined in (8), which is unfortunately unknown since the coefficient µ 0 is not known).From the numerical results we see that in the case of the test problem the convergence takes place also for relatively small values of β (namely for β ∈ {0.4,0.2}) but the observed reduction rate is not related to the value of β any more and that a too small value of β does not improve the accuracy of the final result.

CONCLUDING REMARKS
In this paper we studied a family of two parametric iterative regularization methods based on the Gauss-Newton method for solving ill-posed problems in Hilbert spaces.It turned out that under suitable assumptions about the test function the family of methods has a geometric rate of convergence near the exact solution.The two numerical examples (one in a finite-, the other in an infinite-dimensional Hilbert space) confirm the validity of the obtained results.Of course, from the practical point of view there still remains the problem of choosing a good test function.Likewise, finding a proper regularization parameter α and, especially, a suitable compromise between the values of the parameters α and ε is a fragile problem.It is well known that mathematically efficient methods do not necessarily result in efficient computer programs.On the other hand, the worst case bounds of a method may be too pessimistic and in practice it may perform better than expected.Therefore the determination of the parameters on the basis of the heuristic considerations certainly deserves further investigation as well.

Table 1 .
Numerical results for the exact D k

Table 2 .
Numerical results for an approximate D k