f18-6

806 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).necessary. (For “unsticking” procedures, see [10].) The uniqueness of the solution is also not well understood, although for two-dimensional images of reasonablecomplexity it is believed to be unique. Deterministic constraints can be incorporated, via projection operators, into iterativemethodsoflinearregularization. Inparticular,rearrangingtermssomewhat,we can write the iteration (18.5.21) as /hatwideu (k+1)=(1−/epsilon1λH)·/hatwideu(k)+/epsilon1AT·(b−A·/hatwideu(k))( 18.5.27 ) If the iteration is modified by the insertion of projectionoperators at each step /hatwideu(k+1)=(P1P2···P m)[(1−/epsilon1λH)·/hatwideu(k)+/epsilon1AT·(b−A·/hatwideu(k))] (18.5.28 ) (or, instead of Pi’s, the Tioperators of equation 18.5.26),then it can be shown that the convergence condition (18.5.22) is unmodified, and the iteration will converge to minimize the quadratic functional (18.5.6) subject to the desired nonlinear deterministic constraints. See [7]for references to more sophisticated, and faster converging, iterations along these lines. CITED REFERENCES AND FURTHER READING: Phillips, D.L. 1962, Journal of the Association for Computing Machinery , vol. 9, pp. 84–97. [1] Twomey, S. 1963, Journalof the Association for Computing Machinery , vol. 10, pp. 97–101. [2] Twomey, S. 1977, Introduction to the Mathematics of Inversion in Remote Sensing and Indirect Measurements (Amsterdam: Elsevier). [3] Craig,I.J.D., andBrown,J.C. 1986, InverseProblemsinAstronomy (Bristol, U.K.: Adam Hilger). [4] Tikhonov,A.N., andArsenin, V.Y. 1977, Solutions of Ill-Posed Problems (NewYork: Wiley). [5] Tikhonov, A.N., and Goncharsky, A.V. (eds.) 1987, Ill-Posed Problems in the Natural Sciences (Moscow: MIR). Miller, K. 1970, SIAM Journal on Mathematical Analysis , vol. 1, pp. 52–74. [6] Schafer, R.W., Mersereau, R.M., and Richards, M.A. 1981, Proceedings of the IEEE , vol. 69, pp. 432–450. Biemond, J., Lagendijk, R.L., and Mersereau, R.M. 1990, Proceedings of the IEEE , vol. 78, pp. 856–883. [7] Gerchberg, R.W., and Saxton, W.O. 1972, Optik, vol. 35, pp. 237–246. [8] Fienup, J.R. 1982, Applied Optics , vol. 15, pp. 2758–2769. [9] Fienup, J.R., and Wackerman, C.C. 1986, Journal of the Optical Society of America A , vol. 3, pp. 1897–1907. [10] 18.6 Backus-Gilbert Method TheBackus-Gilbertmethod [1,2](see,e.g., [3]or[4]forsummaries)differsfrom other regularization methods in the nature of its functionals AandB.F o r B, the method seeks to maximize the stabilityof the solution /hatwideu(x)rather than, in the first instance, its smoothness. That is, B≡Var [/hatwideu(x)] ( 18.6.1 ) 18.6Backus-GilbertMethod 807Sample 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).is used as a measure of how much the solution /hatwideu(x)varies as the data vary within their measurement errors. Note that this variance is not the expected deviation of /hatwideu(x)from the true u(x)— that will be constrained by A— but rather measures the expected experiment-to-experiment scatter among estimates /hatwideu(x)if the whole experiment were to be repeated many times. ForAtheBackus-Gilbertmethodlooksattherelationshipbetweenthesolution /hatwideu(x)and the true function u(x), and seeks to make the mapping between these as close to the identity map as possible in the limit of error-free data. The method is linear, so the relationship between /hatwideu(x)andu(x)can be written as /hatwideu(x)=/integraldisplay /hatwideδ(x, x/prime)u(x/prime)dx/prime(18.6.2 ) for some so-called resolution function oraveraging kernel /hatwideδ(x, x/prime). The Backus- Gilbert method seeks to minimize the width or spreadof/hatwideδ(that is, maximize the resolving power). Ais chosen to be some positive measure of the spread. WhileBackus-Gilbert’sphilosophyisthusratherdifferentfromthatofPhillips- Twomey and related methods, in practice the differences between the methods are less than one might think. A stablesolution is almost inevitably bound to be smooth: The wild, unstable oscillations that result from an unregularized solution are always exquisitely sensitive to small changes in the data. Likewise, making /hatwideu(x)close to u(x)inevitably will bring error-free data into agreement with the model. Thus AandBplay roles closely analogous to their corresponding roles in the previous two sections. TheprincipaladvantageoftheBackus-Gilbertformulationis thatit givesgood control over just those properties that it seeks to measure, namely stability and resolving power. Moreover,in the Backus-Gilbertmethod, the choice of λ(playing its usual role of compromise between AandB) is conventionally made, or at least caneasilybemade, beforeanyactualdataareprocessed. One’suneasinessatmaking aposthoc,andthereforepotentiallysubjectivelybiased,choiceof λisthusremoved. Backus-Gilbert is often recommended as the method of choice for designing, and predicting the performance of, experiments that require data inversion. Let’s see how this all works. Starting with equation (18.4.5), ci≡si+ni=/integraldisplay ri(x)u(x)dx +ni (18.6.3 ) and building in linearity from the start, we seek a set of inverse response kernels qi(x)such that /hatwideu(x)=/summationdisplay iqi(x)ci (18.6.4 ) is the desired estimator of u(x). It is useful to define the integrals of the response kernels for each data point, Ri≡/integraldisplay ri(x)dx (18.6.5 ) 808 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).Substituting equation (18.6.4) into equation (18.6.3), and comparing with equation (18.6.2), we see that /hatwideδ(x, x/prime)=/summationdisplay iqi(x)ri(x/prime)( 18.6.6 ) We can require this averaging kernel to have unit area at every x, giving 1=/integraldisplay /hatwideδ(x, x/prime)dx/prime=/summationdisplay iqi(x)/integraldisplay ri(x/prime)dx/prime=/summationdisplay iqi(x)Ri≡q(x)·R (18.6.7 ) whereq(x)andRare each vectors of length N, the numberof measurements. Standard propagation of errors, and equation (18.6.1), give B=Var [/hatwideu(x)] =/summationdisplay i/summationdisplay jqi(x)Sijqj(x)=q(x)·S·q(x)( 18.6.8 ) where Sijisthecovariancematrix(equation18.4.6). Ifonecanneglectoff-diagonal covariances (as when the errors on the ci’s are independent), then Sij=δijσ2 i is diagonal. We now need to define a measure of the width or spread of /hatwideδ(x, x/prime)at each valueof x. While manychoicesarepossible,Backus andGilbertchoosethesecond moment of its square. This measure becomes the functional A, A≡w(x)=/integraldisplay (x/prime−x)2[/hatwideδ(x, x/prime)]2dx/prime =/summationdisplay i/summationdisplay jqi(x)Wij(x)qj(x)≡q(x)·W(x)·q(x)(18.6.9 ) wherewe havehereusedequation(18.6.6)anddefinedthe spreadmatrix W(x)by Wij(x)≡/integraldisplay (x/prime−x)2ri(x/prime)rj(x/prime)dx/prime(18.6.10 ) The functions qi(x)are now determinedby the minimization principle minimize: A+λB=q(x)·/bracketleftbig W(x)+λS/bracketrightbig ·q(x)( 18.6.11 ) subject to the constraint (18.6.7) that q(x)·R=1. The solution of equation (18.6.11) is q(x)=[W(x)+λS]−1·R R·[W(x)+λS]−1·R(18.6.12 ) (Reference [4]gives an accessible proof.) For any particular data set c(set of measurements ci), the solution /hatwideu(x)is thus /hatwideu(x)=c·[W(x)+λS]−1·R R·[W(x)+λS]−1·R(18.6.13 ) 18.7Maximum EntropyImageRestoration 809Sample 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).(Don’t let this notation mislead you into inverting the full matrix W(x)+λS.Y o u only need to solve for some ythe linear system (W(x)+λS)·y=R, and then substitute yinto both the numeratorsand denominatorsof 18.6.12or 18.6.13.) Equations (18.6.12) and (18.6.13) have a completely different character from thelinearlyregularizedsolutionsto(18.5.7)and(18.5.8). Thevectorsandmatricesin(18.6.12)allhavesize N,thenumberofmeasurements. Thereisnodiscretizationof theunderlyingvariable x,soMdoesnotcomeintoplayatall. Onesolvesadifferent N×Nset of linear equations for each desired value of x. By contrast, in (18.5.8), onesolvesan M×Mlinearset,butonlyonce. Ingeneral,thecomputationalburden of repeatedly solving linear systems makes the Backus-Gilbert method unsuitablefor other than one-dimensional problems. How does one choose λwithin the Backus-Gilbert scheme? As already mentioned,youcan(insomecases should)makethechoice beforeyouseeanyactual data. For a given trial value of λ, and for a sequence of x’s, use equation (18.6.12) tocalculate q(x); thenuse equation(18.6.6)to plottheresolutionfunctions /hatwideδ(x, x /prime) as a function of x/prime. These plots will exhibit the amplitude with which different underlying values x/primecontribute to the point /hatwideu(x)of your estimate. For the same valueof λ, also plotthefunction/radicalbig Var [/hatwideu(x)]usingequation(18.6.8). (Youneedan estimate of your measurement covariance matrix for this.) As you change λyou will see very explicitly the trade-off between resolution and stability. Pick the value that meets your needs. You can even choose λto be a function of x,λ=λ(x), in equations (18.6.12)and (18.6.13), should you desire to do so. (This is one benefit of solving a separate set of equations for each x.) For the chosen value or values of λ, you now have a quantitative understandingof your inversesolutionprocedure. This can proveinvaluableif — once youare processing real data — you need to judge whether a particular feature, a spike or jump for example, is genuine, and/or is actually resolved. The Backus-Gilbert method has foundparticularsuccessamonggeophysicists,whouseittoobtaininformationaboutthestructureoftheEarth(e.g.,densityrunwithdepth)fromseismictraveltimedata. CITED REFERENCES AND FURTHER READING: Backus, G.E., and Gilbert, F. 1968, Geophysical Journal of the Royal Astronomical Society , vol. 16, pp. 169–205. [1] Backus, G.E., and Gilbert, F. 1970, Philosophical Transactions of the Royal Society of London A, vol. 266, pp. 123–192. [2] Parker, R.L. 1977, Annual Review of Earth and Planetary Science , vol. 5, pp. 35–64. [3] Loredo, T.J., and Epstein, R.I. 1989, Astrophysical Journal , vol. 336, pp. 896–919. [4] 18.7 Maximum Entropy Image Restoration Above, we commented that the association of certain inversion methods with Bayesian arguments is more historical accident than intellectual imperative. Maximum entropy methods , so-called, are notorious in this regard; to summarize these methods without some, at least introductory, Bayesian invocations would be to serve a steak without the sizzle, or a sundae without the cherry. We should