Improving Tikhonov Regularization with Linearly Constrained Optimization: Application to the Inverse Epicardial Potential Solution ILIAS IAKOVIDIS

AND

Institute of Biomedical

RAMESH

Engineeting,

M. GULRAJANI University of Montreal,

Montreal, Quebec H3C 3J7, Canada Received 6 December

1991; revised 30 March 1992

ABSTRACT Two methods to improve on the accuracy of the Tikhonov regularization technique commonly used for the stable recovery of solutions to ill-posed problems are presented. These methods do not require a priori knowledge of the properties of the

solution or of the error. Rather they exploit the observed properties of overregularized and underregularized Tikhonov solutions so as to impose linear constraints on the sought-after solution. The two methods were applied to the inverse problem of electrocardiography using a spherical heart-torso model and simulated inner-sphere (epicardial) and outer-sphere (body) potential distributions. It is shown that if the overregularized and underregularized Tikhonov solutions are chosen properly, the two methods yield epicardial solutions that are not only more accurate than the optimal Tikhonov solution but also provide other qualitative information, such as correct position of the extrema, not obtainable using ordinary Tikhonov regularization. A heuristic method to select the overregularized and underregularized solutions is discussed.

1.

INTRODUCTION

It is of great interest in cardiology to noninvasively recover information about cardiac electric events using measurements of the electric potential on the body surface. One way to achieve this is to compute the electric potentials on the epicardial surface. Calculating epicardial potentials using as data the body surface potentials constitutes the so-called inverse problem of electrocardiography. This particular inverse problem has been studied extensively for the past two decades, and an excellent review has been written by Rudy and Messinger-Rapport [ll]. The main difficulty in recovering a satisfactory inverse solution is the inherent instability of the inverse problem that results from the smoothing property of the intervening conducting medium (the body). This instabilMATHEMATICAL

BIOSCIENCES

112:55-80

(1992)

OElsevier Science Publishing Co., Inc., 1992 655 Avenue of the Americas, New York, NY 10010

55 0025-5564/92/$5.00

56

ILIAS IAKOVIDIS

AND

RAMESH

M. GULRAJANI

