18.6 Backus-Gilbert Method 815 necessary.(For"unsticking"procedures,see [101.)The uniqueness of the solution is also not well understood,although for two-dimensional images of reasonable complexity it is believed to be unique. Deterministic constraints can be incorporated,via projection operators,into iterative methods of linear regularization.In particular,rearranging terms somewhat we can write the iteration (18.5.21)as +1)=(1-eH田）.)+AT.(b-A.) (18.5.27) If the iteration is modified by the insertion of projection operators at each step +1=(P1P2.Pm)[1-AHD·)+AT.(b-A.】 (18.5.28) (or,instead of Pi's,the Ti operators 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. 9 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.Journal of the Association for Computing Machinery.vol.10,pp.97-101.[2] 9 Twomey,S.1977,Introduction to the Mathematics of Inversion in Remote Sensing and Indirect 9 Measurements (Amsterdam:Elsevier).[3] 是 Craig,I.J.D.,and Brown,J.C.1986,Inverse Problems in Astronomy(Bristol,U.K.:Adam Hilger). 4 Tikhonov,A.N.,and Arsenin,V.Y.1977,Solutions of Il-Posed Problems (New York:Wiley).[5] 6 Tikhonov,A.N.,and Goncharsky,A.V.(eds.)1987,IlI-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 /EEE.vol.69. Pp.432-450. Biemond,J.,Lagendijk,R.L.,and Mersereau,R.M.1990,Proceedings of the /EEE,vol.78. pp.856-883.7 Numerica 10621 Gerchberg,R.W.and Saxton,W.O.1972.Optik,vol.35.pp.237-246.[8] 431 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, (outside Recipes pp.1897-1907.[101 North Software. 18.6 Backus-Gilbert Method The Backus-Gilbert method [1.2](see,e.g.,[3]or [4]for summaries)differs from other regularization methods in the nature of its functionals A and B.For B,the method seeks to maximize the stability of the solution u(x)rather than,in the first instance,its smoothness.That is, B≡Var[a(x] (18.6.1)
18.6 Backus-Gilbert Method 815 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). necessary. (For “unsticking” procedures, see [10].) The uniqueness of the solution is also not well understood, although for two-dimensional images of reasonable complexity it is believed to be unique. Deterministic constraints can be incorporated, via projection operators, into iterative methods of linear regularization. In particular, rearranging terms somewhat, we can write the iteration (18.5.21) as u(k+1) = (1 − λH) · u(k) + AT · (b − A · u(k) ) (18.5.27) If the iteration is modified by the insertion of projection operators at each step u(k+1) = (P1P2 ···Pm)[(1 − λH) · u(k) + AT · (b − A · u(k) )] (18.5.28) (or, instead of Pi’s, the Ti operators 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, Journal of 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., and Brown, J.C. 1986, Inverse Problems in Astronomy (Bristol, U.K.: Adam Hilger). [4] Tikhonov, A.N., and Arsenin, V.Y. 1977, Solutions of Ill-Posed Problems (New York: 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 The Backus-Gilbert method [1,2] (see, e.g., [3] or [4] for summaries) differs from other regularization methods in the nature of its functionals A and B. For B, the method seeks to maximize the stability of the solution u(x) rather than, in the first instance, its smoothness. That is, B ≡ Var[u(x)] (18.6.1)
816 Chapter 18.Integral Equations and Inverse Theory is used as a measure of how much the solution (z)varies as the data vary within their measurement errors.Note that this variance is not the expected deviation of (z)from the true u(x)-that will be constrained by A-but rather measures the expected experiment-to-experiment scatter among estimates u(z)if the whole experiment were to be repeated many times. For A the Backus-Gilbert method looks at the relationship between the solution ()and the true function u(),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 (x)and u(x)can be written as 81 (x)= 5(x,r')u(r)dr (18.6.2) for some so-called resolution function or averaging kernel 6(r,').The Backus- Gilbert method seeks to minimize the width or spread of (that is,maximize the resolving power).A is chosen to be some positive measure of the spread. While Backus-Gilbert's philosophy is thus rather different from that of Phillips- Twomey and related methods,in practice the differences between the methods are less than one might think.A stable solution is almost inevitably bound to be smooth:The wild.unstable oscillations that result from an unregularized solution g2 are always exquisitely sensitive to small changes in the data.Likewise,making (z)close to u(x)inevitably will bring error-free data into agreement with the model.Thus A and B play roles closely analogous to their corresponding roles 豆芝② 9 in the previous two sections The principal advantage of the Backus-Gilbert formulation is that it gives good control over just those properties that it seeks to measure,namely stability and resolving power.Moreover,in the Backus-Gilbert method,the choice of A(playing 61 its usual role of compromise between A and B)is conventionally made,or at least can easily be made,before any actual data are processed.One's uneasiness at making a post hoc,and therefore potentially subjectively biased,choice of A is thus removed. 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), Numerica 10621 431 ci=si+ni=ri()u()dx+ni (18.6.3) Recipes and building in linearity from the start,we seek a set of inverse response kernels qi()such that )=∑ga(z)c (18.6.4) is the desired estimator of u(z).It is useful to define the integrals of the response kernels for each data point, R:三T()da (18.6.5)
816 Chapter 18. Integral Equations and Inverse Theory Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). is used as a measure of how much the solution u(x) varies as the data vary within their measurement errors. Note that this variance is not the expected deviation of u(x) from the true u(x) — that will be constrained by A — but rather measures the expected experiment-to-experiment scatter among estimates u(x) if the whole experiment were to be repeated many times. For A the Backus-Gilbert method looks at the relationship between the solution u(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 u(x) and u(x) can be written as u(x) = δ (x, x )u(x )dx (18.6.2) for some so-called resolution function or averaging kernel δ (x, x ). The BackusGilbert method seeks to minimize the width or spread of δ (that is, maximize the resolving power). A is chosen to be some positive measure of the spread. While Backus-Gilbert’s philosophy is thus rather different from that of PhillipsTwomey and related methods, in practice the differences between the methods are less than one might think. A stable solution 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 u(x) close to u(x) inevitably will bring error-free data into agreement with the model. Thus A and B play roles closely analogous to their corresponding roles in the previous two sections. The principal advantage of the Backus-Gilbert formulation is that it gives good control over just those properties that it seeks to measure, namely stability and resolving power. Moreover, in the Backus-Gilbert method, the choice of λ (playing its usual role of compromise between A and B) is conventionally made, or at least can easily be made, before any actual data are processed. One’s uneasiness at making a post hoc, and therefore potentially subjectively biased, choice of λ is thus removed. 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 = 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 u(x) = i qi(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 ≡ ri(x)dx (18.6.5)
18.6 Backus-Gilbert Method 817 Substituting equation(18.6.4)into equation(18.6.3),and comparing with equation (18.6.2),we see that x,x)=∑g(ar(x) (18.6.6) We can require this averaging kernel to have unit area at every z,giving 1=,r=∑9a回r(dr=∑oR三qR(I86.7) where g(z)and R are each vectors of length N,the number of measurements. Standard propagation of errors,and equation (18.6.1),give B=Var(z=∑∑9(o)Sg(a)=q(a)·s·q(a) (18.6.8) ij where Si;is the covariance matrix (equation 18.4.6).If one can neglect off-diagonal server 2 covariances (as when the errors on the ci's are independent),then Sj=j? (Nort University is diagonal. Press. THE We now need to define a measure of the width or spread of 6(,')at each value ofr.While many choices are possible,Backus and Gilbert choose the second moment of its square.This measure becomes the functional A, Program A三w(x)= (d'-x)26,x']2dr OF SCIENTIFIC( (18.6.9) =∑∑9()Ww(r()三q()w(·q) 6 where we have here used equation(18.6.6)and defined the spread matrix W()by W()三(x'-x)2rn(rr(x'）dr (18.6.10) Numerical 10621 The functions gi(z)are now determined by the minimization principle idge.org 431 Recipes minimize:A+λB=q(x)·W(x)+λS·q(x) (18.6.11) (outside subject to the constraint (18.6.7)that q(x).R =1. North The solution of equation (18.6.11)is [W(z)+λS-1.R q)=R·w国)+S-1R (18.6.12) (Reference [4]gives an accessible proof.)For any particular data set c (set of measurements ci),the solution (z)is thus aa-图+SR (18.6.13)
18.6 Backus-Gilbert Method 817 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). Substituting equation (18.6.4) into equation (18.6.3), and comparing with equation (18.6.2), we see that δ (x, x ) = i qi(x)ri(x ) (18.6.6) We can require this averaging kernel to have unit area at every x, giving 1 = δ (x, x )dx = i qi(x) ri(x )dx = i qi(x)Ri ≡ q(x) · R (18.6.7) where q(x) and R are each vectors of length N, the number of measurements. Standard propagation of errors, and equation (18.6.1), give B = Var[u(x)] = i j qi(x)Sij qj (x) = q(x) · S · q(x) (18.6.8) where Sij is the covariance matrix (equation 18.4.6). If one can neglect off-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 δ (x, x ) at each value of x. While many choices are possible, Backus and Gilbert choose the second moment of its square. This measure becomes the functional A, A ≡ w(x) = (x − x) 2[δ (x, x )]2dx = i j qi(x)Wij (x)qj (x) ≡ q(x) · W(x) · q(x) (18.6.9) where we have here used equation (18.6.6) and defined the spread matrix W(x) by Wij (x) ≡ (x − x) 2ri(x )rj (x )dx (18.6.10) The functions qi(x) are now determined by the minimization principle minimize: A + λB = q(x) · W(x) + λS · 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 u(x) is thus u(x) = c · [W(x) + λS] −1 · R R · [W(x) + λS] −1 · R (18.6.13)
818 Chapter 18.Integral Equations and Inverse Theory (Don't let this notation mislead you into inverting the full matrix W(z)+AS.You only need to solve for some y the linear system (W(z)+AS).y =R,and then substitute y into both the numerators and denominators of 18.6.12 or 18.6.13.) Equations(18.6.12)and (18.6.13)have a completely different character from the linearly regularized solutions to(18.5.7)and(18.5.8).The vectors and matrices in (18.6.12)all have size N,the number of measurements.There is no discretization of the underlying variable z,so M does not come into play at all.One solves a different N x N set of linear equations for each desired value of z.By contrast,in(18.5.8), one solves an M x M linear set,but only once.In general,the computational burden of repeatedly solving linear systems makes the Backus-Gilbert method unsuitable 81 for other than one-dimensional problems. How does one choose A within the Backus-Gilbert scheme?As already mentioned,you can (in some cases should)make the choice before you see any actual data.For a given trial value of A.and for a sequence of x's,use equation(18.6.12) to calculate q();then use equation(18.6.6)to plot the resolution functions 6(,' as a function of x'.These plots will exhibit the amplitude with which different underlying values z'contribute to the point ()of your estimate.For the same value of A,also plot the function Var()]using equation (18.6.8).(You need an estimate of your measurement covariance matrix for this. As you change A you will see very explicitly the trade-off between resolution and stability.Pick the value that meets your needs.You can even choose A to be a function of A =A(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 z.)For the chosen value or values of A,you now have a quantitative understanding of your 三兰号∽6 inverse solution procedure.This can prove invaluable if-once you are processing OF SCIENTIFIC 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 found particular success among geophysicists,who use it to obtain information about the structure of the Earth(e.g.,density run with depth)from seismic travel time data 、复公 CITED REFERENCES AND FURTHER READING: Backus,G.E.,and Gilbert,F.1968,Geophysical Journal of the Royal Astronomical Society, vol.16,pp.169-205.I1) Numerica 10621 Backus,G.E..and Gilbert,F.1970,Philosophical Transactions of the Royal Society of London 431 A,vol.266,pp.123-192.2 Recipes 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] (outside Software. 18.7 Maximum Entropy Image Restoration Above.we commented that the association of certain inversion methodsbreak 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
818 Chapter 18. Integral Equations and Inverse Theory Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). (Don’t let this notation mislead you into inverting the full matrix W(x) + λS. You only need to solve for some y the linear system (W(x) + λS) · y = R, and then substitute y into both the numerators and denominators of 18.6.12 or 18.6.13.) Equations (18.6.12) and (18.6.13) have a completely different character from the linearly regularized solutions to (18.5.7) and (18.5.8). The vectors and matrices in (18.6.12) all have size N, the number of measurements. There is no discretization of the underlying variable x, so M does not come into play at all. One solves a different N × N set of linear equations for each desired value of x. By contrast, in (18.5.8), one solves an M ×M linear set, but only once. In general, the computational burden of repeatedly solving linear systems makes the Backus-Gilbert method unsuitable for other than one-dimensional problems. How does one choose λ within the Backus-Gilbert scheme? As already mentioned, you can (in some cases should) make the choice before you see any actual data. For a given trial value of λ, and for a sequence of x’s, use equation (18.6.12) to calculate q(x); then use equation (18.6.6) to plot the resolution functions δ (x, x ) as a function of x . These plots will exhibit the amplitude with which different underlying values x contribute to the point u(x) of your estimate. For the same value of λ, also plot the function Var[u(x)] using equation (18.6.8). (You need an 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 understanding of your inverse solution procedure. This can prove invaluable if — once you are 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 found particular success among geophysicists, who use it to obtain information about the structure of the Earth (e.g., density run with depth) from seismic travel time data. 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 methodsbreak 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