f18-4

18.4InverseProblemsandtheUse ofA PrioriInformation 795Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).18.4 Inverse Problems and the Use of A Priori Information Later discussion will be facilitated by some preliminary mention of a couple of mathematical points. Suppose that uis an “unknown” vector that we plan to determine by some minimization principle. Let A[u]>0andB[u]>0be two positive functionals of u, so that we can try to determine uby either minimize: A[u]or minimize: B[u]( 18.4.1 ) (Ofcoursethese will generallygivedifferentanswers for u.) As anotherpossibility, nowsupposethat we want to minimize A[u]subject to the constraint thatB[u]have someparticularvalue,say b. ThemethodofLagrangemultipliersgivesthevariation δ δu{A[u]+λ1(B[u]−b)}=δ δu(A[u]+λ1B[u]) = 0 ( 18.4.2 ) where λ1is a Lagrange multiplier. Notice that bis absent in the second equality, since it doesn’t depend on u. Next, suppose that we change our minds and decide to minimize B[u]subject to the constraint that A[u]have a particular value, a. Instead of equation (18.4.2) we have δ δu{B[u]+λ2(A[u]−a)}=δ δu(B[u]+λ2A[u]) = 0 ( 18.4.3 ) with, this time, λ2the Lagrange multiplier. Multiplying equation (18.4.3) by the constant 1/λ 2, and identifying 1/λ 2with λ1, we see that the actual variations are exactly the same in the two cases. Both cases will yield the same one-parameter family of solutions, say, u(λ1).A s λ1varies from 0to∞, the solution u(λ1) varies along a so-called trade-off curve between the problem of minimizing Aand the problem of minimizing B. Any solution along this curve can equally well be thought of as either (i) a minimization of Afor some constrained value of B, or (ii) a minimization of Bfor some constrained value of A, or (iii) a weighted minimization of the sum A+λ1B. Thesecondpreliminarypointhastodowith degenerate minimizationprinciples. In the example above, now suppose that A[u]has the particular form A[u]=|A·u−c|2(18.4.4 ) forsomematrix Aandvector c.I fAhasfewerrowsthancolumns,orif Ais square but degenerate (has a nontrivial nullspace, see §2.6, especially Figure 2.6.1), then minimizing A[u]willnotgive a unique solution for u. (To see why, review §15.4, and note that for a “design matrix” Awith fewer rows than columns, the matrix AT·Ain the normal equations 15.4.10 is degenerate.) However, if we add any multiple λtimes a nondegeneratequadraticform B[u], forexample u·H·uwithH a positive definite matrix, then minimization of A[u]+λB[u]willlead to a unique solution for u. (The sum of two quadratic forms is itself a quadratic form, with the second piece guaranteeing nondegeneracy.) 796 Chapter18. IntegralEquationsandInverseTheorySample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).We can combine these two points, for this conclusion: When a quadratic minimizationprincipleiscombinedwithaquadraticconstraint,andbotharepositive,onlyoneofthetwoneedbenondegeneratefortheoverallproblemtobewell-posed. We are now equipped to face the subject of inverse problems. The Inverse ProblemwithZeroth-OrderRegularization Supposethat u(x)is someunknownorunderlying( ustands for bothunknown and underlying!) physical process, which we hope to determine by a set of N measurements ci,i=1,2,...,N. The relation between u(x)and the ci’s is that eachcimeasuresa(hopefullydistinct)aspectof u(x)throughitsownlinearresponse kernel ri, and with its own measurement error ni. In other words, ci≡si+ni=/integraldisplay ri(x)u(x)dx+ni (18.4.5 ) (compare this to equations 13.3.1 and 13.3.2). Within the assumption of linearity, this is quite a general formulation. The ci’s might approximate values of u(x)at certain locations xi, in which case ri(x)would have the form of a more or less narrowinstrumentalresponsecenteredaround x=xi. Or,the ci’smight“live”inan entirelydifferentfunctionspacefrom u(x),measuringdifferentFouriercomponents ofu(x)for example. Theinverseproblem is,giventhe ci’s,the ri(x)’s,andperhapssomeinformation about the errors nisuch as their covariance matrix Sij≡Covar [ni,n j]( 18.4.6 ) how do we find a good statistical estimator of u(x), call it /hatwideu(x)? It should be obvious that this is an ill-posed problem. After all, how can we reconstruct a whole function /hatwideu(x)from only a finite number of discrete values ci? Yet,whetherformallyorinformally,we dothis all thetimeinscience. We routinely measure “enough points” and then “draw a curve through them.” In doing so, weare making some assumptions, either about the underlying function u(x), or about thenatureofthe responsefunctions r i(x), orboth. Our purposenowis toformalize these assumptions, and to extendour abilities to cases where the measurements and underlying function live in quite different function spaces. (How do you “draw a curve” through a scattering of Fourier coefficients?) We can’t really want every point xof the function /hatwideu(x). We do want some large number Mof discrete points xµ,µ=1,2,...,M, where Mis sufficiently large, and the xµ’s are sufficiently evenly spaced, that neither u(x)norri(x)varies muchbetweenany xµandxµ+1. (Hereandfollowingwe will use Greekletters like µto denote values in the space of the underlying process, and Roman letters like i to denote values of immediate observables.) For such a dense set of xµ’s, we can replace equation (18.4.5) by a quadrature like ci=/summationdisplay µRiµu(xµ)+ni (18.4.7 ) where the N×MmatrixRhas components Riµ≡ri(xµ)(xµ+1−xµ−1)/2( 18.4.8 ) 18.4InverseProblemsandtheUse ofA PrioriInformation 797Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).(or anyothersimple quadrature— it rarelymatters which). We will view equations (18.4.5) and (18.4.7) as being equivalent for practical purposes. How do you solve a set of equations like equation (18.4.7) for the unknown u(xµ)’s? Here is a bad way, but one that contains the germ of some correct ideas: Form a χ2measure of how well a model /hatwideu(x)agrees with the measured data, χ2=N/summationdisplay i=1N/summationdisplay j=1/bracketleftBigg ci−M/summationdisplay µ=1Riµ/hatwideu(xµ)/bracketrightBigg S−1 ij/bracketleftBigg cj−M/summationdisplay µ=1Rjµ/hatwideu(xµ)/bracketrightBigg ≈N/summationdisplay i=1/bracketleftBigg ci−/summationtextM µ=1Riµ/hatwideu(xµ) σi/bracketrightBigg2(18.4.9 ) (compare with equation 15.1.5). Here S−1is the inverse of the covariance matrix, and the approximate equality holds if you can neglect the off-diagonalcovariances, with σi≡(Covar [i, i])1/2. Now you can use the method of singular value decomposition (SVD) in §15.4 to find the vector /hatwideuthat minimizes equation (18.4.9). Don’t try to use the method of normal equations; since Mis greater than Nthey will be singular, as we already discussed. The SVD process will thus surely find a large number of zero singular values,indicativeofahighlynon-uniquesolution. Amongtheinfinityofdegeneratesolutions (most of them badly behaved with arbitrarily large /hatwideu(x µ)’s) SVD will select the one with smallest |/hatwideu|in the sense of /summationdisplay µ[/hatwideu(xµ)]2aminimum (18.4.10 ) (look at Figure 2.6.1). This solution is often called the principal solution .I t is a limiting case of what is called zeroth-order regularization , corresponding to minimizing the sum of the two positive functionals minimize: χ2[/hatwideu]+λ(/hatwideu·/hatwideu)( 18.4.11 ) in the limit of small λ. Below, we will learn how to do such minimizations, as well as more general ones, without the ad hocuse of SVD. Whathappensifwedetermine /hatwideubyequation(18.4.11)withanon-infinitesimal value of λ? First, note that if M/greatermuchN(manymore unknownsthan equations),then uwill often have enough freedom to be able to make χ2(equation 18.4.9) quite unrealisticallysmall, if not zero. Inthe languageof §15.1,the numberof degreesof freedom ν=N−M, which is approximately the expected value of χ2when νis large, is being driven down to zero (and, not meaningfully, beyond). Yet, we know that for the trueunderlying function u(x), which has no adjustable parameters, the numberofdegreesoffreedomandtheexpectedvalueof χ2shouldbeabout ν≈N. Increasing λpullsthesolutionawayfromminimizing χ2infavorofminimizing /hatwideu·/hatwideu. From the preliminarydiscussion above, we can view this as minimizing /hatwideu·/hatwideu subject to the constraint thatχ2have some constant nonzero value. A popular choice,infact,istofindthatvalueof λwhichyields χ2=N,thatis,togetaboutas much extra regularization as a plausible value of χ2dictates. The resulting /hatwideu(x)is calledthe solution of the inverse problem with zeroth-order regularization . 798 Chapter18. IntegralEquationsandInverseTheorySample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).best agreement (independent of smoothness)best smoothness(independent of agreement) best solutions Better SmoothnessBetter Agreementachievable solutions Figure 18.4.1. Almost all inverse problem methods involve a trade-off between two optimizations: agreementbetweendataandsolution,or “sharpness ”ofmappingbetweentrueandestimatedsolution(here denoted A), and smoothness or stability of the solution (here denoted B). Among all possible solutions, shown here schematically as the shaded region, those on the boundary connecting the unconstrainedminimum of Aand the unconstrained minimum of Bare the“best”solutions, in the sense that every other solution is dominated by at least one solution on the curve. The value Nis actually a surrogate for any value drawn from a Gaussian distribution with mean Nand standard deviation (2N)1/2(the asymptotic χ2 distribution). One might equally plausibly try two values of λ, one giving χ2= N+( 2N)1/2, the other N−(2N)1/2. Zeroth-orderregularization,thoughdominatedbybettermethods,demonstrates most of thebasic ideas that areused in inverseproblemtheory. In general,thereare two positive functionals, call them AandB. Thefirst,A, measures something like the agreementof a modelto the data (e.g., χ2), or sometimes a related quantity like the“sharpness ”of the mapping between the solution and the underlying function. When Aby itself is minimized, the agreement or sharpness becomes very good (oftenimpossiblygood),but the solutionbecomes unstable,wildly oscillating,or inother ways unrealistic, re flecting that Aalone typically de fines a highly degenerate minimization problem. That is where Bcomes in. It measures somethinglike the “smoothness ”of the desired solution, or sometimes a related quantity that parametrizes the stability of the solutionwith respect tovariationsin thedata, orsometimes a quantityre flecting a priorijudgments about the likelihood of a solution. Bis called the stabilizing functional orregularizingoperator . Inanycase, minimizing Bby itself is supposed to give a solution that is “smooth”or“stable”or“likely”—and that has nothing at all to do with the measured data. 18.5LinearRegularizationMethods 799Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).The single central idea in inverse theory is the prescription minimize: A+λB (18.4.12 ) for various values of 0<λ< ∞along the so-called trade-off curve (see Figure 18.4.1),andthentosettle ona “best”valueof λbyoneoranothercriterion,ranging from fairly objective (e.g., making χ2=N) to entirely subjective. Successful methods, several of which we will now describe, differ as to their choices of Aand B, as to whether the prescription (18.4.12) yields linear or nonlinear equations, as to their recommendedmethod for selecting a finalλ, and as to their practicality for computer-intensive two-dimensional problems like image processing. They also differ as to the philosophical baggage that they (or rather, their proponents) carry. We have thus far avoided the word “Bayesian. ”(Courts have consistently held that academic license does not extendto shouting “Bayesian”in a crowdedlecture hall.) But it is hard, nor have we any wish, to disguise the fact that Bhas somethingto do with a prioriexpectation,or knowledge,of a solution, while Ahas something to do with a posteriori knowledge. The constant λadjudicates a delicate compromisebetweenthe two. Some inversemethodshaveacquireda more Bayesian stamp than others, but we think that this is purely an accident of history. An outsider looking only at the equations that are actually solved, and not at theaccompanyingphilosophicaljusti fications,wouldhaveadif ficulttimeseparatingthe so-called Bayesian methods from the so-called empirical ones, we think. The next three sections discuss three different approaches to the problem of inversion, which have had considerable success in different fields. All three fit within the general framework that we have outlined, but they are quite different indetail and in implementation. CITED REFERENCES AND FURTHER READING: Craig,I.J.D., andBrown,J.C. 1986, InverseProblemsinAstronomy (Bristol, U.K.: Adam Hilger). Twomey, S. 1977, Introduction to the Mathematics of Inversion in Remote Sensing and Indirect Measurements (Amsterdam: Elsevier). Tikhonov, A.N., and Arsenin, V.Y. 1977, Solutions of Ill-Posed Problems (NewYork: Wiley). Tikhonov, A.N., and Goncharsky, A.V. (eds.) 1987, Ill-Posed Problems in the Natural Sciences (Moscow: MIR). Parker, R.L. 1977, Annual Review of Earth and Planetary Science , vol. 5, pp. 35–64. Frieden, B.R. 1975, in Picture Processing and Digital Filtering , T.S. Huang, ed. (New York: Springer-Verlag). Tarantola, A. 1987, Inverse Problem Theory (Amsterdam: Elsevier). Baumeister,J.1987, StableSolutionofInverseProblems (Braunschweig,Germany:Friedr.Vieweg & Sohn) [mathematically oriented]. Titterington, D.M. 1985, Astronomy and Astrophysics , vol. 144, pp. 381–387. Jeffrey, W., and Rosner, R. 1986, Astrophysical Journal , vol. 310, pp. 463–472. 18.5 Linear Regularization Methods What we will call linear regularization is also called the Phillips-Twomey method[1,2], theconstrained linear inversion method [3], themethod of regulariza- tion[4], andTikhonov-Miller regularization [5-7]. (It probablyhas other names also,