The error principle in the solution of incompatible equations by Tikhonov regularization

THE ERROR PRINCIPLE IN THE SOLUTION OF INCOMPATIBLE EQUATIONS BY TIKHONOV REGULARIZATION* V. A. MOROZOV Moscow (Received 29 April 1972; reuised 10 January 1973) THE construction of approximate solutions on the basis of a principle whereby the regularization parameter is selected according to the size of the error is discussed. While the treatment is in the context of the finite-dimensional case, all the propositions are in such a form that they hold for any linear operators acting in Hilbert spaces. We examine the possibility of constructing approximate solutions when incompatible equations are solved by Tikhonov regularization [ 1 ] , the construction being based on selection of the parameter according to the size of the error. We outlined in [24] a similar method for finding approximate solutions in the case of compatible operator equations. The present regularization scheme is quite general and includes the case when the stabilizing functional does not define the norm. The regularization method was discussed under such conditions for the first time in [5]. While our treatment refers to the case when all the spaces are finite-dimensional, our propositions are given in such a form that they also hold for more general cases, notably, when we are concerned with linear closed operators acting in Hilbert spaces. The present paper is based on the author’s report to the All-Union Conference on Computational Methods of Linear Algebra, held at Novosibirsk, in March 1970. 1. Let A :R, + R, and L :R, -+R,be matrix operators, acting in Euclidean spaces R,,R, and R, of n, m, and s dimensions respectively, and let II- I(,,, (I - II,,,, II - II8 be the vector norms in these spaces. We assume that the quadratic form

is positive definite, i.e., there exists a number 7 >

0 such that, for any u E R,

Ilull*, .27lldln. *Zh.vpchisl..Mat.

mat. Fiz., 13, 5, 1099-1111, 1973

1

U.S.S.R.

Vol.

13. No.

5-A

(1)

2

V.A. Morozov

Notice that (1) only demands that the over-all quadratic form defined by the matrix operators A and L be positive definite, while allowing that each individual form may be degenerate. Consider the problem of finding a vector uOE R,, satisfying the conditions mA

=

inf IlAu - film =IIAuO - film, UER*

m,’ = inf II&z - gll, =IILu, - gll,,

(2)

U&VA NA

‘{U

=

R,: mA =llAU - fs,),

where f E R,, g E R, are vectors. It can easily be seen that problem (2) includes the case discussed in [5]. The statement (2) of the problem is clearly desirable in a number of cases of importance in the theory of the method of least squares. For instance, we arrive at problem (2) when attempting to find a solution of the conditional system of equations Au=f,

(3)

for which the linear form Lu deviates least from a previously assigned vector g. In particular, if the vector g = Lu*, where U* is a ftied vector, obtained e.g., by direct measurements of the required vector U, problem (2) amounts to finding a vector ue which, in some specific sense, best corresponds both to “indirect” and to “direct” measurements. We can also treat the system (3) as a system of preassigned linear connections satisfied by a required vector-parameter (and originating e.g., from the desire to satisfy given laws of conservation [6]). The solution of problem (2) is then a vector u. which best corresponds to a system of conditional equations Lu = g

and satisfies the linear relationships (3). Following [ 1] , we introduce the parametric functional (a > 0)

Rt[ul = IMu -

fllmZ+ allLu -

glls

and consider the variational problem

%[ul,

UER,

-

min.

(4)

Error principle in the solution of incompatible equations

It can be shown [5] that, under condition (l), the solution u. of problem (2) and the solution uo of problem (4) exist and are uniquely determined, given any f and g and any parameter a > 0 by the system of linear equations aL’(Lu,

- g) + A’(An, - f) == 0,

where the subscript T denotes transposition.

q”(a)

2. Wesetp(a) = Mu, - ftl,, = @p,[ua].

y(a) = IILu, - gll,,

$(a) = p2(a-f+

We shall establish below the limits of these functions, and also the behaviour of the reduction family of vectors & as a + 0 + and a -+ ~0. By analogy with Eqs. (2), we set NL ={u E R,: m, =Ih

mL = inf IlLu - gll,,

- gll,),

UER,

m,’

=

inf IlAu - film. UENL

It can easily be seen that m, < mA’, ntL G mL';thevector u,, mL =

IlLh -

satisfying the conditions

mA’= Mum - flL,

gll,,

is then uniquely determined, like the vector u. . Lemma I

We have

lim II@,- aoIIA,L=

a-o*

0,

lim p(a) = lim rp” (a) = mAt C&+0-+

(5)

S-+0+-

lim 7 (a) = m,'. CS+O+ Proof: Using the extremal properties of the vectors ug and ua, we get the chain of inequalities mA2G p"(a)G q(a) =p"(a) + ay"(a)G &[uo] = ma2 + c6mt'* G p'(a)4 amL'2, from which it follows that

mA

G

p(a),

p(a) G

mA

+

Ia

m,',y(a) G mL'.Hence,

4

V. A. Morozov

lim p(a) = m,, CA-O+

(6)

limy (a)

(7)

and also

C&+0+

< m=‘.

We will show that the first of Eqs. (5) follows from Eqs. (6) and (7). In fact, from Eqs. (l), (6), and (7), the family ua is bounded. Let cz, be an arbitrary sequence of positive numbers, having limit zero as 7~+ 00. Since the sequence u, = ua,, is also bounded, we can isolate from it a convergent subsequence. For simplicity, this will be assumed to be the same as the sequence u,. We set u^= lim u,. It then follows from Eq. (6) that u^E Na, and from Eq. (7) that IlLU - .&YrnL’. Hence i is a solution of problem (2), and since it is unique, u = z&J. Since the sequence a,, + 0,as n -+ 00, is arbitrary, and using what has been proved, we find that the first of (5) holds; the remaining relationships of (5) are a simple consequence of the first. The lemma is proved. Let us examine the behaviour of these functions as a * 00. In fact, Lemma 2

lim p (a) = mA’,

lim llua - rzmllA,I.= 0,

DLtm

a-m

lim y (a) = rnL,

a-too

lim [cp(a) - mA” - am=‘] = 0. a-c-m

Proo_f We have

p”(a) + af(a)

= @a[&1 G aL[uml

= m, ‘2 + ccmL2 < mA” + ayZ (a). Hence limp (a) < mA’, a+-

Since mL G llLua - gll = r(a),

we have

limy (a) G mL. a-bm

03)

Error principle in the solution of incompatible equations

In conjunction with the earlier inequalities, this gives lim T (a) = mL, a-m

lim p(a) G m,‘. a-PC.2

On using the same arguments as when proving Lemma 1, we obtain from these last expressions lim llua - UmliA,L= 0, a+m i.e., the first of Eqs. (8) has been proved. The other relationships of Eqs. (8) follow from the first. The lemma is proved. Corollary 1

Let &I= LU*;then, obviously, m, = 0 and the relationships (8) become lim p (a) = lim cp”(a) = mA’, a+m a*m

lim y(a) = 0. .z-rm

(8’)

If the rank of L is n, we have ma’ = IIAu’ - fll,,,. The relationships (8’) also hold when the vector g = 0 and L has arbitrary rank; then, na’ = Ilf]],. Corollary 2

It follows from Lemmas 1 and 2 that the values of the function p (a), 0 C a < 00, exhaust the interval (ma, M*‘) where ma =Cnn’. Lemma 3

If m, = m,‘, then u, = u. = u,. Pro05 Using the extremal properties of the vectors ua, uo and u,,we have

ax[u,l

G %[zL]

= ma2 -t

= rn*‘?.+ ccrnL2 =zm*‘z + allLu, - gll*Z

dlLu, - gll,2 G IIAu, - fllm2+ allh,

- gll.” = @,,[~a].

Hence CD,[ u,] = a,,[ u, 1. Since the solution of problem (4) is unique, u, = u,. Since lim ua = uo, we have u. = u-. The lemma is proved. a-Al+ The following lemma is proved similarly.

6

K A. Morozov

Lemma 4 If mL = m,‘, then ua = u. = u,. The conditions of Lemmas 3 and 4 are satisfied, e.g., if a vector I$ exists such that g = Cu’and Au’ = j(notice that U* is uniquely defined by these conditions). It is easily seen that then, ma = m,’ = m, = m,’ = 0, while the regularization vectors ua = U’. This case refers to the situation when the mathemati~l model Au = j is ideally matched to the results of direct me~urements of the required vector. 3. Next take the case when the matrix A and vector j”are specified approximately, i.e., instead of A and f we know a matrix Ah and vector j5 such that II(A. - AJ&I

< ~llr&4, z.,

Ilf - AlI c 6.

The quantities h and 6 are assumed to be sufficiently small. We set iiin = inf lldhu - fall, UERn It is easily seen that 6ZAd fld& limiiiA 0-d

jji~’ = inf IJAhU- jJ,

=11&Z, - j811m.

UENL

- j&, < I(

m&

mA

+

8 +

hlluollA,

L.

COnSeqUently,

o=(h,6).

We then find, in the same way as above, that (9) i.e., lim.riia’ < m,‘. o-to

(10)

In addition, it follows from relation (9) that IL&i, ll,G’,, jl~~~~l~~~lgl~~ and hence A‘L G Cz, where the C+ are ~dep~~dent of o (we assume that o is suf~ciently IIK%?li small). Then, since

we have

From relations (10) and (lo’), we have lim iiia’ = mA’. O-+0

Error principle in the solution of incompatible equations

We have thus proved Lemma 5

Given any o we have ?Z, < \\A& - fa\lm and in addition,

limEA<

mar

We set p(a) problem

Iim ?ZiA’= ma’. O-*0

O-CO

= llAhEa - fall,,,, where & is the solution of the variational

6&l = II&u

u = R, - min.

- fAl,n2+ all-h - dl,“,

(11)

Let us now assume that ma < m.*‘.

(12)

It follows from this inequality and Lemma 5 that quantities pU,can be found, satisfying, for sufficiently small (3 the conditions lim p. =

mn,

llAhuo-faLsp,~

64'.

(13)

a-0

Lemma 6

Let conditions (12) and (13) be satisfied. The equation P (4

= pu

(14)

then has a unique solution a = 01(o) > 0 for sufficiently small o. For, the values of p (a) 0 < a < 00 exhaust the interval (%A) f&l). Since, obviously, pUE (tZ_+ iii*‘), Eq. (14) is solvable. To show that Eq. (14) is uniquely solvable, we only need to show that the function 0 (a) is monotonically strictly increasing. In fact, for the vectors & we have Euler’s equation c&(X,

- g) + Ah=(A,,fi, - fa) = 0,

while the vector derivative di.?, / da is given by the equation

(aL'L+

A.‘A,)!$ =- L’(L&

- g).

V. A. Morozov

8

By definition, p”(a)

= II&u”, - f~]],,,‘,whence dp(a)

P(a)da=

du”,

(da’

AT (A

&-A)) . n

Using the equations written above for 6, and dZ, / da, we get

4 (4

p(a)-

da

= a ( (aL’L

+

Ah’Ah) -’

XL’(Lu”,-g),L’(L~,-g)),~O. Let us show that the strict inequality in fact holds for all a > 0.Assume the contrary, and let a > 0 be such that do (a) / da=O. Then, LT(L& - g) = 0 and hence Ahr (A,,& - f6) = 0, i.e., the vector G, belongs simultaneously to the set NL and the set NAh. Hence ]]A,& - fall,,, = iiia < %,‘< /IA&, - fallm,i.e., iii* = EA’ for all sufficiently small (3. But this is impossible, since, by Lemma 4, the assumption (12) implies the inequality Si, < Zia’. This proves the strict inequality and hence the lemma. Note. It may be shown in the same way that the function $(a) is also strictly increasing, while T(a) is strictly decreasing, under condition (12). These statements were obtained by another method in [4] for a less general case. 4. Let us now prove our main proposition. Theorem 1

Let a = a(o) be defined in accordance with the error principle (14). Then, lim IIU”,- u~II~,~= lim II& - &In = 0. o-to ,a-+0

(15)

praof: Using the extremal properties of the vector &, we get s,[ &] G za[ uO] < poz + amL’2 and hence p,,’ -I- a]]& - g]].” ( pm2+ amL’2. Hence limIIL& - gll, d mL’. 0-M

Further, using the triangle inequality, we have lM& -- fll m G pu + 6 + hll &lla, I,.

(16)

9

Error principle in the solution of incompatible equations

It follows from relation (16) that ]IZ,llA, 2 < con&. From this and Eq. (13),

Since, on the other hand, mA G llAiia

-

fll,,

we

have

1imllA u”, - film = m,. o-to

(17)

Using arguments similar to those in the proof of Lemma 1, we find that Eq. (15) follows from relations (16) and (17). The theorem is proved. 5. Again let conditions

(12) and (13) be satisfied. We define a set D, =

ilAhU - fbllmG pa}

{U = R,,:

and consider the following extremal problem:

to find

a vector u. E R, such that inf IlLu - gil, =IILU, USDli

- gll..

It is easily shown that problem (18) has at least one solution.

Theorem 2 If a = a (o) is defined in accordance with condition

(14), the vector

C., is a

solution of problem (18). u”, E D, obviously holds. On the other hand, for every vector

For, the inclusion

u=R,we

have from(l1):

~,[u”,]&r,[u].

On settinghere

is any solution of problem (18)) we also find that JILU”, vector

a = a(o), g]ls < ]]Lu,, -