ity of the inverse problem manifests itself in the form of large oscillations in the values of the inverse solution even in the presence of very small errors in data and geometry. Several techniques to overcome this instability and to restore a continuous dependence of the solution on the data have been reported in the literature [9, 141. Here we will consider a widely used technique introduced by Tikhonov [13], namely, the zero-order Tikhonov regularization method. The accuracy of inverse solutions obtained by this method (Tikhonov solution) depends on the value chosen for a regularization parameter. The optimal value of this parameter that results in an Z2-closest solution to the actual epicardial potentials has been the subject of many studies [2, 91. As in practice the actual epicardial potentials are not known, the value of the regularization parameter is usually calculated by some empirical method [2]. It is also clear that if any a priori knowledge about the data or geometry errors is available, for example in the form of the statistical characteristics of the data noise, it can be used to recover a quasi-optimal value of the regularization parameter. Moreover, a priori information about the solution can also be used to impose constraints on the sought-after solution and consequently to even improve on the accuracy of the optimal Tikhonov solution [7]. The objectives of this paper are (1) to derive constraints in the form of linear bounds on the values of the solution without any a priori knowledge of either the data or geometry errors or of any particular characteristics of the epicardial potentials, (2) to present a method of calculating the inverse solution satisfying such constraints and to show that it is possible to obtain more accurate solutions than the unconstrained Tikhonov method can provide, and (3) to test this method on a concentric spheres model representing the heart-torso system. In the first part of this work we discuss the Tikhonov zero-order regularization method and derive a particular closed form of the Tikhonov solution. We examine the question of the existence and uniqueness of an optimal value of the regularization parameter and show how the data and geometry errors affect the solution. We also briefly present the relationship between the Tikhonov regularization method and the quasi-solution method introduced by Ivanov 161. In particular, it is shown that the two methods are equivalent, and the relationship between their respective parameters is presented. In the second part we discuss the properties of the Tikhonov solution that correspond to different values of the regularization parameter. On the basis of these properties we describe the derivation of linear bounds on the values of the sought-after solution. These bounds form the solution set to which the search for a solution is restricted. Specific parameters needed to define such a solution set are given heuristically on the basis

IMPROVING

TIKHONOV

REGULARIZATION

57

of available simulation and experimental data. Next we present two methods of obtaining an inverse solution from the previously constructed solution set and discuss the strengths and weaknesses of each method. It is shown that both methods can result in more accurate solutions than the optimal Tikhonov solution, provided that the solution set is chosen properly. Finally, we present a concentric spheres model that represents the heart-torso system. The inverse problem is then solved numerically using both the Tikhonov method and our own proposed method, and the accuracy in each case is tested with analytically obtained values. In what follows we treat the inverse problem in its discrete form; that is, we assumed that the epicardial potentials are represented by an n-vector & E R” and the body surface potentials by an m-vector & E R”, m 2 n. The operator that maps the epicardial potentials to body surface potentials is represented by an m x n transfer matrix A, and the inverse problem takes the form A&

= 4~.

(1)

Thus the inverse problem is reduced to solving the above matrix equation for &, with A being an ill-conditioned matrix. Finally, we denote by II.11 the 1’ norm in R” and by Ilx - yll the Z* distance between the n-vectors x and y. 2.

TIKHONOV

REGULARIZATION

In this section we consider the zero-order Tikhonov regularization method [14], where the solution to (1) is given as an unconstrained minimum in R” of the functional

for some value of a regularization parameter t > 0. It has been shown 1141 that for a given positive value of the regularization parameter t there exists a unique minimum; that is, there exists a C&such that

We call the solution c$~ a Tikhonov solution, and the subscript will indicate the dependence of c$, on the choice of the parameter t. It can be easily shown that +l has the form

(4)

58

ILIAS IAKOVIDIS

AND

RAMESH

M. GULRAJANI

where AT and Z denote the transpose of A and the n x n identity matrix, respectively. Thus, using Tikhonov zero-order regularization, the problem of recovering a meaningful solution to (1) is reduced to finding an optimal value of the regularization parameter t. Here a value fort of the regularization parameter is called optimal if the corresponding Tikhonov solution $Lloptsatisfies

The above definition of the optimal parameter requires a priori knowledge of the epicardial vector &, which, in practice, is unknown. Therefore, some algorithm (estimator) must be used that estimates a quasi-optimal value of the parameter t without knowledge of &. The goal of an estimator is to select a quasi-optimal value of t that lies close to tort. One such method for estimating a quasi-optimal value of t, based on a knowledge of the standard deviation of the noise level affecting the signal and the level of error in approximating A, is known as the generalized discrepancy criterion [3, 91. Such information is not always available in practical applications, and therefore other empirical methods for choosing the parameter t without knowledge of the characteristics of the noise affecting the data have been introduced. An empirical method called the composite residual and smoothing operator (CRESO) was presented in [2] and consists of finding the smallest t > 0 that results in the first relative maximum of the function C(t)

=

d

llc$,l12 +2td,l14,112.

The function C(t) is a derivative of the function B(t), where

B(t) = tll(b,l12- IlAd+ - 4112. and it follows that the quasi-optimal value t, calculated by this method corresponds to the first point where B(t) changes concavity. It was shown in [2] that CRESO results in values of t for which the corresponding Tikhonov solution 4, was more accurate than Tikhonov solutions using values of t obtained by two other empirical methods. Note that in the above discussion about the estimators it is implicitly assumed that an optimal parameter top, exists, which is not necessarily true. Next we briefly examine this question of the existence and uniqueness of topt. Let ATA be an IZx it symmetric and positive definite matrix. Denote by A, > A, > ***> A, the eigenvalues of ATA and by Ul,U 2,. . . , u, the corresponding orthonormalized system of eigenvectors.

IMPROVING The

TIKHONOV

59

REGULARIZATION

vector A’& can then be represented

as

A’~~ = ~ PjUj.

(5)

j=l

Substituting (5) into (4) and using the fact that the eigenvalues j.Q < 112< .‘. < /-L, of (ATA + tZ)F’ satisfy 1 pj = Aj + t 9

j=1,2

,.*., n,

it follows that

In the case when some of the eigenvalues of ATA are zero, Aj = 0, I< j G IZ,or equivalently when the normalized eigenvectors ui, u2,. . . , ul form an orthonormal basis of R’ c R”, we consider an extension uI+l,..., u,, such that u~,u~,..., u, forms an orthonormal basis of R” and take pj=O, j=l+l,..., II. Now, from the representation of the vector of epicardial potentials,

4E =

L

(yi"j7

(7)

j=l

and noting that the usual inner product in R” is independent of the choice of an orthonormal basis, we obtain the following form for the relative error of the solution C&with respect to the epicardial vector &

In the presence of noise in the body surface data and/or geometry errors, the coefficients pi are perturbed from their true values pi’. The eigenvalues Ai, on the other hand, are perturbed from their true values A; only in the presence of geometry errors. Let pj = pi’ + ej and Aj = A; + hi. It is immediately evident from (8) that for a zero epicardial vector, that is ‘Ye= 0 for all j, and data noise and/or geometry errors such that ej # 0 for some j, the function RE(t) goes to zero monotoni-

60

ILIAS IAKOVIDIS

AND

RAMESH

M. GULRAJANI

tally as t tends to infinity, and hence there does not exist a finite value for the optimal regularization parameter top*. Similarly, if the data noise and/or geometry errors ej result in pj such that aj /3, < 0 for all j, then again there does not exist a finite value of top*. Thus, the existence and uniqueness of topt depends on the characteristics of the data and geometry errors as well as on &. However, on the basis of experimental and simulation data where & was known (e.g., see [2], [81, [121), it can be heuristically shown that there exists a unique top;. Furthermore, it was observed that RE(t) is strictly concave up on the interval 10,y,), which contains topt, and strictly concave down and an increasing function for t > y,,, bounded from above by 1, because it is known [91that lim c#+= 0. f+” Now let us examine whether the Tikhonov solution c#+can recover exactly a nonzero epicardial vector & in the presence of nonzero error vectors (ej)TZ1 and (hi),?=,. Note that in the absence of errors in data, geometry, and discretization, the coefficients aj, pi’, and the eigenvalues A; are related as follows: ffj/4; = p;,

j=1,2

,..., n.

(9)

Setting RE(t) = 0, Equations (8) and (9) yield t+hj-;=O

when oj # 0 J

(10)

and ej = 0

when crj =O.

(11)

In a practical situation it is very unlikely that the vectors (ej);= ,,(hj)yz, will satisfy (10) and (11); hence RE(t) > 0 for all t > 0. Because even the is not exactly equal to &, a natural optimal Tikhonov solution & question that immediately aris& is whether it is possible to improve the accuracy of 4, without any additional information about the data and geometry errog or about the solution itself. In this paper we deal with exactly this question. We show that an improvement can be attained if the solution is sought in a specific solution set whose construction is based on certain properties of an overregularized Tikhonov solution 4 lo, where t, > topt. These properties of an overregularized solution and the construction of the solution set will be described in detail in the following section.

IMPROVING

TIKHONOV

61

REGULARIZATION

In order to see clearly how a proper choice of a solution set can lead to an improved Tikhonov solution, we first demonstrate the equivalence between a Tikhonov solution and a quasi-solution as defined by Ivanov [6]. The quasi-solution & of problem (1) is defined as &:

llA&

- 4Bl12 = min lM4

@EXR

-

&l12,

(14

where X, is a set of vectors whose norm does not exceed R, X, = { 4 E R”I lIdI< R}.

It was shown by Ivanov [6] that the quasi-solution

has the form

(13) where Ai, uj, and fij are defined as before, and t > 0 is the solution of

(14) Thus, we clearly see that the Tikhonov zero-order solution 4, and a quasi-solution & are identical provided that t and R satisfy relation (14). The above observation defines yet another method of choosing the parameter t given that a sufficiently accurate estimate of II&l] is known. An analytic expression for the quasi-solution & in the case of spherical geometry, together with bounds on II& - &II, has been derived by Iakovidis and Martin [4]. It is also worthwhile to note that by using representation (6) of the solution +t, it is not necessary to invert the matrix (ATA + tZ) to obtain #+ for different values of the parameter t. Consequently, the computational time needed to calculate $* for many different values of the parameter t can be significantly reduced. Now let us return to the question of the proper choice of a solution set X. From the above discussion it follows that the optimal Tikhonov solution is a solution of the constrained least squares problem given by Equation (12) for some choice of the solution set X, that represents a ball in R” with radius R. Suppose additional regional information is available that can be used to impose further constraints on the solution. In this way the ball X, will be “deformed,” for example, into an n-dimensional parallelogram, by using linear bounds on the values of the solution. Provided that these bounds are chosen correctly, the solution from this new solution set will be more accurate than the

62

ILIAS IAKOVIDIS

Tikhonov solution. obtaining regional on the values of information about 3.

AND

RAMESH

M. GULRAJANI

In the following section we present a method for information about the solution in the form of bounds the solution with no additional data or statistical the errors or the solution.

CONSTRUCTION

OF A SOLUTION

SET

In the previous section it was shown that the Tikhonov solution 4, cannot recover the epicardial values & in the presence of error vectors (ej);2=,,(hj)~= 1 that do not satisfy (10) and (11). Consequently, in what follows it is assumed that RE(t), defined as in (8), is bounded from below, RE( t) > 6 > 0

for all t > 0.

In other words, there exists a S neighborhood

of &,

(15) that does not contain a Tikhonov solution 4, for any value of t. In the case of a unique value of topt we can take 6 = RE(t,,,) > 0. Thus, we observe that to obtain a solution to the inverse problem (1) within B, requires either the use of a different method or knowledge of some characteristics of the solution. These characteristics lead to the formulation of a solution set X from which the solution is to be selected. For example, if both upper and lower bounds are known for each value +$ of the epicardial vector &, a solution set X can be constructed as X={~ER”la,~~‘~b,,i=1,2

)...) n}.

(16)

If the bounds ai and bi are chosen correctly, that is, if X fl B, = V # 0, then each element $,, E I/ is an example of a solution whose relative error is less than that of the optimal Tikhonov solution. A solution set of the form (16) was considered in [7], where the bounds ai and b, were chosen on the basis of known epicardial values &. The inverse solution was taken to be an element of X that minimizes Ml, and it was shown that when narrow bounds around the epicardial values were used the relative error was significantly improved in comparison to the unconstrained Tikhonov solution 4,. These types of narrow bounds on the solution are not available in a practical setting. Therefore, it is desirable to develop methods that would use a posteriori information about the solution, that is, information that can be extracted by examining the characteristics of inverse solutions that are calculated either by different methods or for different values of the regularization parame-

63

IMPROVING TIKHONOV REGULARIZATION

ter. Here we present an empirical method for constructing a solution set X of the form (16) that exploits certain properties that Tikhonov solutions exhibit for different values of the regularization parameter t. From the form (6) of the Tikhonov solution c#+we can observe that the errors ej in the coefficients pj are magnified for small values of the parameter t and small values of Ai by the factor ( Aj + t)- ’ , which results in large oscillations in the solution 4,. Therefore, one of the first effects that an increased value of parameter t has on the solution is the reduction of oscillations associated with the small eigenvalues of A%. Thus, an increase in the value of the parameter t results in smoother Tikhonov solutions with reduced amplitudes and gradients. Moreover, these smoother solutions tend to recover more accurately the distribution of the positive and negative values of & or, equivalently, to recover more accurately the zero line of c#+. Such dependence of the Tikhonov solution has been observed in experimental data [121, where the underregularized solution c#+,,,, corresponding to some t, G topt, recovered more accurately the magnitudes and location of the extrema whereas the overregularized solution c#+, corresponding to some t, > t opt, was more accurate in recovering the zero line of &. We use the above property of an overregularized solution to construct the solution set X. First, let us denote by S, t,, such that IS,1 < IS,uI. Now, for E and K positive, such that K > E, we have X={~ERRlai~~igb,,i=1,2

,..., n},

where bi = K

when c#(> E > 0,

ai = - K,

b,=-E

when 4i0