f18-7

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 810 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).also comment in passing that the connection between maximum entropy inversion methods, considered here, and maximum entropy spectral estimation, discussed in§13.7, is rather abstract. For practical purposes the two techniques, though both namedmaximum entropy method orMEM, are unrelated. Bayes’Theorem,whichfollowsfromthestandardaxiomsofprobability,relates the conditional probabilities of two events, say AandB: Prob(A|B)=Prob(A)Prob(B|A) Prob(B)(18.7.1 ) HereProb (A|B)is the probabilityof AgiventhatBhas occurred,andsimilarly for Prob(B|A), while Prob (A)and Prob (B)are unconditional probabilities. “Bayesians” (so-called) adopt a broader interpretation of probabilities than do so-called “frequentists.” To a Bayesian, P(A|B)is a measure of the degree of plausibilityof A(givenB)onascalerangingfromzerotoone. Inthisbroaderview, AandBneed not be repeatable events; they can be propositions or hypotheses. The equations of probability theory then become a set of consistent rules for conducting inference [1,2]. Since plausibility is itself always conditioned on some, perhaps unarticulated, set of assumptions, all Bayesian probabilities are viewed as conditional on some collective background information I. Suppose His some hypothesis. Even before there exist any explicit data, a Bayesian can assign to Hsome degree of plausibility Prob (H|I), called the “Bayesian prior.” Now, when some data D1comes along, Bayes theorem tells how to reassess the plausibility of H, Prob(H|D1I)=Prob(H|I)Prob (D1|HI) Prob(D1|I)(18.7.2 ) The factor in the numerator on the right of equation (18.7.2) is calculable as the probability of a data set giventhe hypothesis (compare with “likelihood” in §15.1). The denominator,called the “priorpredictiveprobability”of the data, is in this case merely a normalization constant which can be calculated by the requirement thatthe probability of all hypotheses should sum to unity. (In other Bayesian contexts, the prior predictive probabilities of two qualitatively different models can be used to assess their relative plausibility.) If some additional data D 2comes along tomorrow, we can further refine our estimate of H’s probability, as Prob(H|D2D1I)=Prob(H|D1I)Prob(D2|HD 1I) Prob(D2|D1I)(18.7.3 ) Using the product rule for probabilities, Prob (AB|C)=Prob(A|C)Prob (B|AC), we find that equations (18.7.2) and (18.7.3) imply Prob(H|D2D1I)=Prob(H|I)Prob (D2D1|HI) Prob(D2D1|I)(18.7.4 ) which shows that we would have gotten the same answer if all the data D1D2 had been taken together. 18.7Maximum EntropyImageRestoration 811Sample 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).FromaBayesianperspective,inverseproblemsareinferenceproblems [3,4]. The underlying parameter set uis a hypothesis whose probability, given the measured data values c, and the Bayesian prior Prob (u|I)can be calculated. We might want to report a single “best” inverse u, the one that maximizes Prob(u|cI)=Prob(c|uI)Prob(u|I) Prob(c|I)(18.7.5 ) over all possible choices of u. Bayesian analysis also admits the possibility of reporting additional information that characterizes the region of possible u’s with high relative probability, the so-called “posterior bubble” in u. Thecalculationoftheprobabilityofthedata c,giventhehypothesis uproceeds exactlyasinthemaximumlikelihoodmethod. ForGaussianerrors,e.g.,itisgivenby Prob(c|uI)=e x p ( −1 2χ2)∆u1∆u2···∆uM (18.7.6 ) whereχ2is calculated from uandcusing equation (18.4.9), and the ∆uµ’s are constant,small ranges ofthe componentsof uwhoseactual magnitudeis irrelevant, because they do not depend on u(compare equations 15.1.3 and 15.1.4). In maximum likelihood estimation we, in effect, chose the prior Prob (u|I)to beconstant. Thatwasaluxurythatwecouldaffordwhenestimatingasmallnumberof parameters from a large amount of data. Here, the number of “parameters” (components of u) is comparable to or larger than the number of measured values (components of c); weneedto have a nontrivial prior, Prob (u|I), to resolve the degeneracy of the solution. In maximum entropy image restoration, that is where entropycomes in. The entropy of a physical system in some macroscopic state, usually denoted S, is the logarithm of the number of microscopically distinct configurations that all have the same macroscopic observables (i.e., consistent with the observed macroscopicstate). Actually, we will find it useful to denote the negativeof the entropy, also called the negentropy ,b yH≡−S(a notation that goes back to Boltzmann). In situations where there is reason to believe that the a prioriprobabilities of the microscopic configurationsareallthesame(thesesituationsarecalled ergodic),then the Bayesian priorProb (u|I)foramacroscopic state with entropy Sis proportional toexp(S)orexp(−H). MEM uses this concept to assign a prior probability to any given underlying functionu. For example [5-7], suppose that the measurement of luminance in each pixel is quantized to (in some units) an integer value. Let U=M/summationdisplay µ=1uµ (18.7.7 ) be the total numberof luminancequantain the whole image. Then we can base our “prior” on the notion that each luminance quantum has an equal a priorichance of beinginanypixel. (See [8]foramoreabstractjustificationofthisidea.) Thenumber of ways of getting a particular configuration uis U! u1!u2!···uM!∝exp/bracketleftBigg −/summationdisplay µuµln(uµ/U)+1 2/parenleftBigg lnU−/summationdisplay µlnuµ/parenrightBigg/bracketrightBigg (18.7.8 ) 812 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).Here the left side can be understood as the number of distinct orderings of all the luminance quanta, divided by the numbers of equivalent reorderings withineach pixel, while the right side follows by Stirling’s approximation to the factorial function. Taking the negativeof the logarithm, and neglecting terms of order logU in the presence of terms of order U, we get the negentropy H(u)= M/summationdisplay µ=1uµln(uµ/U)( 18.7.9 ) From equations (18.7.5),(18.7.6),and (18.7.9)we now seek to maximize Prob(u|c)∝exp/bracketleftbigg −1 2χ2/bracketrightbigg exp[−H(u)] ( 18.7.10 ) or, equivalently, minimize: −ln [Prob(u|c)]=1 2χ2[u]+H(u)=1 2χ2[u]+M/summationdisplay µ=1uµln(uµ/U) (18.7.11 ) Thisoughttoremindyouofequation(18.4.11),orequation(18.5.6),orinfactanyof ourpreviousminimizationprinciplesalongthe lines of A+λB, whereλB=H(u) is a regularizing operator. Where is λ? We need to put it in for exactly the reason discussed following equation (18.4.11): Degenerate inversions are likely to be able to achieve unrealistically small values of χ2. We need an adjustable parameter to bringχ2intoitsexpectednarrowstatisticalrangeof N±(2N)1/2. Thediscussionat the beginningof §18.4showed that it makes nodifferencewhich term we attach the λto. Forconsistencyinnotation,weabsorbafactor2into λandputitontheentropy term. (Anotherwaytoseethenecessityofanundetermined λfactoristonotethatit is necessaryif ourminimizationprincipleis tobeinvariantunderchangingtheunitsin which uis quantized, e.g., if an 8-bit analog-to-digitalconverteris replaced by a 12-bit one.) We can now also put “hats” back to indicate that this is the procedure for obtaining our chosen statistical estimator: minimize: A+λB=χ 2[/hatwideu]+λH(/hatwideu)=χ2[/hatwideu]+λM/summationdisplay µ=1/hatwideuµln(/hatwideuµ)(18.7.12 ) (Formally, we might also add a second Lagrange multiplier λ/primeU, to constrain the total intensity Uto be constant.) Itisnothardtoseethatthenegentropy, H(/hatwideu),isinfactaregularizingoperator, similar to /hatwideu·/hatwideu(equation 18.4.11) or /hatwideu·H·/hatwideu(equation 18.5.6). The following of its properties are noteworthy: 1. When Uis held constant, H(/hatwideu)is minimizedfor /hatwideuµ=U/M =constant, so it smoothsinthesenseoftryingtoachieveaconstantsolution,similartoequation (18.5.4). Thefactthattheconstantsolutionisaminimumfollowsfromthefact that the second derivative of ulnuis positive. 18.7Maximum EntropyImageRestoration 813Sample 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).2. Unlike equation (18.5.4), however, H(/hatwideu)islocal, in the sense that it does not difference neighboringpixels. It simply sums some function f, here f(u)=ulnu (18.7.13 ) overallpixels;it isinvariant,infact,underacompletescramblingofthepixels in an image. This form implies that H(/hatwideu)is not seriously increased by the occurrence of a small number of very bright pixels (point sources) embedded in a low-intensity smooth background. 3.H(/hatwideu)goes to infinite slope as any one pixel goes to zero. This causes it to enforcepositivityoftheimage,withoutthenecessityofadditionaldeterministic constraints. 4. The biggest difference between H(/hatwideu)and the other regularizingoperators that we have met is that H(/hatwideu)is not a quadratic functional of /hatwideu, so the equations obtainedbyvaryingequation(18.7.12)are nonlinear . Thisfact is itself worthy of some additional discussion. Nonlinear equations are harder to solve than linear equations. For image processing, however, the large number of equations usually dictates an iterativesolutionprocedure,evenforlinearequations,sothepracticaleffectofthenonlinearity is somewhat mitigated. Below, we will summarize some of the methods that are successfully used for MEM inverse problems. Forsomeproblems,notablytheprobleminradio-astronomyofimagerecovery from an incomplete set of Fourier coefficients, the superior performance of MEM inversioncan be, in part, traced to the nonlinearityof H(/hatwideu). One way to see this [5] is to considerthe limit of perfect measurements σi→0. In this case the χ2term in the minimizationprinciple (18.7.12)gets replacedby a set of constraints, each withits ownLagrangemultiplier,requiringagreementbetweenmodelanddata; that is, minimize:/summationdisplay jλj/bracketleftBigg cj−/summationdisplay µRjµ/hatwideuµ/bracketrightBigg +H(/hatwideu)( 18.7.14 ) (cf.equation18.4.7). Setting the formalderivativewith respectto /hatwideuµto zerogives ∂H ∂/hatwideuµ=f/prime(/hatwideuµ)=/summationdisplay jλjRjµ (18.7.15 ) or defining a function Gas the inverse function of f/prime, /hatwideuµ=G /summationdisplay jλjRjµ  (18.7.16 ) Thissolutionis onlyformal,sincethe λj’s mustbefoundbyrequiringthatequation (18.7.16)satisfy all the constraints built into equation(18.7.14). However,equation(18.7.16)doesshowthecrucialfactthatif Gislinear,thenthesolution /hatwideucontainsonly a linear combination of basis functions R jµcorresponding to actual measurements j. This is equivalent to setting unmeasured cj’s to zero. Notice that the principal solution obtained from equation (18.4.11) in fact has a linear G. 814 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).In the problem of incomplete Fourier image reconstruction, the typical Rjµ has the form exp(−2πikj·xµ), wherexµis a two-dimensionalvector in the image space and kµis a two-dimensional wave-vector. If an image contains strong point sources, then the effect of setting unmeasured cj’s to zero is to produce sidelobe ripples throughout the image plane. These ripples can mask any actual extended,low-intensityimagefeatureslying betweenthe pointsources. If, however,the slope ofGis smaller for small values of its argument,largerfor large values, then ripples in low-intensity portions of the image are relatively suppressed, while strong point sources will be relativelysharpened(“superresolution”). This behavioron the slope ofGis equivalent to requiring f /prime/prime/prime(u)<0.F o rf(u)=ulnu, we in fact have f/prime/prime/prime(u)=−1/u2<0. In more picturesque language, the nonlinearity acts to “create” nonzerovalues for the unmeasured ci’s, so as to suppress the low-intensity ripple and sharpen the point sources. Is MEM Really Magical? How unique is the negentropyfunctional (18.7.9)? Recall that that equation is based on the assumption that luminance elements are a prioridistributed over the pixelsuniformly. Ifweinsteadhadsomeotherpreferred aprioriimageinmind,one with pixel intensities mµ, then it is easy to show that the negentropybecomes H(u)=M/summationdisplay µ=1uµln(uµ/mµ)+constant (18.7.17 ) (theconstantcanthenbeignored). All therest ofthe discussionthengoesthrough. Morefundamentally,anddespitestatementsbyzealotstothecontrary [7], there is actually nothing universal about the functional form f(u)=ulnu. In some otherphysicalsituations (forexample,the entropyofanelectromagneticfieldin the limit of many photons per mode, as in radio-astronomy) the physical negentropy functional is actually f(u)=−lnu(see[5]for other examples). In general, the question,“Entropyofwhat?” is notuniquelyanswerablein anyparticularsituation. (Seereference [9]foranattemptat articulatingamoregeneralprinciplethatreduces to one or another entropy functional under appropriate circumstances.) The four numbered properties summarized above, plus the desirable sign for nonlinearity, f/prime/prime/prime(u)<0, are all as true for f(u)=−lnuas forf(u)=ulnu.I n fact these properties are shared by a nonlinear function as simple as f(u)=−√u, which has no information theoretic justification at all (no logarithms!). MEM reconstructions of test images using any of these entropy forms are virtually indistinguishable [5]. By all available evidence, MEM seems to be neither more nor less than one usefullynonlinearversionofthegeneralregularizationscheme A+λBthatwehave by now considered in many forms. Its peculiarities become strengths when applied to the reconstruction from incomplete Fourier data of images that are expectedto be dominated by very bright point sources, but which also contain interesting low-intensity, extended sources. For images of some other character, there is no reason to suppose that MEM methods will generally dominate other regularization schemes, either ones already known or yet to be invented. 18.7Maximum EntropyImageRestoration 815Sample 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).AlgorithmsforMEM The goal is to find the vector /hatwideuthat minimizes A+λBwhere in the notation of equations (18.5.5), (18.5.6), and (18.7.13), A=|b−A·/hatwideu|2B=/summationdisplay µf(/hatwideuµ)( 18.7.18 ) Compared with a “general” minimization problem, we have the advantage that we can compute the gradients and the second partial derivative matrices (Hessian matrices) explicitly, ∇A =2 (AT·A·/hatwideu−AT·b)∂2A ∂/hatwideuµ∂/hatwideuρ=[ 2AT·A]µρ [∇B]µ=f/prime(/hatwideuµ)∂2B ∂/hatwideuµ∂/hatwideuρ=δµρf/prime/prime(/hatwideuµ)(18.7.19 ) Itisimportanttonotethatwhile A’ssecondpartialderivativematrixcannotbestored (its size is the square of the number of pixels), it can be applied to any vector by first applying A, thenAT. In the case of reconstruction from incomplete Fourier data, or in the case of convolutionwith a translationinvariantpoint spreadfunction, these applications will typically involve several FFTs. Likewise, the calculation of the gradient ∇Awill involve FFTs in the application of AandAT. While some success has been achieved with the classical conjugate gradient method ( §10.6), it is often found that the nonlinearity in f(u)=ulnucauses problems. Attempted steps that give /hatwideuwith even one negative value must be cut in magnitude,sometimessoseverelyastoslowthesolutiontoacrawl. Theunderlying problem is that the conjugate gradient method develops its information about the inverseoftheHessianmatrixabitatatime,whilechangingitslocationinthesearch space. When a nonlinearfunction is quite different from a pure quadratic form, the old information becomes obsolete before it gets usefully exploited. Skilling and collaborators [6,7,10,11] developed a complicated but highly suc- cessful scheme, wherein a minimum is repeatedly sought not along a single search direction,butinasmall-(typicallythree-)dimensionalsubspace,spannedbyvectorsthat are calculated anew at each landing point. The subspace basis vectors are chosen in such a way as to avoid directions leading to negative values. One of the most successful choices is the three-dimensional subspace spanned by the vectors with components given by e (1) µ=/hatwideuµ[∇A]µ e(2) µ=/hatwideuµ[∇B]µ e(3) µ=/hatwideuµ/summationtext ρ(∂2A/∂/hatwideuµ∂/hatwideuρ)/hatwideuρ[∇B]ρ/radicalBig/summationtext ρ/hatwideuρ([∇B]ρ)2−/hatwideuµ/summationtext ρ(∂2A/∂/hatwideuµ∂/hatwideuρ)/hatwideuρ[∇A]ρ/radicalBig/summationtext ρ/hatwideuρ([∇A]ρ)2 (18.7.20 ) (Intheseequationsthereisnosumover µ.) Theformofthe e(3)hassomejustification if one views dot products as occurring in a space with the metric gµν=δµν/uµ, chosen to make zero values “far away”; see [6]. 816 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).Withinthethree-dimensionalsubspace,thethree-componentgradientandnine- component Hessian matrix are computed by projection from the large space, andthe minimum in the subspace is estimated by (trivially) solving three simultaneous linear equations, as in §10.7, equation (10.7.4). The size of a step ∆/hatwideuis required to be limited by the inequality /summationdisplay µ(∆/hatwideuµ)2//hatwideuµ<(0.1to0.5)U (18.7.21 ) Because the gradient directions ∇Aand∇Bare separately available, it is possible tocombinetheminimumsearchwithasimultaneousadjustmentof λsoasfinallyto satisfy the desired constraint. There are various further tricks employed. A less general, but in practice often equally satisfactory, approach is due to Cornwell and Evans [12]. Here, noting that B’s Hessian (second partial derivative) matrix is diagonal, one asks whether there is a useful diagonal approximation to A’s Hessian, namely 2AT·A.I fΛµdenotes the diagonal components of such an approximation, then a useful step in /hatwideuwould be ∆/hatwideuµ=−1 Λµ+λf/prime/prime(/hatwideuµ)(∇A+λ∇B)( 18.7.22 ) (again compare equation 10.7.4). Even more extreme, one might seek an approx- imation with constant diagonal elements, Λµ=Λ, so that ∆/hatwideuµ=−1 Λ+λf/prime/prime(/hatwideuµ)(∇A+λ∇B)( 18.7.23 ) SinceAT·Ahas something of the nature of a doubly convolved point spread function, and since in real cases one often has a point spread function with a sharp central peak, even the more extreme of these approximations is often fruitful. Onestarts with a rough estimate of Λobtained from the A iµ’s, e.g., Λ∼/angbracketleftBigg/summationdisplay i[Aiµ]2/angbracketrightBigg (18.7.24 ) An accurate value is not important, since in practice Λis adjusted adaptively: If Λ is too large, then equation(18.7.23)’ssteps will be too small (that is, larger steps in thesamedirectionwillproduceevengreaterdecreasein A+λB). IfΛis toosmall, thenattemptedstepswilllandinanunfeasibleregion(negativevaluesof /hatwideuµ),orwill resultinanincreased A+λB. Thereisanobvioussimilaritybetweentheadjustment ofΛhere and the Levenberg-Marquardt method of §15.5; this should not be too surprising, since MEM is closely akin to the problem of nonlinear least-squares fitting. Reference [12]also discusses how the value of Λ+λf/prime/prime(/hatwideuµ)can be used to adjust the Lagrangemultiplier λso as to convergeto the desired value of χ2. All practical MEM algorithms are found to require on the order of 30 to 50 iterations to converge. This convergence behavior is not now understood in any fundamental way. “Bayesian” versus “Historic”Maximum Entropy Several more recent developments in maximum entropy image restoration go under the rubric “Bayesian” to distinguish them from the previous “historic” methods. See [13]for details and references. 18.7Maximum EntropyImageRestoration 817Sample 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).•Better priors: We already noted that the entropy functional (equation 18.7.13) is invariant under scrambling all pixels and has no notion ofsmoothness. The so-called “intrinsic correlation function” (ICF) model (Ref. [13], where it is called “New MaxEnt”) is similar enough to the entropy functional to allow similar algorithms, but it makes the values ofneighboring pixels correlated, enforcing smoothness. •Better estimation of λ: Above we chose λto bringχ 2into its expected narrowstatistical rangeof N±(2N)1/2. This in effectoverestimates χ2, however,sincesomeeffectivenumber γofparametersarebeing“fitted”in doing the reconstruction. A Bayesian approach leads to a self-consistentestimate of this γand an objectively better choice for λ. CITED REFERENCES AND FURTHER READING: Jaynes, E.T. 1976, in Foundations of Probability Theory, Statistical Inference, and Statistical Theories of Science , W.L. Harper and C.A. Hooker, eds. (Dordrecht: Reidel). [1] Jaynes,E.T.1985,in Maximum-Entropy andBayesianMethodsinInverseProblems ,C.R.Smith and W.T. Grandy, Jr., eds. (Dordrecht: Reidel). [2] Jaynes, E.T. 1984, in SIAM-AMS Proceedings , vol. 14, D.W. McLaughlin, ed. (Providence, RI: American Mathematical Society). [3] Titterington, D.M. 1985, Astronomy and Astrophysics , vol. 144, 381–387. [4] Narayan, R., and Nityananda, R. 1986, Annual Review of Astronomy and Astrophysics , vol. 24, pp. 127–170. [5] Skilling, J., and Bryan, R.K. 1984, Monthly Notices of the Royal Astronomical Society , vol. 211, pp. 111–124. [6] Burch, S.F., Gull, S.F., and Skilling, J. 1983, Computer Vision, Graphics and Image Processing , vol. 23, pp. 113–128. [7] Skilling,J.1989,in MaximumEntropyandBayesianMethods ,J.Skilling,ed.(Boston:Kluwer).[8] Frieden, B.R. 1983, Journal of the Optical Society of America , vol. 73, pp. 927–938. [9] Skilling,J.,andGull,S.F.1985,in Maximum-EntropyandBayesianMethodsinInverseProblems , C.R. Smith and W.T. Grandy, Jr., eds. (Dordrecht: Reidel). [10] Skilling, J. 1986, in Maximum Entropy andBayesian Methods in AppliedStatistics , J.H. Justice, ed. (Cambridge: Cambridge University Press). [11] Cornwell, T.J., and Evans, K.F. 1985, Astronomy and Astrophysics , vol. 143, pp. 77–83. [12] Gull,S.F.1989,in MaximumEntropyandBayesianMethods ,J.Skilling,ed.(Boston:Kluwer).[13]