u = u,, (u, gIla, i.e., the

u”, is in fact a solution of (18). The theorem is proved. It helps to clarify the meaning of the solution obtained by

the regularization method (11) when the choice of the regularization on the so-called error principle (14).

parameter

a is based

The following is also of interest.

Theorem 3 Let U, be the set of solutions of problem (18). Then, lim sup {max(iiu o-+0 u=Lln Apart from slight modifications,

- UoltA,L,11u- &ll,,)) = 0.

the proof is the same as the proof of Theorem

1.

10

V. A. Morozov

There was a previous discussion in [3-8J of an error method similar to Eq. (18), for the case of compatible equations. 6. Let us next consider some aspects of the numerical realization of the above ~go~t~s. It is easy to see that the important thing is to be able to determine effectively the root of equation (14). Unfortunately, the function p (a) can change its curvature, so that it is extremely difficult to apply Newton’s method to the solution of Eq. (14). But we shall show that a slight modi~cation of this function, as proposed in [9], leads in general to a function with a second derivative of fixed sign, the root of the modifkatii n equation being connected in a simple way with the root of Eq. (14). Notice that the method of proof used below is different from that used in [9]. In fact, consider the problem of minimizing the functional !I?,,[ u] = ]ILu - g]Isz + h/Au - fllm2, h > 0: Yv~]u”]= inf Yr,EUl,

