Andrew Polar (Andrei Poluektov)

Published in Computational Mathematics and Mathematical Physics Vol. 29, No. 6, pp. 95-97, 1989 (Great Britain).

A METHOD OF SUCCESSIVE APPROXIMATIONS FOR SELECTING THE REGULARIZATION PARAMETER IN THE NUMERICAL SOLUTION OF ILL-POSED PROBLEMS

An algorithm is proposed for the numerical solution of ill-posed problems, which can be reduced to the system of linear algebraic equations with inaccurately specified operator and vector. The algorithm is based on the method of successive approximations. It is proved that the algorithm converges to the required solution.

The numerical solution of great many ill-posed problems by regularization reduces to solving a system of linear algebraic equations with inaccurately specified matrix A (of order n*m) and vector u (of dimension n), of the form

(1)

subject to an additional condition for selecting the regularization parameter a /1/:

(2)

where I is the identity matrix of the same order as ATA, za is an m-dimentional solution vector of (1) for the given parameter values, h>0, delta>0, mu>=0 are constants characterizing the accuracy of the initial data and the degree of incompatibility of the initial equation Az=u, and the square of norm in (2) is defined as the scalar product of itself of a vector of the appropriate dimension (from now on indices m and n, indicating dimentionality, will be omitted).
Many well-known algorithms /2, pp.219, 220; 3/ for solving (1), (2) are based on the properties of the auxiliary functions ro(a)=||Aza-u|| and gamma(a)=||za||, which are corollaries of the following analytical expressions, established by using the singular decomposition A=P1LP2 (see /4/):

(3)

where lambdai2 is the i-th singular value of the matrix ATA, gi is the i-th component of the vector P1Tu and P1 and P2 orthogonal operators. The main disadvantage of these methods is the need to compute derivatives roa' and gammaa'. In this paper we consider an algorithm in which both these and many other procedures involving the computation of derivatives are not needed.
Consider the function



Among its obvious properties are the following. Since gamma(a) is a decreasing function and ro(a) is increasing function over the entire interval of increasing a, and moreover the constants h,delta,mu are such that the functions intersect at a single point (denote it by as) we can write

(4)

if we assume that phi(a) increases monotonically, it is easy to see that phi(a) is a contractive operator and equation a=phi(a) can be solved by successive approximations. to prove this assertion, consider an arbitrary alfa0 smaller than alfaS. Since phi(a) is increasing, it follows that



and by (4) we have phi(a0)=a1 greater than a0, so a0 smaller than a1 , smaller than aS. Similarly, we see that a1 smaller than a2 , smaller than aS and the sequence aj=phi(aj-1), j=1,2,... converges from below to aS; if the initial approximation has been greater than solution, the convergence would have been from above.
To prove that phi(a) is indeed an increasing function, we need only prove that the following auxiliary function psi(a):



where



(the notation lambdai and gi is as in (3)) is decreasing, since



and 2h+(gamma(a))-1(2delta+mu) is a monotone increasing function.
since by assumption psi' non-negative,

(5)



Obviously, a necessary condition for (5) to be true is that



Multiplying and dividing the fractions with denominator (lambdai2+a)2 on either side of the inequality by lambdai2+a and then subtracting



from both sides, we obtain



which is just the Cauchy inequality. This is easy proved by substituting



Thus (5) is true and phi(a) is indeed monotone increasing function.
Remark. Since the equation a=phi(a) has an irrelevant root a=0, the initial approximation should be chosen to be positive, even if (1) is solvable for a=0.
Our algorithm is useful for a more general form of (1):



where C is preassigned symmetric, positive-definite matrix. In this case one need only apply the substitution used in /5/.
To illustrate the rate of convergence, we present data from a numerical solution of the Wiener-Hopf equation. The operator A and vector u were constructed from experimental data. The investigation concerned the operation of the diesel generator as a linear dynamic object. We present only data pertaining to the algorithm proposed above. Tables 1 and 2 demonstrate two computations of the regularization parameter by successive approximation. The initial approximations in the examples were selected to be one order of magnitude less and two orders of magnitude greater than the unknown value. Apart from the already mentioned auxiliary function ro(a) and gamma(a), the table also lists values of the function



which approaches unity from above or below.
Table 1
# alfa ro gamma ksy alfa
1 0.0010 0.0556 3.9118 4.9299 0.0049
2 0.0049 0.1345 3.0477 1.6715 0.0082
3 0.0082 0.1717 2.7385 1.2060 0.0993
4 0.0993 0.1862 2.6319 1.0799 0.0107
5 0.0107 0.1923 2.5896 1.0335 0.0110
6 0.0110 0.1949 2.5717 1.0145 0.0112
7 0.0112 0.1962 2.5639 1.0064 0.0113
7 0.0113 0.1967 2.5605 1.0027 0.0113


Table 2
# alfa ro gamma ksy alfa
1 1.0000 0.8940 1.0330 0.1232 0.1232
2 0.1232 0.4348 1.6416 0.3331 0.0410
3 0.0410 0.3054 2.0021 0.5414 0.0222
4 0.0222 0.2516 2.2368 0.7107 0.0158
5 0.0158 0.2233 2.3910 0.8396 0.0132
6 0.0132 0.2092 2.4776 0.9198 0.0122
7 0.0112 0.2024 2.5209 0.9633 0.0117
8 0.0117 0.1994 2.5411 0.9838 0.0115
9 0.0115 0.1983 2.5503 0.9923 0.0114

REFERENCES
1. TIKHONOV A.N. and ARSENIN V.YA., Methods of Solving Ill-Posed Problems, Nauka, Moscow, 1986.
2. MOROZOV V.A., Regular methods for Solving Ill-Posed Problems, Nauka, Moscow, 1987.
3. LEVIN A.M., On a modification of Newton's method and the method of secants and their application in regularizing algorithms. Zh. Vychisl. Mat. Mat. Fiz., 28, 8, 1123-1134, 1988.
4. GORDONOVA V.I. and MOROZOV V.V., On the regularization technique. Zh. Vychisl. Mat. Mat. Fiz., 9,3, 673-675, 1969.