A problem frequently encountered in the earth sciences requires deducing physical parameters of the system of interest from measurements of some other (hopefully) closely related physical quantity. The obvious example in seismology (either surface reflection seismology or crosswell seismic tomography) is the use of measurements of sound wave traveltime to deduce wavespeed distribution in the earth and then subsequently to infer the values of other physical quantities of interest such as porosity, water or oil saturation, permeability, etc. The author presents and discusses some general ideas about iterative nonlinear output least-squares methods. The main result is that, if it is possible to do forward modeling on a physical problem in a way that permits the output (i.e., the predicted values of some physical parameter that could be measured) and the first derivative of the same output with respect to the model parameters (whatever they may be) to be calculated numerically, then it is possible (at least in principle) to solve the inverse problem using the method described. The main trick learned in this analysis comes from the realization that the steps in the model updates may have to be quite small in some cases for the implied guarantees of convergence to be realized

Berryman, J.G.

UNT Digital Library

DISCLAIMER This document was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor the University of California nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise, does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or the University of Califomia. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or the University of Caliiomia, and shall not be used for advertising or product endorsement purposes. Nonlinear Least Squares and Regularization James G. Berryman ABSTRACT I present and discuss some general ideas about iterative nonlinear output least-squares methods. The main result is that, if it is possible to do forward modeling on a physical problem in a way that permits the output ( i e . ,  the predicted values of some physical parameter that could be measured) and the first derivative of the same output with respect to the model parameters (whatever they may be) to be calculated numerically, then it is possible (at least in principle) to solve the inverse problem using the method described. The main trick learned in this analysis comes from the realization that the steps in the model updates may have to be quite small in some cases for the implied guarantees of convergence to be realized. INTRODUCTION A problem frequently encountered in the earth sciences, and in other physical and biomedical sciences as well, requires deducing physical parameters of the system of interest from measure- ments of some other (hopefully) closely related physical quantity. The obvious example in seismology (either surface reflection seismology or crosswell seismic tomography) is the use of measurements of sound wave traveltime to deduce wavespeed distribution in the earth and then subsequently to infer the values of other physical quantities of interest such as porosity, water or oil saturation, permeability, etc. Many of the problems of interest can be formulated in a way such that the measured quantities (often called the “output)’ for reasons that will become clear in a moment) may be compared to predicted values of those same quantities and the resulting observed discrepancies then used to make “improvements)’ in the system model parameters of real interest. Predicted values are obtained by forward modeling based on some assumed model of the physical quantities of interest, and these predicted values are clearly the “output” of such a forward modeling code. Comparisons between predicted output and measured output parameters may be done in a variety of ways, but a common method is the output least-squares method: discrepancies are squared and summed, and some numerical procedure is established to  reduce the overall least-squares error in the output quantity. The trick in all this is to find a method that actually does guarantee convergence of the output least-squares functional to zero, or at least to a small number whose size is comparable to that expected from the sums of squares of the measurement errors. The purpose of this paper is to show that such a procedure can essentially always be constructed as long as one additional condition is present: If the nonlinear functional of the model parameters used to compute the outputs can be differentiated once with the respect to each of the model parameters of interest (whether this derivative is taken analytically or numerically does not appear to  be important), then an evolution equation can be found that will essentially guarantee that a sequence of models gradually improving the agreement between the measured data and the predicted data can be found in a systematic way. The only caveat is that the step size of the improvement from one iteration to the next will be problem dependent and may be rather small in some applications of interest. NONLINEAR OUTPUT LEAST-SQUARES METHOD In an abstract setting, let {d i }  for i = 1,. . . , m be a set of m measurements. This data may be organized into a data m-vector d. Let s be an n-vector containing the n parameters for the modellsolution space to be constructed. Then, let N i ( S )  be a (possibly nonlinear) function of the model parameter vector s such that Ni(s) = d; in infinite precision forward modeling of the physical problem of interest. If the data d; and the functions Ni are known, then the inversion problem is to find an approximate s that satisfies the data ( i .  e., has as small a data misfit as possible) in some sense. If life were simple, the operators Ni would be analytically invertible, so one might think it could be possible to write the solution to the inverse problem as s = Nrl(di)s  However, the result is certainly a strange looking formula: While it is perfectly sensible that the domain of each function Ni is a set of n model parameters, and that the range is a single datum, the postulated inverse function with multiple model outputs from single datum inputs is clearly 1 not sensible. In most cases, there will simply not be sufficient information contained in a single datum to determine multiple model parameters. For strictly linear problems, there might exist an appropriate operator X acting on the full data m-vector such that X(d) = s produces the model n-vector. However, even for linear problems, the existence of such an operator is not necessarily guaranteed since the sizes m and n of the two vector spaces will generally differ. This difficult situation leads to the introduction of pseudoinverses in linear problems or to the methods to be described now in nonlinear problems. So, instead of using some hypothetical analytical approach, the desired solution to the inverse problem will make simultaneous use of all (or at  least most) of the data while determining some sort of optimal fit to a chosen set of model parameters. This situation is common in data analysis and often leads to formulation of least-squares methods for this type of inversion problem. An approach to this problem based on nonlinear least-squares inversion considers the non- linear function ’ which is simply the sum of the squares of the data residuals T i  = d; - Ni. I formulate an “evolution” equation for the model vector s = s ( q ) ,  where q is an evolution parameter - treated as a continuous generalization of an iteration number. In particular, I must ultimately discretize this parameter in my final numerical method, in which case the particular values of q = 0, 1 , 2 , .  . . will be exactly the iteration numbers. But for purposes of explaining the method, it will prove useful to treat q initially as a continuous parameter. The thrust of the nonlinear least-squares method is to produce a model s that minimizes the least-squares error function F(s) .  One way to guarantee [see, for example, Jeffrey and Rosner (1986) and Lu and Berryman (1991)l that a local minimum is achieved is to pose the evolution equation for s(q) in a way that guarantees the value of F(s )  decreases monotonically as q increases (or at  each iteration step for the discretized problem). Taking the derivative of F ( s )  with respect to the evolution parameter, the chain rule gives dN;  ds j  677 j=1 i=l dsj dq n m  -- aF(S)  - -2 E p i  - Ni(S)]--. It is desired that the evolution equation for s be chosen to guarantee that (2) is 5 0. It is easy to see that (2) will always be negative or zero if the evolution equations for the model parameters s j  are chosen to be where 7(q) is a positive parameter to be determined. Then, 2 SO each term in the sum over j in (4) is a square quantity and must be positive or vanish identi- cally. The choice (3) of evolution then clearly guarantees the desired monotonically decreasing behavior of the total squared error functional as q 3 00. To implement this procedure, a choice of discretization must be made together with a choice for y(q) .  The obvious choice of discretization for the evolution equation (3) is where I have taken the finite step size for the evolution to be A7 = 1. In the continuous evolution problem, the infinitesimal changes in the evolution parameter guarantee similarly infinitesimal changes in the model vector s ( q )  and therefore this makes the choice of y(7) largely arbitrary. In contrast, for the discretized problem, the evolution of s is finite at  each step and care must be taken not to violate the desired condition that the least-squares functional should decrease at  each step. Such violations may occur if the step size is taken too large. Reconsidering (l), I find that, by keeping only those terms proportional to the first and second powers of the components of As, I have m m After substituting ( 5 ) ,  the parameter ~ ( q )  can now be chosen so that the right-hand side of (6) decreases as much as possible at  each step of the iteration scheme. The optimum choice is easily shown to be y(b + 1) = aN. aN. aN. (7) E;=, [EL [di - N~(s(~))I*] aN. Ej,jt [ET=l[dit- ~ i t ( s ( k ) ) ~ g  C Z ~  e~ C$=l[dilt- ~it#(s(b))]%] ' since this minimizes the right-hand side of (6). It will prove enlightening to compare the procedure just presented with the well-known method of conjugate gradients (Hestenes and Stiefel, 1952; Fomel, 1996) for a linear operator such that N j ( S )  = (LS)i. Then, model updates are obtained using #+I) = p + y(k+')U(k+l), (8) where, for the discrete iterative problem, I put the iteration numbers in superscripts. The new vector u ( ~ + ~ )  is the (somehow) known update direction in the model space and 7(k+') is a parameter used to optimize the step size. The updated residual vector is again given by (k+l)LU(k+'). (9) r(k+l)  = d - L s ( k + l )  - r(k) - 7 The magnitude of the residual vector is easily shown to decrease most at each step of the iteration sequence if the optimization parameter satisfies 3 Equation (10) is exactly analogous to the formula (7) obtained in the nonlinear least-squares problem if the components of the residual are given by .lk) = d; - N;(s(')), (11) the matrix elements of the linear operator are and the update direction vector satisfies where the components of the right-hand side of (13) are given by (3) evaluated at With this identification, it becomes clear that the nonlinear least-squares method outlined above is one natural generalization of the well-known conjugate gradients technique. REGULARIZATION Now it may happen, due to a combination of poor conditioning of the linear operator L and unfortunate choices of starting guess of the model vector s, that the denominator of the right- hand side of (10) may be very small or vanish to numerical accuracy. Similar circumstances can arise in the nonlinear least-squares problem in cases of sparse or irregularly sampled data. When such circumstances arise in practice, it may be necessary to regularize the method by adding an additional constraint equation to the least-squares functional (1). Regularization is a well-known technique often associated with the names of Tikhonov and Arsenine (1976), Levenberg (1944), and Marquardt (1963,1970), among others. The regularization constraint usually takes the form of a quadratic functional of the model vector. One typical choice is p ( s  - I)T(s - 3) where p (often called the damping parameter) is some small positive constant and I is some value of the model vector that the final solution should not deviate from too much. Other typical choices of regularization constraint are based on differentials of the model taking the form (Ds)~(Ds), where Ds might be either a simple gradient or a Laplacian of the model - assuming that s is some simple physical quantity. If s is a more complicated vector of model parameters which is not easily given a simple physical interpretation, then some other choice of regularization constraint might be needed. It is preferable to avoid using regularization if possible, because such techniques tend to modify the entire spectrum of the operator to be inverted and therefore tend to degrade reso- lu t ion. For the linear least-squares problem, I can carry the analysis further to understand how the iteraton scheme modifies the model estimate at each step. Assuming that the linear operator L is a matrix of dimensions m x n, where m and n are usually not equal, it is helpful to use L and its transpose LT to form a square, symmetric matrix and the associated eigenvalue problem 4 Assuming that the rn-vectors +q and the n-vectors 4q are normalized to unity, the singular value decomposition of the matrix L may now be written in terms of these eigenfunctions as The sum is taken over T terms, where T is the rank of L. according to Now, each model estimate may be expanded in terms of the appropriate eigenfunctions t q=l where the aik)s are the expansion coefficients for the kth iteration. Similarly, the data vector may also be expanded as where the 6,s are constant coefficients. Both the model vector and the data vector expansion might in addition include a term from the null space of L, but for simplicity I will ignore this possibility for the present purposes. With these definitions of the various coefficients, I find that equation (10) becomes and the update equation for the a q s  becomes It follows from (19) that the iteration process h& converged when As the coefficients approach convergence, I discover that the denominator of (18) can become quite small. If1 assume that the eigenvectors with the largest eigenvalues converge most quickly, then after some number of iterations the main contributions to the denominator will be from the terms associated with the smallest eigenvalues, and these contributions are proportional to the fourth power of these small eigenvalues. If this happens, some type of regularization may be required to obtain useful results. Consider the simplest type of regularization involving a model vector constraint so that the modified objective function becomes: m n i= l  j=1 5 . The expansion of the constraint vector is given by r s = CC5,Qq. q=1 The equation for the modified model update is then given by Clearly, the modified iteration sequence for aq has converged when showing that the limiting relation is just a linear combination of the starting value and the result for pure least-squares as seen in (20). If the damping parameter p is small, these two values will of course be quite close. But, in general, for any nonzero value of the damping parameter, it is inevitable that the resolution will suffer due to the fact that the coefficients cannot approach their optimum value (20) but are actually constrained away from it by the regularization procedure. This is why regularization should be avoided, or at least minimized, as much as possible. SUMMARY AND SUGGESTIONS FOR FUTURE WORK The main result obtained in this paper is this: If it is possible to do forward modeling on a physical problem in a way that permits the output (Le., the predicted values of some physical parameter that could be measured) and the first derivative of that same output with respect to model parameters (whatever they may be) to be calculated numerically, then (at least in principle) it is possible to solve the inverse problem using the method described. The main trick learned in this analysis comes from the realization that the steps in the model updates may have to be quite small in some cases for the implied guarantees of convergence to be realized. We have concentrated on least-squares methods in the presentation in order to keep the analysis simple. However, it is clear that the general method presented could be applied as well to any norm of the data discrepancies that possesses the same necessary qualities, the most important of these being the existence of a first derivative with respect to the model parameters. One interesting example of such an alternative is Huber’s hybrid of I1 and Zz [see Clearbout (1996)l. Other choices include special cases of the more general Z, norm that also have a first derivative. These other methods will be discussed in future work. ACKNOWLEDGMENTS I want to thank Liliana Borcea and George Papanicolaou for numerous discussions of inverse problems and the use of nonlinear output least-squares methods. Work of JGB was performed under the auspices of the U. S. Department of Energy by the Lawrence Livermore National Laboratory under contract No. W-7405-ENG-48 and supported specifically by the Geosciences Research Program of the DOE Office of Energy Research within the Office of Basic Energy Sciences, Division of Engineering and Geosciences. 6 REFERENCES Claerbout, J., 1996, Conjugate-direction Huber regression: SEP-92, 229-236. Fomel, S., 1996, Least-square inversion with inexact adjoints. Method of conjugate directions: A tutorial: SEP-92, 253-266. Hestenes, M. R., and E. Stiefel, 1952, Methods of conjugate gradients for solving linear systems: J .  Res. NBS49,  409-436. Jeffrey, W., and R. Rosner, 1986, Optimization algorithms: Simulated annealing and neural network processing: Astrophys. J. 310, 473-481. Levenberg, K., 1944, A method for the solution of certain non-linear problems in least squares: Quart. Appl. Math. 2, 164-168. Lu, S.-Y., and J. G. Berryman, 1991, Inverse scattering, seismic traveltime tomography, and neural networks: Intern. J. Imaging Sys. Tech. 2, 112-118. Marquardt, D. W., 1963, An algorithm for least-squares estimation of nonlinear parameters: SIAM J.  Appl. Math. 11, 431-441. Marquardt, D. W., 1970, Generalized inverses, ridge regression, biased linear estimation, and nonlinear estimation: Technometrics 12, 591-612. Tikhonov, A., and V. Arshnine, Mdthodes de Rdsolution de Problimes Mal Posds, Editions Mir, Moscow, 1976. -- . . -  . .  . .  _ . _ .  . _  _ .  . .. Recycled Recyclable . . .  . . . .. . .  . . .  . . .  . . .  . .  ... . .  . .. .. . . . . ... - .. . .  . 

Nonlinear least squares and regularization

https://digital.library.unt.edu/ark:/67531/metadc672961/m2/1/high_res_d/238461.pdf

Nonlinear least squares and regularization

Abstract

Similar works

Full text

Available Versions

UNT Digital Library