Vol.20,No.6,pp.1964–1977
c1999SocietyforIndustrialandAppliedMathematics
ANONLINEARPRIMAL-DUALMETHODFORTOTAL
VARIATION-BASEDIMAGERESTORATION∗
TONYF.CHAN†,GENEH.GOLUB‡,ANDPEPMULET§
Abstract.Wepresentanewmethodforsolvingtotalvariation(TV)minimizationproblemsinimagerestoration.Themainideaistoremovesomeofthesingularitycausedbythenondif-ferentiabilityofthequantity|∇u|inthedefinitionoftheTV-normbeforeweapplyalinearizationtechniquesuchasNewton’smethod.Thisisaccomplishedbyintroducinganadditionalvariableforthefluxquantityappearinginthegradientoftheobjectivefunction,whichcanbeinterpretedasthenormalvectortothelevelsetsoftheimageu.Ourmethodcanbeviewedasaprimal-dualmethodasproposedbyConnandOverton[APrimal-DualInteriorPointMethodforMinimizingaSumofEuclideanNorms,preprint,1994]andAndersen[Ph.D.thesis,OdenseUniversity,Denmark,1995]fortheminimizationofasumofEuclideannorms.Inadditiontopossessinglocalquadraticconvergence,experimentalresultsshowthatthenewmethodseemstobegloballyconvergent.
Keywords.imagerestoration,totalvariation,primal-dual,conjugategradientAMSsubjectclassifications.68U10,65F10,65K10PII.S1064827596299767
1.Introduction.Duringsomephasesofthemanipulationofanimagesomerandomnoiseandblurringisusuallyintroduced.Thepresenceofthisnoiseandblurringmakesthelatterphasesoftheimageprocessingdifficultandinaccurate.
Thealgorithmsfornoiseremovalanddeblurringhavebeenmainlybasedonleastsquares.TheoutputoftheseL2-basedalgorithmswillbeacontinuousfunction,whichcannotobviouslybeagoodapproximationtoouroriginalimageifitcontainsedges.Toovercomethisdifficultyatechniquebasedontheminimizationofthetotalvariationnorm(TV-norm)subjecttosomenoiseconstraintsisproposedin[19],whereatimemarchingschemetosolvetheassociatedEuler–Lagrangeequationswasalsoutilized.Sincethismethodcanbeslowduetostabilityconstraintsinthetimestepsize,anumberofalternativemethodshavebeenproposedintheliterature[23],[7],[16].
OneofthedifficultiesinsolvingtheEuler–Lagrangeequationsisthepresenceofahighlynonlinearandnondifferentiableterm,whichcausesconvergencedifficultiesfor
bytheeditorsFebruary26,1996;acceptedforpublication(inrevisedform)November
7,1997;publishedelectronicallyMay25,1999.AnearlierversionofthisworkappearedasTechnicalreport95-43,UCLAComputationalandAppliedMathematics,UniversityofCalifornia,LosAngeles,CA,1995.
http://www.siam.org/journals/sisc/20-6/29976.html
†DepartmentofMathematics,UniversityofCalifornia,LosAngeles,405HilgardAvenue,LosAngeles,CA90095-1555(chan@math.ucla.edu).TheworkofthisauthorwassupportedbygrantsNSFDMS-9626755andONR-N00014-96-1-0277andNSFInternationalProgramU.S.-SpainResearchProjectTotalVariationMethodsinImageProcessing.
‡ComputerScienceDepartment,StanfordUniversity,Stanford,CA94305-9025(golub@sccm.stanford.edu).TheworkofthisauthorwassupportedinpartbyNSFgrantCCR-9505393andNSFInternationalProgramU.S.-SpainResearchProjectTotalVariationMethodsinImageProcessing.
§DepartmentofMathematics,UniversityofCalifornia,LosAngeles,405HilgardAvenue,LosAngeles,CA90095-1555(mulet@math.ucla.edu).Currentaddress:DepartamentdeMatem`aticaAplicada,UniversitatdeVal`encia,Dr.Moliner,50,46100Burjassot,Spain(mulet@uv.es,mulet@gata.matapl.uv.es).TheworkofthisauthorwassupportedbyDGICYTgrantsEX9428990695andPB94-0987andNSFInternationalProgramU.S.-SpainResearchProjectTotalVaria-tionMethodsinImageProcessing.
1964
∗Received
APRIMAL-DUALMETHODFORTVIMAGERESTORATION
1965
Newton’smethodevenwhencombinedwithaglobalizationtechniquesuchasalinesearch.Theideaofournewalgorithmistoremovesomeofthesingularitycausedbythenondifferentiabilityoftheobjectivefunctionbeforeweapplyalinearizationtech-niquesuchasNewton’smethod.Thisisaccomplishedbyintroducinganadditionalvariableforthefluxquantityappearinginthegradientoftheobjectivefunction,whichcanbeinterpretedastheunitnormaltothelevelsetsoftheimagefunction.Ourmethodcanbeviewedasaprimal-dualmethodasproposedbyConnandOverton[9]andAndersen[3]fortheminimizationofasumofEuclideannorms.Wecautionthatouruseofthenameprimal-dualisbasedonadualityprincipleappliedtotheTV-norm(tobeexplainedinmoredetailinsection4)andshouldnotbeconfusedwiththepopularalgorithmsinlinearandnonlinearprogrammingwiththesamename.Experimentalresultsshowthatthenewmethodisgloballyconvergent,whereastheprimalNewtonmethodhasasmalldomainofconvergence.Itishopedthatthenewapproachcanbeappliedtoothergeometry-basedPDEmethodsinimagerestoration,suchasanisotropicdiffusion[18],affineinvariantflows[20],andmeancurvatureflows[2],sincethesamesingularitycausedby|∇u|occursintheseformulationsaswell.
Theorganizationofthispaperisasfollows:Insection2weintroducetheproblemandthenonlinearequationsassociatedwithitanddiscusshowtosolvethem.Insection3wepresentournewlinearizationtechniqueforthe(unconstrained)Tikhonovregularizationformoftheproblem.Insection4weshowthatthenewlinearizationcorrespondstoaclassicalprimal-dualformulationforsomeconvexproblems.Finally,insection5wepresentsomenumericalresultsforthedenoisingcase.
2.Totalvariationregularization.Animagecanbeinterpretedasarealfunc-tiondefinedonΩ,aboundedandopendomainofR2(forsimplicitywewillassumeΩtobearectangle),orasuitablediscretizationofthiscontinuousimage.
Ourinterestisinrestoringanimagewhichiscontaminatedwithnoiseorblurorboth.Therestorationprocessshouldrecovertheedgesoftheimage.Letusdenotebyztheobservedimageandbyutherealimage.ThemodelofdegradationweassumeisKu+n=z,wherenisaGaussianwhitenoiseandKisa(known)linearbluroperator(usuallyaconvolutionoperator).
Ingeneral,theproblemKu=z,withKacompactoperator,isillposed,soitisnotworthsolvingthisequation(oradiscretizationofit),forthedataintheusualapplicationsisinexactandthesolutionwouldbehighlyoscillatory.Butifweimposeacertainregularityconditiononthesolutionu,thentheproblembecomeswellposed.Wecanconsidertworelatedtechniquesofregularization:Tikhonovregularizationandnoiselevelconstrainedregularization.
Tikhonovregularizationconsistsofsolvingtheunconstrainedoptimizationprob-lem(2.1)
1
minαR(u)+||Ku−z||2L2u2
forsomefunctionalRwhichmeasurestheirregularityofuinacertainsenseanda
suitablychosencoefficientαwhichwillbalancethetradeoffbetweenagoodfittothedataandaregularsolution.
Anotherapproachconsistsofsolvingthefollowingconstrainedoptimizationprob-lem:(2.2)
minR(u)
u
2
subjectto||Ku−z||2L2=σ.
1966
TONYF.CHAN,GENEH.GOLUB,ANDPEPMULET
Hereweseekasolutionwithminimumirregularityfromallcandidateswhichmatchtheknownnoiselevel.
ExamplesofregularizationfunctionalsthatcanbefoundintheliteratureareR(u)=||u||,||∆u||,||∇u·∇u||,where∇isthegradientand∆istheLaplacian;see[21]and[14].Thedrawbackofusingthesefunctionalsisthattheydonotallowdiscontinuitiesinthesolution,andsinceweareinterestedinrecoveringfeaturesoftheimage,theyarenotsuitableforourpurposes.
In[19],itisproposedtouseasregularizationfunctionaltheso-calledTV-norm:
2(2.3)u2|∇u|dxdy=TV(u)=x+uydxdy.
Ω
Ω
TheTV-normdoesnotpenalizediscontinuitiesinuandthusallowsustorecover
theedgesoftheoriginalimage.ForsimplicityweuseinthissectiontheTikhonovformulationoftheproblem.Treatmentoftheconstrainedproblem(2.2)canbefoundin[8].Hencetherestorationproblemcanbewrittenas
122αu2(2.4)dxdy,minx+uy+(Ku−z)u2Ωthatis,(2.5)
minα|||∇u|||L1
u
1
+||Ku−z||2L2,2
TheEuler–Lagrangeequationforthisproblem,assuminghomogeneousNeumann
boundaryconditions,is
∇u(2.6)+K∗(Ku−z),0=−α∇·
|∇u|whereK∗istheadjointoperatorofKwithrespecttotheL2innerproduct.Thisequationisnotwelldefinedatpointswhere∇u=0,duetothepresenceoftheterm1/|∇u|.AcommonlyusedtechniquetoovercomethisdifficultyistoslightlyperturbtheTV-normfunctionaltobecome
(2.7)|∇u|2+βdxdy,
Ω
2|∇u|=u2x+uy.
whereβisasmallpositiveparameter.In[1]itisshownthatthesolutionsofthese
perturbedproblemsconvergetothesolutionof(2.4)whenβ→0.Sonowtheproblemis
1
|∇u|2+βdxdy+||Ku−z||2minα(2.8)L2,u2Ωwhichcanbeshowntobestrictlyconvexifα>0andKisinvertible,whichwewill
assumehenceforth.ItscorrespondingEuler–Lagrangeequationis
∇u(2.9)+K∗(Ku−z)=g(u).0=−α∇·|∇u|2+β|∇u|+β
Themaindifficultythatthisequationposesisthelinearizationofthehighlynonlinearterm−∇·(√∇u2).
APRIMAL-DUALMETHODFORTVIMAGERESTORATION
1967
Anumberofmethodshavebeenproposedtosolve(2.9).Rudin,Osher,andFatemi[19]usedatimemarchingschemetoreachthesteadystateoftheparabolicequationut=−g(u)withinitialconditionu=z:
∇u(2.10)−K∗(Ku−z),u(x,0)=z(x).ut=α∇·|∇u|2+βdoesnotworksatisfactorilyinthesensethatitsdomainofconvergenceisverysmall.Thisisespeciallytrueiftheregularizingparameterβissmall.Ontheotherhand,ifβisrelativelylarge,thenthistermiswellbehaved.Soitisnaturaltouseacontinuationprocedurestartingwithalargevalueofβandgraduallyreducingittothedesiredvalue.Chan,Chan,andZhouproposedsuchanapproachin[7].Althoughthismethodislocallyquadraticallyconvergent,thechoiceofthesequenceofsubproblemstosolveiscrucialforitsefficiency.Theauthorshavenotsucceededinfindingafullysatisfactoryselectionprocedure,althoughsomeheuristicscanbeused.
3.Anewlinearizationbasedonadualvariable.Weproposehereabetter
∇utechniquetolinearizetheterm∇·(|∇u|).Thistechniquebearssomesimilaritytotechniquesfromsomespecificnonlinearprimal-dualoptimizationmethods[3],[9]andgivesabetterglobalconvergencebehaviorthanthatoftheusualNewton’scontinua-tionmethod.Inthenextsectionwegiveajustificationoftheprimal-dualnatureofthealgorithmthatwepresentbelow.
Themethodisbasedonthefollowingsimpleobservation.Whilethesingularityandnondifferentiabilityofthetermw=∇u/|∇u|isthesourceofthenumericalproblems,witselfisusuallysmoothbecauseitisinfacttheunitnormalvectortothelevelsetsofu.Thenumericaldifficultiesariseonlybecausewelinearizeitinthewrongway.
Theideaofthenewmethodistointroduce(3.1)
w=|∇u|2+β∇uStartingwithu0=z,ateachstepuk+1isobtainedasthesolutionofthelineardifferential-convolutionequation(2.11),whosecoefficientsarecomputedfromuk.Thismethodisrobustbutonlylinearlyconvergent.
Duetothepresenceofthehighlynonlinearterm∇·(√∇u2),Newton’smethod
|∇u|+β
Thismethodcanbeslowlyconvergentduetostabilityconstraints.VogelandOman[23]proposedthefollowingfixedpointiterationtosolvetheEuler–Lagrangeequation:
k+1∇u(2.11)−α∇·+K∗(Kuk+1−z)=0.
k2|∇u|+βasanewvariableandreplace(2.9)bythefollowingequivalentsystemofnonlinear
PDEs:(3.2)
−α∇·w+K∗(Ku−z)=0,
w|∇u|2+β−∇u=0.
Wecanthenlinearizethis(u,w)system,forexample,byNewton’smethod.Thisapproachissimilartothetechniqueofintroducingafluxvariableinthemixedfinite
1968
TONYF.CHAN,GENEH.GOLUB,ANDPEPMULET
Table3.1
ComparisonofthenumberofiterationsrequiredbyNewton’smethodtosolvef(x)=0andg(x)=0,fora=0.9999,fordifferentβ(horizontally)anddifferentinitialguessesx0(vertically).A∗meansthatthecorrespondingiterationfailedtoconverge.
Newton’siterationforg(x)=010−210−310−410−57643565246643565435753467Newton’siterationforf(x)=010−210−310−410−5*966**84**1075**96***85***x0\\β
1234510−110−1elementmethod[4].Aswewillseeinthenextsection,theseequationsare,precisely,theequationsofasaddlepointproblemforwhichtheproblem(2.4)istheprimalproblemandwisthevariableofthedualproblem.
Forcompletenesswecomparethelinearizationoftheusystem,
1∇u∇uT
−α∇·I−(3.3)∇+K∗Kδu=−g(u),2|∇u||∇u|tothelinearizationofthe(w,u)system,w∇uT
δwf(w,u)|∇u|−I−|∇u|∇
=−(3.4).δug(w,u)−α∇·K∗K
Equation(3.4)canbesolvedbyfirsteliminatingδwandsolvingtheresultingequation
forδu:
w∇uT1
(3.5)I−∇+K∗Kδu=−g(u).−α∇·
|∇u||∇u|Afterδuisobtainedwecancomputeδwby
1w∇uT∇u(3.6)δw=I−∇δu−w+.
|∇u||∇u||∇u|
Wenotethatthecostperiterationofournewlinearizationtechniqueisonlyslightly
higherthanforthestandardNewton’smethod(3.3),becausethemaincostisthesolutionofthedifferential-convolutionequations(3.3)and(3.5)forδu.
Themotivationisthatthe(w,u)systemissomehowbetterbehavedthantheusystem.Althoughatthispointwedonothaveacompletetheorytosupportthis,wewillnowgiveascalarexamplethatcanexplainthebetterconvergencebehaviorofthenewapproach.equationsWecompareNewton’smethodappliedtotheequivalent
22+β)=0(whichresembles(3.1))andg(x)=ax+β−x=0f(x)=a−(x/x
(whichresemblesw|∇u|2+β−∇u=0),wherea≈1andβ≈0.InFigure3.1wecanseethatglooksmore“linear”thanfovermuchofthex-axis.Inparticular,Newton’smethodappliedtofdivergeswhentheinitialguessisnotclosetothesolution,whereasitconvergeswhenappliedtogandforanypositiveinitialguess.ThisisconfirmedbythenumericalresultsshowninTable3.1.WebelievethatthereasonwhythealgorithmpresentedhereshowssuchadramaticconvergenceimprovementoverthestandardNewton’smethodispreciselythisbetterlinearization.
APRIMAL-DUALMETHODFORTVIMAGERESTORATION
x 10−31969
10.80.60.40.20g(x)=0.9999*sqrt(x^2+1e−4)−xf(x)=0.9999−x/sqrt(x^2+1e−4)f(x)−0.2−0.4−0.6−0.8−10g(x)0.511.522.533.544.55Fig.3.1.Plotoff(x)andg(x).They-axishasbeenscaledbyafactorof10−3.4.Duality.Thegoalofthissectionistoshowthatproblem(2.4)istheprimalproblemforacorrespondingdualproblemandthatthesolutionforbothischarac-terizedbypreciselythesystemofequations(3.2)involvingboththeprimalandthedualvariables.
Ascanbeseenin[12],forinstance,theTV-normadmitsthefollowingdualformulation:
∞
(4.1)−u∇·wdxdy:w∈(C0(Ω))2,w∞≤1,TV(u)=sup
Ω
∞
(Ω)isthesetofsmoothfunctionsonΩthatvanishontheboundary∂Ω.whereC0
Itisfairlyeasytoseethatthisextends(2.3)tonondifferentiableu.Actually,forubelongingtotheSobolevspaceW1,1(Ω),anargumentusingintegrationbypartsshowsthatthesolutionwtotheproblem
u∇·wdxdysup
|w|≤1
Ω
satisfiesinthelimit
w(x)=
∇u(x)|∇u(x)|
foralmostallx∈Ωsuchthat∇u(x)=0or,equivalently,(4.2)
|∇u|w−∇u=0
almosteverywhere.Thisispreciselythesecondequationin(3.2).
Withthisformulation,problem(2.4)canbewrittenas(4.3)
minsupΦ(u,w),
u
w∞≤1
1970where
TONYF.CHAN,GENEH.GOLUB,ANDPEPMULET
Φ(u,w)=
Ω
1
−αu∇·w+(Ku−z)2
2
dxdy.
Byusingargumentsofdualitytheoryforconvexprogramming(see[11],forin-stance),wehavethatasolution(u∗,w∗)of(4.3)isasaddlepointforΦ(u,w),thatis,(4.4)
minsupΦ(u,w)=Φ(u∗,w∗)=
u
w∞≤1
supminΦ(u,w),
w∞≤1
u
and,moreover,w∗isasolutionofthedualproblem(4.5)where(4.6)
Ψ(w)=minΦ(u,w).
u
supΨ(w),
w∞≤1
SinceΦisquadraticinu,wecaneasilyderivethatthisminimummustsatisfythefollowingequation:(4.7)
−α∇·w+K∗(Ku−z)=0,
whereK∗denotestheadjointoperatorofK.Thisispreciselythefirstequationin(3.2).
5.Numericalresults.DenotebyCandCstandardfinitedifferencediscretiza-tionsofthelinearoperators
∇u∇uT1
I−∇+K∗K,(5.1)−α∇·2|∇u||∇u|
w∇uT1
(5.2)I−∇+K∗K−α∇·
|∇u||∇u|thatappearin(3.3)and(3.5),respectively.Duetothesizeofthesematricesweuse
iterativemethodstosolvethediscretizationof(3.3)and(3.5).
ThematrixCissymmetricandpositivedefinitewhenthenaturalassumptionsα>0andKbeinginvertiblehold(wewillassumesohenceforth),thustheconjugategradientcanbeappliedtosolvethediscretizationof(3.3).Itcanbeshownthat
˜=1(C+CT)ispositivethematrixCispositivedefinite(i.e.,itssymmetrizationC2definite)ifα>0,Kisinvertible,and|wi|≤1atallgridpoints(see[8]).However,itisgenerallynonsymmetric,thustheconjugategradientmethodcannotbeapplied
˜,i.e.,useanapproximatetoit.OurapproachistoreplaceCbyitssymmetrizationC
˜canberegardedasadiscretizationofNewton’smethod.ThesymmetrizedmatrixC
theoperator:
1w∇uT+∇uwT1
I−∇+K∗K.−α∇·
|∇u|2|∇u|FromthetheoryoftheconvergenceoftheapproximateNewton’smethod,it
˜convergesfollowsthattheconvergenceoftheresultingiterationisquadratic,sinceC
∇u
locallytoC(sincew→|∇u|),atleastlinearly;see[15,Section5.4.1].
APRIMAL-DUALMETHODFORTVIMAGERESTORATION
1971
u=z,βsucc=βcurr=||∇z||2∞,ρ<1whileβcurr>βtarget
[unew,success]=newton(u,z,βcurr)ifsuccess,thenβsucc=βcurr
βcurr=βcurr∗ρ%decreaseβcurr
ρold=ρ,chooseρ∈(0,ρold)%decreaseρu=unewelse
βcurr=βsucc%increaseβcurr
ρold=ρ,chooseρ∈(ρold,1)%increaseρendend
Fig.5.1.Pseudocodeforcontinuationalgorithm.βcurristheparameterβfortheproblemtobesolvedatthecurrentiteration,βsuccistheβusedinthelastsuccessfulsolutionofaNewton’smethod,ρisthereductionfactorforβcurr.Theinitialchoiceofβcurr=||∇z||2∞isbasedupontheheuristicthatthedomainofconvergenceoftheNewton’smethodforalargeβ(relativetotheterm|∇z|2appearingintheformulationoftheTV-norm)shouldbelarge.
˜canbeensuredifweimposetheAsmentionedabove,thenonsingularityofC
condition|wi|≤1∀ithroughouttheprocess.Thiscanbeachievedinexpensivelybychoosingasteplengths>0suchthat(5.3)
s=ρsup{α:|wi+αδwi|<1∀i},
ρ∈(0,1),
whereitisassumedthat|wi|<1∀i.Actually,ourexperiencerevealsthatthischoiceofsteplengthforwalreadyproducesagloballyconvergentalgorithm,evenforsmallβandwithouttheuseofalinesearchforthevariableu.However,toensureglobalconvergenceif(5.3)werenotused,asteplengthprocedureontheprimal
˜issymmetricpositivedefinite,δu=−(C˜)−1g(u)isvariablesumustbeused.SinceC
adescentdirectionforψ,theobjectivefunctionof(2.8).Webrieflyshownextthatthisdirectionisgradientrelated.Thecosineoftheangleθbetweeng=∇ψandδu
δuTg|
iscosθ=|δug,where·denotestheEuclideannorm.Byusingthesymmetry˜andstandardlinearalgebraproperties,wededuceofC(5.4)(5.5)(5.6)
˜−1)g2˜)˜−1g|λmin(Cλmin(C|gTC|δuTg|
,≥==cosθ=
˜−1gg˜−1)g2˜)δugCλmax(Cλmax(C
˜)=minvTCv,˜λmin(C
v=1
˜)=maxvTCv.˜λmax(C
v=1
˜thatItfollowseasilyfromtheNeumannboundaryconditionsandthedefinitionofC
˜)=λmin(K∗K)>0,bytheassumptionofthenonsingularityofK.Ontheλmin(C
otherhand,itcanbeshownthatthereexistsaconstantM>0,whichdependsonthegridsizeusedinthediscretization,suchthat
M˜≤√+λmax(K∗K)v2,vTCv
β1972
TONYF.CHAN,GENEH.GOLUB,ANDPEPMULET
[u,success]=newton(u,z,β)success=0w=0
fori=1:#iterations
computeupdatedirection(δu,δw)iflinesearchonurequested,thencomputespfrom(5.8)else
sp=1end
iflinesearchonwrequested,thencomputesdfrom(5.3)else
sd=1end
u=u+spδuw=w+sdδw
ifstoppingcriterionissatisfied,thensuccess=1returnendend
Fig.5.2.Pseudocodeforprimal-dualNewtonalgorithmwithpossiblelinesearches.TheprimalNewtonalgorithmisbasicallyidentical,withtheexceptionthatthedualvariableswanddualupdatesδwdonotappear.Thevariableslinesearchonurequestedandlinesearchonwrequestedareglobalvariablessuppliedbytheusertousethecorrespondinglinesearches.
M˜)≤√+λmax(K∗K);therefore,from(5.4)so(5.6)impliesthatλmax(C
β
(5.7)cosθ≥
λmin(K∗K)
>0∀u,M√+λmax(K∗K)
β
thusδuisagradientrelateddescentdirection.Hence,alinesearchbasedona
sufficientdescentofψcouldbeusedtoensureglobalconvergence.Thelinesearchprocedurethatwehaveusedconsistsofabacktrackingalgorithmthatselectsassteplengththefirstnumbersinthesequenceofpowersof2ofdecreasingexponent11
,4,...thatsatisfiestheinequality1,2(5.8)
ψ(u+sδu)≤ψ(u)+ηs∇ψ(u)Tδu,
η∈(0,1).
Inadditiontothelinesearchprocedure,wemayalsouseacontinuationprocedureontheparameterβtoimprovetheglobalconvergencebehaviorofthealgorithm.Thisprocedureconsistsofselectingasufficientlylargevalueofβandgraduallyreducingittothedesiredvalueβtarget.ThepseudocodesforthecontinuationandtheNewton’salgorithmsaregiveninFigures5.1and5.2,respectively.
OurfirstexperimentconsistsofthecomparisonoftheprimalNewtonandtheprimal-dualNewtonmethodsunderthefollowingcircumstances:
1.Continuationonβandnolinesearch.
APRIMAL-DUALMETHODFORTVIMAGERESTORATION
1973
2.Continuationonβandlinesearchonu.
3.Nocontinuationonβanduselinesearchonu.
4.Nocontinuationanduselinesearchonuandwfortheprimal-dualmethod.Fig.5.3.(Left)Originalimage,256×256pixels.(Right)Noisyimage,SNR≈1,||g||=2.07.
Theoriginalimage,whichis256×256pixelsandhasdynamicrange[0,255],appearsinFigure5.3(a).AGaussianwhitenoisewithvarianceσ2≈1200isaddedtoit,resultingintheimagedisplayedinFigure5.3(b)withSNR=(||u−u||L2/σ)≈1.Figure5.4depictsthesolutionobtainedbytheprimal-dualNewtonmethod.
WehavesettheparameterαintheTikhonovformulationtotheinverseoftheLagrangemultiplieryieldedbyapreviousrunoftheconstrainedproblemsolver,inthiscaseα=1.18.Theparameterβhasbeensetto0.01.Furthermore,wehaveusedtruncatedversionsofNewton’salgorithmbasedontheconjugategradientmethodwithincompleteCholeskyaspreconditioner.Thestoppingcriterionforthe(outer)Newton’siterationisarelativedecreaseofthenonlinearresidualbyafactorof10−4.Thestoppingcriterionforthenthinnerlineariterationisarelativedecreaseofthelinearresidualbyafactorofηn,wherewefollowthesuggestionof[15,Eq.6.18]andset
0.1ifn=0,
ηn=(5.9)
min(0.1,0.9||gn||2/||gn−1||2),wheregndenotesthegradientoftheobjectivefunction,i.e.,theright-handsideofthediscretizationof(3.5)or(3.3),atthenthiteration.InTable5.1wecomparetheprimal-dualandtheprimalversionsofNewton’smethodfortheexperimentsdescribedbelow.
Table5.1
Comparisonofprimalandprimal-dualNewtonmethods.
Primal-dualNewtonNWT
Continuation,nolinesearchContinuation,linesearchonuContinuation,linesearchonu&wNocontinuation,linesearchonu&w
25251312
CG1861875158
PrimalNewtonNWT705353
CG265210210
Notconverged
Theconclusionsthatcanbedrawnfromthisexperimentareasfollows:
1974
TONYF.CHAN,GENEH.GOLUB,ANDPEPMULET
Fig.5.4.Denoisedimage,12Newton’siterations,58CGiterations,||g||=5.57×10−5.
•Themostcrucialfactorfortheprimal-dualmethodiscontrollingthedualvariablesviathesteplengthalgorithmappearingin(5.3).Infact,ourexpe-rienceisthatthisalgorithmwiththedualsteplengthisgloballyconvergentfortheparametersαandβinareasonablerange.Alinesearchfortheprimalvariablesalmostalwaysyieldsunitsteplengths.
•Theprimal-dualmethodwiththedualsteplengthalgorithmdoesnotneedcontinuationtoconverge,althoughusingitmightbeslightlybeneficialintermsofwork.
•Theprimal-dualmethodwiththedualsteplengthhasamuchbetterconver-gencebehaviorthantheprimalmethod.
Inoursecondexperiment,wecomparetheprimal-dualNewton,fixedpoint,andtimemarchingmethods.Wehaveimplementedthefixedpoint(2.11)andthetimemarching(2.10)algorithmsusingthesamestandardfinitedifferenceschemeforthediscretizationofthedifferentialoperatorsasfortheprimal-dualmethod.Forthismethod,weusethesteplengthalgorithmforthedualvariables,nocontinuation,andthesameparametersasinthepreviousexperiment.Theseparametersarealsousedforthefixedpointmethod,exceptthatwehaveusedafixedlinearrelativeresidualdecreaseηn=0.1(itis,inthiscase,optimalaccordingtoourexperience).Forthetimemarchingmethodwehaveusedalinesearchbasedonsufficientdecreaseoftheobjectivefunction,inotherwords,anadaptivechoicefortheincrement∆ttomaintaintheiterationstable.Thestoppingcriterionforthetimemarchingmethodisbasedontheiterationcountsincewehavenotbeenabletoachievetheprescribedaccuracyinareasonableamountoftime.InFigures5.5,5.6,and5.7weplottheconvergencehistoryofthisexperiment.
Theconclusionswedrawfromthisexperimentarethefollowing:
•Theprimal-dualalgorithmisquadraticallyconvergent,whereastheothersareatbestlinearlyconvergent.
•Theprimal-dualalgorithmbehavessimilarlytothefixedpointmethodintheearlystages,butinafewiterationscanattainhighaccuracy.
•Thecostperiterationoftheprimal-dualmethodisbetween30and50percentmorethanforthefixedpointiteration.Thememoryrequirementsroughlysatisfythisaswell.
Inthelastexperimentwehavechosenthefollowingparametersfortheprimal-dualmethod:samedegradedimageasinthepreviousexperiments,α=1,β=0.01(nocontinuation);10−8asouteriterationrelativetolerance;10−4asCGrelativetolerance.
APRIMAL-DUALMETHODFORTVIMAGERESTORATION
11975
10100Primal−dual Fixed Point Time marching10−1||gradient||10−210−310−410−501020304050Iterations60708090100Fig.5.5.PlotoftheL2-normofthegradientg(u)oftheobjectivefunctionversusiterationsforthedifferentmethods.
103102Primal−dual Fixed Point Time marchingDifference true image10110010−110−210−301020304050Iterations60708090100Fig.5.6.PlotoftheL2-normofthedifferencebetweenthecurrentiterateandthesolutionfortheproblemcomputedbyNewton’smethodwithhighaccuracyversusiterationsforthedifferentmethods.
Wehavethengenerated50randominitialapproximationszi,i=1,...,50,andobtainedfromtheprimal-dualalgorithmwithinitialguesszianumberofiteratesui0=zi,...,uini.Inalltheseapplicationsofthealgorithmwehaveobtainedconvergenceandthenumberniofiterationshasrangedfrom14to17.Althoughwedonothaveaprooffortheglobalconvergenceofthisalgorithm,theseresultsstrengthenourconjectureoftheglobalconvergenceoftheprimal-dualmethod.
1976
TONYF.CHAN,GENEH.GOLUB,ANDPEPMULET
10010−1% pixels10−2Primal−dual Fixed Point Time marching10−310−401020304050Iterations60708090100Fig.5.7.Plotof#pixelswhichdiffermorethan.001(relatively)fromthesolutionfortheproblemcomputedbyNewton’smethodwithhighaccuracyversusiterationsforthedifferentmethods.
Acknowledgments.WewouldliketothankAndyConnandMichaelOvertonformakingtheirpreprint[9]availabletousandJunZouformanyhelpfulconver-sationswiththefirstauthoronmixedfiniteelements.ThethirdauthorwantstoheartilythankPacoAr´andiga,VicenteCandela,RosaDonat,andAntonioMarquinafortheircontinuoussupport.
REFERENCES
[1]R.AcarandC.R.Vogel,Analysisoftotalvariationpenaltymethodsforill-posedproblems,
InverseProblems,10(1994),pp.1217–1229.
[2]L.Alvarez,P.-L.Lions,andJ.-M.Morel,Imageselectivesmoothingandedgedetectionby
nonlineardiffusionII,SIAMJ.Numer.Anal.,29(1992),pp.845–866.
[3]K.D.Andersen,MinimizingaSumofNorms(LargeScaleSolutionofSymmetricPositive
DefiniteLinearSystems),Ph.D.thesis,OdenseUniversity,Odense,Denmark,1995.
[4]F.Brezzi,Asurveyofmixedfiniteelementmethods,inFiniteElements,ICASE/NASALaRC
Ser.,Springer-Verlag,Berlin,1988,pp.34–49.
[5]P.H.CalamaiandA.R.Conn,Astablealgorithmforsolvingthemultifacilitylocation
probleminvolvingEuclideandistances,SIAMJ.Sci.Statist.Comput.,1(1980),pp.512–526.
[6]P.H.CalamaiandA.R.Conn,Asecond-ordermethodforsolvingthecontinuousmultifacility
locationproblem,inNumericalAnalysis:ProceedingsoftheNinthBiennialConference,Dundee,Scotland,G.A.Watson,ed.,LectureNotesinMath.912,Springer-Verlag,Berlin,1982,pp.1–25.
[7]R.H.Chan,T.F.Chan,andH.M.Zhou,AdvancedSignalProcessingAlgorithms,inPro-ceedingsoftheInternationalSocietyofPhoto-OpticalInstrumentationEngineers,F.T.Luk,ed.,SPIE,1995,pp.314–325.
[8]T.Chan,G.Golub,andP.Mulet,ANonlinearPrimal-DualMethodforTotalVariation-BasedImageRestoration,TechnicalReport95-43,UniversityofCalifornia,LosAngeles,CA,1995.
[9]A.R.ConnandM.L.Overton,APrimal-DualInteriorPointMethodforMinimizinga
SumofEuclideanNorms,preprint,1994.
[10]D.DobsonandF.Santosa,RecoveryofBlockyImagesfromNoisyandBlurredData,Tech-nicalReport94-7,CenterfortheMathematicsofWaves,UniversityofDelaware,Newark,
APRIMAL-DUALMETHODFORTVIMAGERESTORATION
1977
[11][12][13][14][15][16][17][18][19][20][21][22][23]
DE,1994.
I.EkelandandR.Temam,Analyseconvexeetprobl`emesvariationnels,Dunod,Paris,1974.E.Giusti,MinimalSurfacesandFunctionsofBoundedVariations,Birkh¨auser-Verlag,Basel,
Switzerland,1984.
G.GolubandC.vanLoan,MatrixComputations,2nded.,TheJohnsHopkinsUniversity
Press,Baltimore,MD,1989.
C.W.Groetsch,TheTheoryofTikhonovRegularizationforFredholmIntegralEquationsof
theFirstKind,Pitman,Boston,1984.
C.T.Kelley,IterativeMethodsforLinearandNonlinearEquations,FrontiersAppl.Math.
16,SIAM,Philadelphia,PA,1995.
Y.LiandF.Santosa,AnAffineScalingAlgorithmforMinimizingTotalVariationinImage
Enhancement,TechnicalReport12/94,CenterforTheoryandSimulationinScienceandEngineering,CornellUniversity,Ithaca,NY,1994.
M.E.Oman,Fastmultigridtechniquesintotalvariation-basedimagereconstruction,inPre-liminaryProceedingsofthe1995CopperMountainConferenceonMultigridMethods,toappear.
P.PeronaandJ.Malik,Scalespaceandedgedetectionusinganisotropicdiffusion,IEEE
Trans.PatternAnal.Mach.Intelligence,12(1990),pp.629–639.
L.Rudin,S.Osher,andE.Fatemi,Nonlineartotalvariationbasednoiseremovalalgorithms,
Phys.D,60(1992),pp.259–268.
G.SapiroandA.Tannenbaum,Areaandlengthpreservinggeometricinvariantscale-spaces,
inProceedingsThirdEuropeanConferenceonComputerVision,Stockholm,Sweden,Lec-tureNotesinComput.Sci.801,Springer,Berlin,1994,pp.449–458.
A.N.TikhonovandV.Y.Arsenin,SolutionsofIll-PosedProblems,JohnWiley,NewYork,
1977.
C.R.Vogel,Amultigridmethodfortotalvariation-basedimagedenoising,inComputation
andControlIV,Progr.SystemsControlTheory20,Birkh¨auser,Boston,MA,1995,pp.323–331.
C.R.VogelandM.E.Oman,Iterativemethodsfortotalvariationdenoising,SIAMJ.Sci.
Comput.,17(1996),pp.227–238.
因篇幅问题不能全部显示，请点此查看更多更全内容