f18-5

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 theaccompanyingphilosophicaljustifications,wouldhaveadifficulttimeseparatingthe 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, 800 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).sinceitissoobviouslyagoodidea.) Initssimplestform,themethodisanimmediate generalization of zeroth-order regularization (equation 18.4.11, above). As before,thefunctional Aistakentobethe χ 2deviation,equation(18.4.9),butthefunctional Bis replaced by more sophisticated measures of smoothness that derive from first or higher derivatives. For example,supposethat your a prioribelief is that a credible u(x)is not too different from a constant. Then a reasonable functional to minimize is B∝/integraldisplay [/hatwideu/prime(x)]2dx∝M−1/summationdisplay µ=1[/hatwideuµ−/hatwideuµ+1]2(18.5.1 ) since it is nonnegative and equal to zero only when /hatwideu(x)is constant. Here /hatwideuµ≡/hatwideu(xµ), and the second equality (proportionality) assumes that the xµ’s are uniformly spaced. We can write the second form of Bas B=|B·/hatwideu|2=/hatwideu·(BT·B)·/hatwideu≡/hatwideu·H·/hatwideu (18.5.2 ) where/hatwideuis the vector of components /hatwideuµ,µ =1,...,M,Bis the (M−1)×M first difference matrix B= −1100000 ··· 0 0−110000 ··· 0 ......... 0··· 0000 −110 0··· 00000 −11 (18.5.3 ) andHis theM×Mmatrix H=B T·B= 1−100000 ··· 0 −12 −10000 ··· 0 0−12 −1000 ··· 0 ......... 0··· 000 −12 −10 0··· 0000 −12 −1 0··· 00000 −11 (18.5.4 ) Note that Bhas one fewer row than column. It follows that the symmetric H is degenerate; it has exactly one zero eigenvalue corresponding to the valueof a constant function, any one of which makes Bexactly zero. If, just as in §15.4, we write A iµ≡Riµ/σibi≡ci/σi (18.5.5 ) then, using equation (18.4.9), the minimization principle (18.4.12) is minimize: A+λB=|A·/hatwideu−b|2+λ/hatwideu·H·/hatwideu (18.5.6 ) Thiscanreadilybereducedtoalinearsetof normalequations ,just asin §15.4: The components /hatwideuµof the solution satisfy the set of Mequations in Munknowns, /summationdisplay ρ/bracketleftBigg/parenleftbigg/summationdisplay iAiµAiρ/parenrightbigg +λHµρ/bracketrightBigg /hatwideuρ=/summationdisplay iAiµbiµ=1,2,...,M (18.5.7 ) 18.5LinearRegularizationMethods 801Sample 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, in vector notation, (AT·A+λH)·/hatwideu=AT·b (18.5.8 ) Equations (18.5.7) or (18.5.8) can be solved by the standard techniques of Chapter 2, e.g., LUdecomposition. The usual warnings about normal equations being ill-conditioned do not apply, since the whole purpose of the λterm is to cure thatsameill-conditioning. Note,however,thatthe λtermbyitselfisill-conditioned, since it does not select a preferred constant value. You hope your data can at least do that! Althoughinversionofthematrix (AT·A+λH)isnotgenerallythebestwayto solvefor /hatwideu, let usdigress towrite thesolutiontoequation(18.5.8)schematicallyas /hatwideu=/parenleftbigg1 AT·A+λH·AT·A/parenrightbigg A−1·b(schematiconly!) (18.5.9 ) where the identity matrix in the form A·A−1has been inserted. This is schematic not only because the matrix inverse is fancifully written as a denominator, butalso because, in general, the inverse matrix A −1does not exist. However, it is illuminating to compare equation (18.5.9) with equation (13.3.6) for optimal or Wiener filtering, or with equation (13.6.6) for general linear prediction. One sees thatAT·Aplays the role of S2, the signal power or autocorrelation, while λH plays the role of N2, the noise power or autocorrelation. The term in parentheses in equation (18.5.9) is something like an optimal filter, whose effect is to pass the ill-posed inverse A−1·bthrough unmodified when AT·Ais sufficiently large, but to suppress it when AT·Ais small. The above choices of BandHare only the simplest in an obvioussequenceof derivatives. If your a prioribelief is that a linearfunction is a good approximation tou(x), then minimize B∝/integraldisplay [/hatwideu/prime/prime(x)]2dx∝M−2/summationdisplay µ=1[−/hatwideuµ+2/hatwideuµ+1−/hatwideuµ+2]2(18.5.10 ) implying B= −12 −10000 ··· 0 0−12 −1000 ··· 0 ......... 0··· 000 −12 −10 0··· 0000 −12 −1 (18.5.11 ) and H=B T·B= 1−210000 ··· 0 −25 −41000 ··· 0 1−46 −4100 ··· 0 01 −46 −410 ··· 0 ......... 0··· 01 −46 −410 0··· 001 −46 −41 0··· 0001 −45 −2 0··· 00001 −21 (18.5.12 ) 802 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).ThisHhastwozeroeigenvalues,correspondingtothetwoundeterminedparameters of a linear function. Ifyouraprioribeliefis that a quadratic functionis preferable,thenminimize B∝/integraldisplay [/hatwideu/prime/prime/prime(x)]2dx∝M−3/summationdisplay µ=1[−/hatwideuµ+3/hatwideuµ+1−3/hatwideuµ+2+/hatwideuµ+3]2(18.5.13 ) with B= −13 −31000 ··· 0 0−13 −3100 ··· 0 ......... 0··· 00 −13 −310 0··· 000 −13 −31 (18.5.14 ) and now H= 1−33 −100000 ··· 0 −31 0 −12 6 −10000 ··· 0 3−12 19 −15 6 −1000 ··· 0 −16 −15 20 −15 6 −100 ··· 0 0−16 −15 20 −15 6 −10 ··· 0 ......... 0··· 0−16 −15 20 −15 6 −10 0··· 00 −16 −15 20 −15 6 −1 0··· 000 −16 −15 19 −12 3 0··· 0000 −16 −12 10 −3 0··· 00000 −13 −31  (18.5.15 ) (We’ll leave the calculation of cubics and above to the compulsive reader.) Notice that you can regularize with “closeness to a differential equation,” if you want. Just pick Bto be the appropriate sum of finite-difference operators (the coefficients can depend on x), and calculate H=B T·B. You don’t need to know the values of yourboundaryconditions,since Bcan have fewer rows than columns, asabove;hopefully,yourdatawilldeterminethem. Ofcourse,ifyoudoknowsomeboundary conditions, you can build these into Btoo. Withalltheproportionalitysignsabove,youmayhavelosttrackofwhatactual valueofλto try first. A simple trickforat least getting“onthe map”is tofirst try λ=Tr(A T·A)/Tr(H)( 18.5.16 ) where Tr is the trace of the matrix (sum of diagonal components). This choice will tend to make the two parts of the minimization have comparable weights, and you can adjust from there. As for what is the “correct” value of λ, an objective criterion, if you know your errors σiwith reasonable accuracy, is to make χ2(that is, |A·/hatwideu−b|2) equal toN, the number of measurements. We remarked above on the twin acceptable choicesN±(2N)1/2. Asubjectivecriterionistopickanyvaluethatyoulikeinthe 18.5LinearRegularizationMethods 803Sample 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).range 0<λ< ∞, dependingon your relative degreeof belief in the a priorianda posteriori evidence. (Yes, people actually do that. Don’t blame us.) Two-DimensionalProblemsand Iterative Methods Up to now our notation has been indicative of a one-dimensional problem, finding /hatwideu(x)or/hatwideuµ=/hatwideu(xµ). However,allofthediscussioneasilygeneralizestothe problemofestimating a two-dimensionalset ofunknowns /hatwideuµκ,µ=1,...,M,κ = 1,...,K, corresponding, say, to the pixel intensities of a measured image. In this case, equation (18.5.8) is still the one we want to solve. In image processing, it is usual to have the same number of input pixels in a measured “raw” or “dirty” image as desired “clean” pixels in the processed output image,sothematrices RandA(equation18.5.5)aresquareandofsize MK×MK. Ais typically much too large to represent as a full matrix, but often it is either (i) sparse, with coefficients blurring an underlying pixel (i,j)only into measurements (i±few,j±few ),or(ii)translationallyinvariant,sothat A(i,j)(µ,ν )=A(i−µ,j−ν). Both of these situations lead to tractable problems. In the case of translational invariance, fast Fourier transforms (FFTs) are the obvious method of choice. The general linear relation between underlyingfunction and measured values (18.4.7) now becomes a discrete convolution like equation(13.1.1). If kdenotesatwo-dimensionalwave-vector,thenthetwo-dimensionalFFT takes us back and forth between the transform pairs A(i−µ,j−ν)⇐⇒/tildewideA(k)b (i,j)⇐⇒/tildewideb(k)/hatwideu(i,j)⇐⇒/tildewideu(k)(18.5.17 ) We also needa regularizationorsmoothingoperator Bandthederived H=BT·B. One popular choice for Bis the five-point finite-difference approximation of the Laplacian operator, that is, the difference between the value of each point and the average of its fourCartesian neighbors. In Fourier space, this choice implies, /tildewideB(k)∝sin2(πk 1/M )s i n2(πk 2/K ) /tildewideH(k)∝sin4(πk 1/M )s i n4(πk 2/K )(18.5.18 ) In Fourier space, equation (18.5.7) is merely algebraic, with solution /tildewideu(k)=/tildewideA*(k)/tildewideb(k) |/tildewideA(k)|2+λ/tildewideH(k)(18.5.19 ) whereasterisk denotescomplexconjugation. Youcanmakeuse ofthe FFT routines for real data in §12.5. Turn now to the case where Ais not translationally invariant. Direct solution of (18.5.8) is now hopeless, since the matrix Ais just too large. We need some kind of iterative scheme. One way to proceed is to use the full machinery of the conjugate gradient method in §10.6 to find the minimumof A+λB, equation (18.5.6). Of the various methods in Chapter 10, conjugate gradient is the unique best choice because (i) it does not require storage of a Hessian matrix, which would be infeasible here, 804 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).and (ii) it does exploit gradient information, which we can readily compute: The gradient of equation (18.5.6) is ∇(A+λB)=2 [ (AT·A+λH)·/hatwideu−AT·b]( 18.5.20 ) (cf. 18.5.8). Evaluation of both the function and the gradient should of course take advantage of the sparsity of A, for example via the routines sprsaxandsprstx in§2.7. We will discuss the conjugate gradient technique further in §18.7, in the context of the (nonlinear) maximum entropy method. Some of that discussion can apply here as well. The conjugate gradient method notwithstanding, application of the unsophis- ticated steepest descent method (see §10.6) can sometimes produce useful results, particularly when combined with projections onto convex sets (see below). If the solution after kiterations is denoted /hatwideu(k), then after k+1iterations we have /hatwideu(k+1)=[1−/epsilon1(AT·A+λH)]·/hatwideu(k)+/epsilon1AT·b (18.5.21 ) Here/epsilon1isaparameterthatdictateshowfartomoveinthedownhillgradientdirection. The method converges when /epsilon1is small enough, in particular satisfying 0</epsilon1<2 maxeigenvalue (AT·A+λH)(18.5.22 ) There exist complicated schemes for finding optimal values or sequences for /epsilon1, see[7]; or, one can adopt an experimental approach, evaluating (18.5.6) to be sure that downhill steps are in fact being taken. In those image processing problems where the final measure of success is somewhat subjective (e.g., “how good does the picture look?”), iteration (18.5.21) sometimes produces significantly improved images long before convergence is achieved. This probably accounts for much of its use, since its mathematicalconvergenceis extremelyslow. In fact, (18.5.21)can be used with H=0, in which case the solution is not regularizedat all, and full convergencewouldbe disastrous! This is called Van Cittert’s method and goes back to the 1930s. A number of iterations the order of 1000 is not uncommon [7]. Deterministic Constraints: Projections ontoConvex Sets A set of possible underlying functions (or images) {/hatwideu}is said to be convexif, foranytwo elements /hatwideuaand/hatwideubin theset, all thelinearlyinterpolatedcombinations (1−η)/hatwideua+η/hatwideub 0≤η≤1( 18.5.23 ) arealsointheset. Many deterministicconstraints thatonemightwanttoimposeon the solution /hatwideuto an inverse problem in fact define convexsets, for example: •positivity •compact support (i.e., zero value outside of a certain region) 18.5LinearRegularizationMethods 805Sample 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).•known bounds (i.e., uL(x)≤/hatwideu(x)≤uU(x)for specified functions uL anduU). (Inthislastcase,theboundsmightberelatedtoaninitialestimateandits errorbars, e.g.,/hatwideu0(x)±γσ(x), whereγis of order 1 or 2.) Notice that these, and similar, constraintscanbeeitherintheimagespace,orintheFouriertransformspace,or(infact) in the space of any linear transformation of /hatwideu. IfC iis a convexset, then Piis called a nonexpansiveprojectionoperator onto thatsetif(i) Pileavesunchangedany /hatwideualreadyin Ci,and(ii) Pimapsany /hatwideuoutside Cito theclosestelement of Ci, in the sense that |Pi/hatwideu−/hatwideu|≤|/hatwideua−/hatwideu|forall/hatwideuainCi (18.5.24 ) Whilethisdefinitionsoundscomplicated,examplesareverysimple: Anonexpansive projection onto the set of positive /hatwideu’s is “set all negative components of /hatwideuequal to zero.” A nonexpansive projection onto the set of /hatwideu(x)’s bounded by uL(x)≤ /hatwideu(x)≤uU(x)is “set all values less than the lower bound equal to that bound, and set all values greater than the upper bound equal to thatbound.” A nonexpansive projection onto functions with compact support is “zero the values outside of the region of support.” Theusefulnessofthesedefinitionsis thefollowingremarkabletheorem: Let C be the intersection of mconvex sets C1,C 2,...,C m. Then the iteration /hatwideu(k+1)=(P1P2···P m)/hatwideu(k)(18.5.25 ) will converge to Cfrom all starting points, as k→∞. Also, if Cis empty (there is no intersection), then the iteration will have no limit point. Application of this theoremiscalledthe methodofprojectionsontoconvexsets orsometimes POCS[7]. A generalization of the POCS theorem is that the Pi’s can be replaced by a set of Ti’s, Ti≡1+βi(Pi−1)0<βi<2( 18.5.26 ) A well-chosenset of βi’s can acceleratethe convergencetothe intersectionset C. Some inverse problems can be completely solved by iteration (18.5.25)alone! For example, a problem that occurs in both astronomical imaging and X-ray diffraction work is to recover an image given only the modulus of its Fourier transform (equivalent to its power spectrum or autocorrelation) and not the phase. Herethreeconvexsetscanbeutilized: thesetofallimageswhoseFouriertransform has the specified modulus to within specified error bounds; the set of all positive images;andthesetofallimageswithzerointensityoutsideofsomespecifiedregion. In this case the POCS iteration (18.5.25) cycles among these three, imposing each constraint in turn; FFTs are used to get in and out of Fourier space each time theFourier constraint is imposed. The specific application of POCS to constraints alternately in the spatial and Fourier domains is also known as the Gerchberg-Saxton algorithm [8]. While this algorithmis non-expansive,and is frequentlyconvergentin practice, it has not been provedto convergein all cases [9]. In the phase-retrievalproblemmentionedabove, the algorithm often “gets stuck” on a plateau for many iterations before making sudden, dramatic improvements. As many as 104to105iterations are sometimes 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 rB, 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 )