(19)

PO@)= Mu” - film.

(20)

utn*

and we define the error function

It can easily be seen that the solutions u” of problem (19) are connected with the solutions us of problem (5) by r_P= u,, a = 1 /A, while the function p@(h) = ~(1 fh). Hence, if

pkd =A,

p&a) = A,

(211

where A > 0 is a given quantity, such that the first equation is uniquely solvable, then the roots ad, and hA of these equations will obviously be connected by aA = i/ ha. Theorem 4

For all A > 0 the function p. (A) is twice continuously differentiable, and pa’(h) G 0, @l”(k) 2 0. &oof: The solutions zzh of problem (19) are given by Euler’s equation L’(Lu” - g) + L4’(Au;” - f) = 0, the vector derivative du” i d31.by the equation

(22)

Error principle in the solution of incompatible equations

11

and the vector of second derivatives d’u” / dh” by the equation

(L’L + L4’A)

d2U" --=_-2A’Ag. dh2

(23)

Using Eq. (22) we rind from the definition (20) that

p,(h)po’(h) = (Au&4uk-~f)m= = -( (L’L+L4=d)-‘A’(Auh

(d,AT(h”-~f)).

(24)

- f), A’(Az2 - f))n < 0,

i.e., pa’(h) < 0 (here, UI’ = du” / dh). Further, using Eqs. (23) and (24)

But, by the Cauchy inequality,

po’“@)=

