For system of linear algebraic equations
with given approximate (m*n)-matrices A (m>n) and m-dimensional vector u, the solution can often be determined from the equation
(1)where S is an (n*n)-matrix given a priori and a is the regularization parameter. The expected "smoothness" of the solution, together with its sensitivity to errors in the initial data, is known [1] to serve as the basis for using (1).
If formal "smoothness" criterion is defined by the condition
or
(where j is the index of component of za), then the multiplication of za by scalar beta can violate both conditions, although one might think intuitively that smoothness of beta*za is adequate for za and beta plays the role of the scale factor. We shall assume that if, in view of any criterion, the solution beta*za1 is preferable to za2, we can choose the former without violation the assumptions made or a priori concepts.
The parameters a1 and a2 can be distinct. This assertion rests on intuition and serves as an additional criterion for choosing the approximate solution along with other known factors.
Suppose that, in order to choose a , system (1) is supplemented by the condition
(2)where the square of the norm is defined as the scalar product
and delta is a character of the error. Let za1 be a solution of (1) that satisfies (2). Then, if a scalar beta exists such that
then on fixing beta, we can increase a to choose a smoother solution.
Consequently, the numerical problem is the following. For the given system (1), it is required to find the greatest possible parameter a and arbitrary parameter beta such that
(3)If beta is fixed, then F(a,beta) increases monotonically [2] and has a horizontal asymptotic ||u||2. Thus , if the maximum with respect to a is required, then the minimum with respect to beta must obviously serve as the second condition.
To develop a rational algorithm for solving the problem, we shall consider the properties of F(a,beta).
For any fixed a that ensures the solvability of (1), there exists a unique parameter beta0 that can easily be determined from the condition
and is given by
(4)from which it can be seen that beta0 a > or = 1. Taking (1) into account, the right-hand side of (4) can be obtained by elementary transformations of the middle part. The sign of the second derivative d2F(a,beta)/dbeta20 determins beta0 as a local minimum. The fact that, for any a > 0 such that (1) is solvable, there exists a unique parameter beta0 a that minimizes F(a,beta), makes the problem much simpler and enables one to replace (3) by the condition
(5)We shall investigate the function
(6)To simplify the notation, we will omit the superscript a attached to z and beta0, assuming that the components of z are dependent of a, as well as the scalar beta0.
We will differentiate phi(a) with respect to a
One can see from (4) that the factor multiplying beta0' is equal to zero, and so
To determine the properties of phi', we will use the well-known method [3,4] of singular decomposition of the matrix
(7)where P and Q are orthogonal matrices and L is the diagonal matrix of singular numbers. The decomposition (7) will be used exactly in the proof. Having transformed (1) to the form
one can easily show using (7) that
and
where g=PTu. Next, transforming (Az, Az')=(AS-1Sz, AS-1Sz') and (Az', u)=(AS-1Sz',u), we get
where lambdai are the diagonal elements of L, and gj are the components of g.
To determine the sign of the derivative (4), it is sufficient to consider the expression
which, taking into account that
we can transform to the form
(8)
Consider the right-hand side of (8). It is obvious that, on multiplying the items under the summation sign in the appropriate order, one can use grouping to distinguish (n2-n) expressions in the form
(9)
for j=1,2,...,n, k=1,2,...,n and j not equal k . The remaining n expressions for j=k will have another form and their sum is equal to zero. If (9) has a definite sign, then both psi(a) and phi', which is of interest to us, share the same sign. Note that if one of the numbers lambdaj2, lambdak2, gj2, gk2 is equal to zero, the whole expression (9) vanishes for a>0, and so we shall not consider this case. We divide (9) by lambdaj2, lambdak2, gj2, gk2 and multiply by (lambdaj2+a)2(lambdak2+a)2 to get
Next, by means of simple transformations, we can reduce the latter to the form
Consequently psi(a) > or = 0, phi'(a) > or = 0 and phi(a) is non-decreasing function. It follows that if a solution of (1) exists that satisfies (5), then it is unique. The case phi'(a)=0 occurs when all the singular numbers of AS-1 are equal to each other. (AS-1)T(AS-1) is then a diagonal matrix, and such a system of equations does not satisfy the definition of an ill-posed problem. Nevertheless, we must turn our attention to this ideal case in order to better understand role of beta. If L=I in (7), then beta0=1+a and, in addition, the solution beta0 za does not depend on a any longer and is given by beta0 z = (STS)-1ATu, from which it is obvious that the discrepancy (6) is also constant over the whole interval of possible variation of a. Therefore, beta0 eliminates the side effect of applying the smoothing functional, i.e. the reduction of the (absolute values of) components of za. The equation phi(a)=delta may have no solution unless delta belongs to <[phi(0); phi(infinity)], in which case the estimation of delta is incorrect.
The direct algorithm for solving (1) and (5) does not involve any difficulties an can be based, for example, on method of dividing in half, for convergence of which it is suffices that the localization interval of the root of Eq. (5) is specified.
Finally we shall present the data of numerical experiment. The matrices AT and B and the vectors zT and E were given as the input data, so that ATzT=uT constituted the exact system and Az=u or (AT+B)z=uT+E, the approximate system. The parameter delta was estimated by the formula delta=ETE. The given exact system has the unique solution zT=(ATA)-1ATTuT with zero discrepancy in the sense of the best approximation. Th main goal of the experiment was to compare za obtained from (1) by means of the discrepancy condition
(10)
with the solution beta0 z a obtained from the condition
(11)To distinguish between (10) and (11) in the text, we call the former the a-discrepancy and the latter beta-discrepancy, even though (11) depends formally on a. In Fig.1 we present zT and za1 obtained from (10) and beta za2 from (11). The a- and beta-discrepancy functions are both shown in Fig.2 along with the magnitude of delta. The two solutions za1 and beta za2 have norms equal to zT, but the purpose of introducing the beta-discrepancy has been achieved, as beta za2 is better agreement with a priori smoothness requirements then za1. We observe that a1=504 is nearly four times smaller than a2=1908, which naturally ensures better "smoothing" in the latter case. The unit matrix was taken as S in (1). The elements in the main diagonal of ATA lay between 921 and 1054, therefore a1 exceeded them by a factor of two. The computations were carried out on an SM 1800 computer using numerical variables with eight significant digits. It was observed that beta0 za is more stable with respect to errors in the data, and, in particular, it turned out that it was much less sensitive to changes B and E. The trend in Fig.1 was present in various experiments we conducted.
REFERENCES
1. ARSENIN V.Ya. and TIKHONOV A.N., Methods of solving Ill-posed Problems. Nauka, Moscow, 1986.
2. MOROZOV V.A., Regular Methods for Solving Ill-posed Problems. Nauka, Moscow, 1987.
3. GORDONOVA V.I. and MOROZOV V.V., Numerical algorithm for choosing the parameter in the regularization method. Zh. Vychisl. Mat. Mat. Fiz. 13, 539-545, 1973.
4. POLUEKTOV A.R. The method of successive approximation for choosing the regularization paramter for the numerical solution of ill-posed problems. Zh. Vychisl. Mat. Mat. Fiz. 29, 1734-1737, 1989.