-L (A U&‘,AuL PO2 (A)

f)J2 fllhr’l17n2,

whence

po”@)po(h) 2 (uaN,A*(~uk - f))?l =2( (LTL+hATA) -‘ATA ( LTL+hATA) -‘AT(Au” -f)

, A’ (Auk-f)

) ,>O.

The theorem is proved. Note. If the matrix DLis not degenerate, we can set h > 0 in the conditions of the theorem. The process of solving (14) is considerably simplified by using this theorem. 7. Below, g = Lu*, u* E R,, is a given vector, and C = L’L is a positive definite matrix in R,. This case is common when the regularization method is used in practice. The solution of problem (19) is then the solution of the system of equations C(u -

U*) + hA’(Au

- f) = 0,

(25)

where the superscript T denotes tr~sposi~on, The tr~sfo~ations proposed by V. V. Voevodin [ 1Of may be used to reduce the number of arithmetic operations in these computations. As applied to f25), the first stage of the Voevodin tra~sfo~ations consists in performing the following operations: 1) we use the method af square roots to obtain the decomposition c = K’K, where K is an upper t~an~r

matrix,

2) we fmd the solution N of the matrix equation A’ = &?N; 3) we evahrate the vector y’ = Ku”. After performing the non-degenerate change of variables Ku = y, we arrive at the simpler system of equations in the vector

=o.

g-~*+XNT(Ny--~)

(25’)

If gR is a solution of Eq. (253, the solution of Eq. (25) will be obtained by solving the system

where

The transfo~ations

of the second stage are:

4) we use the method of reflections (or rotations) to obtain the decomposition N =:QDRz where Q and R are orthogonal na X 112and n X n matrices respectively, and D is a matrix of the form

IJ=

I3

I

. .. ,

[ 0

where 13 is an upper two-diagonal PIX n matrix, 5) we evahrate the vectors z* = Ry”, rp = Q’J’;setting z = Ry, we arrive at the system of equations z-z*+ilc)‘(Dz-~)

=O;

(25”)

Error principle in the solution of incompatible equations

if zh is the solution of the system (25”), the solution given by yi =

13

yr of the system (25’) will be

RX,

(27)

where

There is a standard sub-program [l l] . written in M-20 computer codes, for obtaining these decompositions (with u’ = 0),and evaluating zA and restoring the solution uh in accordance with Eqs. (26) and (27), for arbitrary h = I/ a, a > 0, Problem (19), with a parameter 3Lasatisfying condition (21), and C = E,u*= 0 was solved in [ 121 in FORTRAN on the BESM-6 computer. In this case the transformations of the first stage are omitted. A feature of these realizations is that the solution zh of (25”) is found explicitly. is We describe below a simple modification, whereby the error vector eh = Dzl- r.p found directly. This modification allows both the values of the function fro’(k) ,and the values of its derivative PC (A), required for numerical solution of the equation p. (A) = A to be evaluated to high accuracy. What we do is multiply Eq. (25”) on the left by the matrix and thus obtain the system of equations

D and set Dz'= a*,

(E+~DD')E=~'-cp,

(28)

the solution of which is the vector EA. To find khawe use the equation PO(~)

=

lIeAL =

A.

(29)

Notice that the solution z1 of system (25”) is connected with the solution er of system (28) by zA

=

--~D'E, + z'.

This is easily seen by substituting Eq. (30) in Eq. (25”) and recalling Eq. (28). In addition, ]I(E+ ADD')-'II < 1uniformly with respect to h > 0. The Equation (29) is solved numerically by Newton’s method of tangents, applied to the equation

(30)

14

V. A, Morozov

(31)

The advantages of taking the equation for 3ta in the form (3 1) were discussed in 19,121. Successive approximations to the root L, of Eq. (3 1) are found from the expression =h

h n+i

n

+P&) A

A-P&) Po’@)

’

where the initial value IQ is chosen in such a way that h, G L. The successive approximations h, are certainly convergent to ha as n -+ 00;this is a consequence of the convexity of the function 1/ p. (h) when L 3 0 (see Note on Theorem 4). With regard to the choice of the initial approximation ho we can obviously always take ho = 0. But it only needs a slight extra computational effort to discover a better approximation to the root ha of Eq. (31). The idea is to construct a simple minorant t”(h) for the function p. (3L)with h 3 0. We then naturally take as no the solution of the equation r(a) = A. We have

hDD’)-‘DD’En, Exfm-

po’(h)po(h) = -((lx+ But ( (E + ADD’) -‘DD’is,

6s) m d


11011”

PO2OJ * 1 + 7~llDll~

Consequently,

PO’(a) 3 -

lIDlIZ 1+

pa(O) =IIW* - qdl,.

PO@) t ?JlDl12

Sol~gt~s~equ~tywith h30, weget pa(h) 3?(h) The required ho is now easily found: O
poW--

= PO(O) (1 +hjlDlf2)-*.

A

AllDl12

-G ha.

Some points that arise in the course of the computations are worth mentioning. Noting that the matrix DD’has the form

Error principle in the solution of incompatibie equations

where 6

15

was ~troduced above, we set

Em),

cp= (cp,VnJ-s, . - * >cpm)

&a =

(En, F,,+t,

&a’=

(2, EL+*, . . . , Earn’),

. . . ,

7

0.

=<6, o,;i,. . G.h’), I

^

where EL, ^Eb’,o*, *‘rpare n-dimensional vectors. It can easily be seen that E&i= 0:’ - Cpj,

&Ail=O,

i>nfl,

and that the vectors EP. and a? are given by the “truncated” systems of equations (32)

The advantages of changing to Eqs. (32) are obvious. It should also be mentioned that the matrix E + A,DtiT is tridiagonal and positive definite. This ensures that the sweep method can be used to solve (32) and that it is numerically stable. The forward run, i.e., evaluation of the pivotal conden~tion coefficients, is obviously the same in both cases. We evaluate pa (h) and po’ (A) from the expressions II)

h

z

n

PO2 (h) = (Ea,Ed

+

(Oi*-(Pi)'7

i=n+i

po’ (A) po (A) = (2,

E^h)n*

With the value obtained for the parameter 3, = h, the solution of system (25”) is given by

or directly by the system of equations (E +

am)

z,q=

l$,‘cp+ z’,

the matrix of which is also tridiagonal and positive definite for all h 3 0. The numerical scheme described has a lot in common with the scheme for computing the lower bounds of functionals in [13]. In view of this, an algorithm can be developed for solving both problems s~ultaneously. T~~ns~~~d by

D. E. Brown

16

V. A. Morozov

REFERENCES 1.

TIKHONOV, A. N., The solution of incorrectly posed problems and the method of regularization, Dokl. Akad. NaukSSSR, 151,3,501-504, 1963.

2.

MOROZOV, V. A., The solution of functional equations by the regularization method, Dokl. Akad. Nauk SSSR, 167, 3,510-512, 1966.

3.

MOROZOV, V. A., The choice of parameter when solving functional equationsby Dokl. Akad. Nauk SSSR, 175,6, 1225-1228, 1961.

4.

MOROZOV, V. A., The error principle in the solution of operator equations by the method of regularization, Zh. vjjchisl.Mat. mat. Fiz., 8, 2, 295-309, 1968.

5.

MOROZOV, V. A., and KIRSANOVA, N. N., An extension of the method of regularization, in: ComputationalMethods and Programming (Vj%hisI. metody i programmirovanie), 14, 40-45, MGU, Moscow, 1970.

6.

ZHIDKOV, E. P., and LUK’YANTSEV, A. F., A method for finding constrained extrema, OIYaI Preprint, Dubna, 1965.

I.

IVANOV, V. K., The approximate solution of operator equations of the first kind, Zh. vjjchisl. Mat. mat. Fiz., 6, 6, 1089-1094, 1966.

8.

MOROZOV, V. A., The solution by regularization of incorrectly posed problems with non-linear unbounded operators, Differents. ur-niya, 16, 8, 1453-1458, 1970.

9.

GORDONOVA, V. I., and MOROZOV, V. A., Numerical algorithms for parameter selection in the method of regularization, Zh. vphisl. Mat. mat. Fiz., 13, 3, 539-545, 1973.

regularization,

10. VOEVODIN, V. A., The method of regularization, Zh. vjjchisl.Mat. mat. Fiz., 9, 3,673-615, 1969. 11. VOLOVICH, V. M., and NOVIKOVA, I. V., A standard program for solving systems of linear algebraic equations by the method of regularization, 33, MGU, Moscow, 1968. 12. Analyses of numerical algorithms for parameter selection in the method of regularization, Scientific Report No. 161-T3(507), VTs MGU, 1972. 13. MOROZOV, V. A., Evaluation of the lower bounds of functionah on the basis of approximate information, Zh. vychisl. Mat. mat. Fiz., 13,4, 1045-1049, 1973.