Kazufumi ITOBangti JINTomoya TAKEUCHI

1 Introduction

We investigate a regularization technique for solving linear inverse problems modeled by

whereg†is the(inaccessible)exact data andu†∈Xrepresents the unknown exact solution,andK:X→Yis a bounded linear operator.Here the spacesXandYare general Banach spaces,and the operatorKcan be an embedding operator(image denoising),a convolution operator(deblurring,scattering)and the Radon transform(computed tomography).The objective is to find an approximationuto the solutionu†from noisy measurementgδ∈Yof the exact datag†.The accuracy of the noisy datagδis measured by the standard L2fidelity functionalwith the noise levelδ.

As is typical for many inverse problems,problem(1.1)suffers from ill-posedness or instability.This poses significant challenges to their accurate yet stable numerical solution in the presence of data noise,which is often the case in practical applications.Often,regularization is applied to find a stable approximate solution.One of the most widely used approaches is known as Tikhonov regularization.It seeks to minimize the following functional

over a closed convex feasible solution setC.The solution to the minimization problem,denoted by(uηin case of the exact datag†),serves as an approximation to the exact solutionu†.Here the(nonnegative)vector-valued penalty functionalψencodes the a priori knowledge,andη·ψ(u)denotes the dot product between the regularization parameter vectorand the penaltyψ(u)=(ψ1(u),ψ2(u))t.The penaltyψis selected to promote desirable features of the sought-for solution,e.g.,edge,sparsity and texture;and often the optimization problem(1.2)is nonsmooth.The(vector)parameterηcompromises the fidelityϕwith the penaltyψ,and its appropriate choice plays a crucial role in obtaining stable yet accurate solutions.Therefore,an automated selection rule and efficient algorithms for determiningηare essential.

One distinct feature of the model(1.2)is that it includes multiple penalties(hence termed as multi-parameter regularization).This is motivated by the following empirical observations.In practice,many objects exhibit multiple distinct features/structures.However,one single penalty generally favors one feature over others,and thus unsuitable for promoting multiple distinct features.For example,total variation(TV)is well suited to reconstructing piecewise constant structures,however,it results in significant staircases in gray regions.One may improve TV-reconstruction by introducing an additional penalty,say L1norm of∆uwhere∆is the Laplacian operator.Hence,a reliable recovery of several distinct features naturally calls for multiple penalties,and it is not surprising that the idea of multi-parameter regularization has been pursued earlier.For instance,in[9]the authors proposed a model to preserve both flat and gray regions in natural images by combining TV with Sobolev smooth penalty.We refer interested readers to[15,17](imaging),[19](microarray data analysis),[18](geodesy)and[13](machine learning)for other interesting applications.

However,a general theory of multi-parameter regularization remains under development[1,4,13,7].In[1]theL-hypersurface was suggested for determining regularization parameters for finite-dimensional linear systems,but without any theoretical justification.In[4],a multi-resolution analysis for ill-posed linear operator equations was analyzed,and some convergence results were established.Lu et al.[13]discussed the discrepancy principle for Hilbert space scales,and derived some error estimates.However,the parameter selection is vastly nonunique due to lack of constraints and thus not directly applicable in practice,for which later a quasi-optimality criterion was suggested(see[14]).Recently,the authors[7]investigated the discrepancy principle and a balancing principle for general convex variational models.However,the nonuniqueness of the discrepancy principle remains unresolved,and further,there is still no theory for the balancing principle for multi-parameter regularization.

The present work extends our earlier work[7],and includes the following essential contributions.We first revisit the balancing principle in[7]from the viewpoint of augmented Tikhonov regularization(see[12]),and establish the equivalence.Then we derive a novel hybrid principle,the balanced discrepancy principle,by incorporating constraints into the augmented approach,which partially resolves the nonuniqueness issue.Further,a priori and a posterior error estimates are derived for both principles.The estimate in Theorem 2.4 was stated in[7]without a proof.Finally,we develop efficient algorithms for implementing these principles,and briefly discuss their properties.

The rest of the paper is organized as follows.In Section 2,we derive the balancing principle and the new hybrid principle,and develop relevant error estimates.In Section 3 we discuss efficient implementations of the two principles.Finally,we provide some numerical results to illustrate the hybrid principle in Section 4.

2 An Augmented Approach

The augmented Tikhonov(a-Tikhonov)regularization is one principled framework for choosing regularization parameters(see[12]).Here we describe the augmented approach for multiparameter models,and derive the balancing principle and a novel balanced discrepancy principle.

2.1 Derivation of the principles

2.1.1 Balancing principle

First we sketch the augmented approach.For the multi-parameter model(1.2),it can be derived analogously from hierarchical Bayesian inference as in[12],and the resulting augmented functionalJ(u,τ,λ)reads

J(u,τ,λ)=τϕ(u,gδ)+λ·ψ(u)+e·(βλ−αlnλ)+β0τ−α0lnτ,

where the vectoreis given bye=(1,1)t.The functionalJ(u,τ,λ)maximizes the posterior probability density function

p(u,τ,λ|gδ)∝p(gδ|u,τ,λ)p(u,τ,λ).

The functionalJ(u,τ,λ)is derived under the assumption that the scalarsλiandτhave Gamma distributions with known parameter pairs.The parameter pairs(α,β)and(α0,β0)are related to the shape parameters in the statistical priors on the prior precisionλiand noise precisionτ,respectively.The special caseβ0=β=0 is known as noninformative prior and customarily adopted in practice.Hence we focus our derivation on this case.Upon lettingthe necessary optimality condition of any minimizer(λi,τ)to the a-Tikhonov functionalJ(u,τ,λ)is given by

Now by rewriting the system with,we arrive at the following system for

The optimality system(2.2)reveals the mechanism of the augmented approach:It selects an optimal regularization parameterηin the model(1.2)by balancing the penaltyψwith the fidelityϕ,from which the term balancing principle follows.We note the term balancing principle here should not be confused with Lepskii’s principle,which is also sometimes called a balancing principle(see[16]).The Lepskii’s principle does require a knowledge of noise level.

Next we characterize(2.2)using the value functionF(η)(see[8])defined by

The functionF(η)is continuous,and it is almost everywhere differentiable,cf.Lemma 2.1.We denote byFηithe partial derivative ofF(η)with respect toηi.The proofis analogous to[8],and hence omitted.

Lemma 2.1The function F(η)is monotone and concave,and hence almost everywhere differentiable.Further,ifit is differentiable,then there holds

Next we provide an alternative characterization of(2.2).First we define the function Φγ(η)by

The necessary optimality condition for Φγ(η),provided thatF(η)is differentiable,reads

which,upon noting Lemma 2.1,is equivalent to

Solving the system with respect toηiyieldsHence,the optimality system of the function Φγcoincides with that of the functionalJ(u,τ,λ).In summary,we have shown our first main result.

Proposition 2.1Let the value function F(η)be differentiable.Then all critical points of the functionΦγare solutions to system(2.2).

Remark 2.1Two remarks on Φγare in order.First,it is very flexible in that the parameterγmay be calibrated to achieve specific desirable properties.Second,by the concavity in Lemma 2.1,F(η)is continuous and thus the problem of minimizing Φγover any bounded and closed region inis well defined.These observations are valid for a general fidelity.

2.1.2 Balanced discrepancy principle

To solve stably and accurately problem(1.1),one should use all prior information,e.g.,the noise levelfor somecm≥1,and other relevant knowledge,whenever it is available.This can be realized by incorporating constraints into the augmented approach,and then deriving the corresponding optimal system.For instance,for the constraintϕ(u,gδ)≤c,the Lagrangian approach gives the following a-Tikhonov functional:

where the unknown scalarµ≥0 is the Lagrange multiplier for the inequality constraintϕ(u,gδ)≤c.Its optimality system reads

Hence the constraint≤cand the balancing principle are both fulfilled:

In the case of one single penalty,identity(2.4)does not provide any additional constraint since the multiplierµis also unknown.We observe that the active constraint,i.e.,is exactly the discrepancy principle(see[5]).The constraint is active under certain conditions(see[10]).Nonetheless,in case of multiple penalties,the discrepancy principle alone cannot uniquely determineη.Hence we include also system(2.4),which might help resolve the nonuniqueness issue.Upon simplification,this yields a new hybrid principle

The principle can be interpreted as the augmented approach with the constraintHence it integrates the classical discrepancy principlecmδwith the balancing principle,and we shall name the new rule(2.5)balanced discrepancy principle.One noteworthy feature of(2.5)is that it does not involve the free parameterγ.

2.2 Error estimates

Now we derive error estimates for(2.3)and(2.5),capitalizing on[3,5–6].We discuss the following three scenarios separately:hybrid principle(2.5),purely balancing principle(2.3)in Hilbert and Banach spaces.These theoretical results partially justify their practical usages.

2.2.1 Balanced discrepancy principle

In this part,we discuss the consistency and an a priori error estimate for the hybrid principle(2.5).To this end,we make the following assumption.

Assumption 2.1There exists a τ-topology such that for anyη>0,the functional Jη(u)is coercive and its level set{u∈C:Jη(u)≤c}for any c>0is compact in τ-topology,and the functionals ϕ and ψiare τ lower semi-continuous.

Remark 2.2Theτ-topology is naturally induced by the penalty functionalψ,and it is not arbitrarily in order to ensure the lower semicontinuity.

Now we can state a consistency result.The line of proofis standard(see[7]),and thus omitted.

Theorem 2.1Let Assumption2.1be fulfilled,andLet the sequence{η(δ)}δbe selected by(2.5).If a subsequence of{η(δ)}δconverges andthen the subsequencecontains a subsequence τ-converging to a·ψ-minimizingsolution of Ku=g†and

Remark 2.3The condition∈(0,1)in Theorem 2.1 is tantamount to the uniform boundedness of the penalties

Next we have the following convergence rate,i.e.,the distance between the approximationand the true solutionu†(in Bregman distance(see[3]))in terms of the noise levelδ.We denote the subdifferential of a convex functionalψ(u)atu†by∂ψ(u†),i.e.,

and the Bregman distancedξ(u,u†)for anyξ∈∂ψ(u†)is defined as

Now we can state a convergence rates result.

Theorem 2.2Let the exact solution u†satisfy the source condition:For any t∈[0,1],there exists a wt∈Y such that

Then for anyη∗determined by the principle(2.5)and with

the following estimate holds:

ProofThe line of proofis again well known,but we include a sketch for completeness.In view of the minimizing property of the approximationand the constraintcmδ,we have

The source condition implies that there exists aξt∗∈∂([t∗,1−t∗]t·ψ(u†))andwt∗∈Ysuch thatFrom this and the Cauchy-Schwarz inequality,we deduce

This shows the desired estimate.

Remark 2.4In Theorem 2.2,the order of convergence relies solely on the constraintwhile the weightt∗is determined by the balancing principle.Hence the reduced system(2.4)does help resolve the vast nonuniqueness issue in the discrepancy principle.

2.2.2 Balancing principle in Hilbert spaces

We derive a posteriori estimates for the balancing principle Φγ(2.3),i.e.,the distance between the approximationand the exact solutionu†in terms of the noise leveland the realized residualWe first treat quadratic regularizationswith linear operatorsLiful filling ker(Li)∩ker(K)={0},i=1,2,and each induces a semi-norm.One typical choice is thatψ1andψ2impose the L2-norm and higher-order Sobolev smoothness,e.g.,We shall utilize a weighted(semi-)norm∥·∥tde fined by

where the weightt≡t(η)∈[0,1]is de fined as before,and byandandClearly,We note that the adjointK∗(and hence)depends on the valuet.

Theorem 2.3Letµ∈(0,1]be fixed,and the exact solution u†satisfy the source condition:For any t∈[0,1],there exists a wt∈Y such thatThen for any parameterη∗selected by(2.3)withthe following estimate holds:

ProofWe decompose the errorand bound the two terms separately.First we estimate the errorIt follows from the optimality conditions foruηandthat

Multiplying the identity withand using the Cauchy-Schwarz and Young’s inequalities give

Next lets=η1+η2.Then we get

Meanwhile,the minimizing property ofη∗to the rule Φγimplies that for any

In particular,we may takeeand arrive at

Next we estimate the approximation erroruη−u†.To this end,we observe

Hence,

Consequently,we deduce from the source condition and the moment inequality(see[5])

where the constantcdepends only on the maximum ofoverFurther,we note the relation

Hence,we deduce

By combining these two estimates,we arrive at the desired inequality.

2.2.3 Balancing principle in Banach space

Lastly,we turn to the balancing principle for general convex regularizationψ.We first recall the following technical lemma(see[11])for single convex regularizationψ.The first estimates the propagation error,and the second plays the role of a triangle inequality.

Lemma 2.2(see[11])Let the exact solution u†satisfy the following source condition:Thereexists a w∈Y such that K∗w=ξ∈∂ψ(u†),and letThen there hold

Now we can state an estimate for the balancing principle(2.3)in Banach spaces.The estimate has been stated in[7]but without a proof.

Theorem 2.4Let the exact solution u†satisfy the source condition:For any t∈[0,1]there exists a wt∈Y such that

Then for everyη∗selected by(2.3)and withthe following estimateholds:

ProofFor anyt∈[0,1],letψt(u†)=[t,1−t]t·ψ(u†)andξt∈∂ψt(u†),withξtandwtbeing the subgradient and the representer in the source condition,respectively.By Lemma 2.2,we have that forη

whereands=η1+η2.It suffices to bound the terms involving Bregman distance.We first estimate the approximation errorTo this end,observe by the minimizing property of the elementuη,i.e.,

This inequality,the definition ofthe source condition and Lemma 2.2 imply

Next we estimate the termIn view of Lemma 2.2,we have

Meanwhile,the minimizing property ofηto the rule Φγgives that for any

Upon lettingand combining the preceding two inequalities,we get

Now combining these three estimates gives the desired assertion.

The a posteriori error estimate in Theorem 2.4 coincides with that for a priori choice,e.g.,η∼δe,if the realized discrepancyδ∗is of the same order with the exact noise levelδ.

3 Numerical Algorithms

Now we describe the algorithms for numerically realizing the hybrid principle and the balancing principle,i.e.,Broyden’s method and fixed-point algorithm,and discuss their properties.

3.1 Broyden’s method

In practice,the application of the hybrid principle invokes solving the nonlinear system(2.5),which is nontrivial due to its potential nonsmoothness and high degree of nonlinearity.We propose using Broyden’s method(see[2])for its efficient solution(see Algorithm 1 for a complete description).

For the numerical treatment,we reformulate system(2.5)equivalently as

The system is numerically more amenable than(2.4).In Algorithm 1,the JacobianJ0can be approximated by finite difference.Step 7 represents the celebrated Broyden update.The stopping criterion is based on monitoring the residual norm∥T(η)∥.Note that each iteration involves evaluatingT(η),which in turn incurs solving one optimization problem of minimizingJη.Our experiences indicate that it converges fast and steadily,however,a convergence analysis is still missing.

Algorithm 1Broyden’s method for system(2.5)

3.2 Fixed point algorithm

In this part,we describe a fixed point algorithm for computing the minimizer of the rule Φγ.The algorithm was originally introduced in[7],but without any analysis.One basic version is listed in Algorithm 2,where the subscript−irefers to the index different fromi.The stopping criterion at Step 4 can be based on monitoring the relative change of the regularization parameterηor the inverse solution

Algorithm 2Fixed point algorithm for minimizing(2.3)

We shall analyze Algorithm 2.First,we introduce a fixed point operatorTby

We shall also need the next result(see[8,Lemma 2.1 and Corrollary 2.3]).

Lemma 3.1The functionis monotonically decreasing in ηi,and the following relations hold:

We have the next monotone result for the fixed point operatorT.

Proposition 3.1Let the function F(η)be twice differentiable.Then the mapT(η)ismonotone if

ProofLetA(η)=ϕ+η2ψ2andB(η)=ϕ+η1ψ1.By Lemma 3.1,there hold

With the help of these two relations,we deduce

and

where we have used the relationfrom Lemma 2.1.Similarly,we have

Therefore,the Jacobian∇Tof the operatorTis given by

Now Lemma 3.1 implies that

Hence,it suffices to show that the determinant|∇T|>0.By Lemma 2.1,the identity

holds,and thus|∇T|is given by

Hence,the nonnegativity of|∇T|follows from the assumption

This concludes the proof.

4 Numerical Experiments

We now provide some numerical results for the hybrid principle(2.5);and the balancing principle(2.3)has been numerically exemplified in[7]and will not be addressed here.The examples are integral equations of the first kind with kernelk(s,t)and solutionu(t).All the examples are taken from[7].The discretized linear system takes the formKu†=g†.The datag†is then corrupted by noises,i.e.,whereζiare standard Gaussian variables,andεis the relative noise level.

4.1 H1-TV model

Example 4.1Letand the kernelk(s,t)is given byξ(s−t).The true solutionu†exhibits both flat and smoothly varying regions and it is shown in Figure 1,and the integration interval is[−6,6].We adopt two penalties

The numerical results are summarized in Table 1.In the table,the subscripts bdp and opt respectively refer to the hybrid principle and the optimal choice,i.e.,the value giving the smallest error.The single-parameter models are indicated by subscripts h1 and tv,and the regularization parameter shown in Table 1 is the optimal one.The accuracy of the results is measured by the relative L2errorWe observe that the H1-TV model inconjunction with the hybrid principle achieves a smaller error than either H1or TV with the optimal choice,thereby showing the advantages of the H1-TV model.Further,the hybrid principle gives an error fairly close to the optimal one,within a factor of two,and the error decreases as the noise level decreases.

Table 1 Numerical results for Example 4.1

Figure 1 Numerical results for Example 4.1 with ε=5%noise.

Let us briefly comment on the performance of the multi-parameter model.The classical H1model recovers the flat region unsatisfactorily,whereas the TV approach clearly suffers from staircasing effect in the gray region and reduced magnitude in the flat region,cf.Figure 1.In contrast,the H1-TV model preserves the magnitude of the flat region while recovering the gray region excellently.Therefore,the H1-TV model does combine the strengths of both H1and TV models.Finally,we would like to remark that Broyden’s method converges rapidly with the convergence achieved in five iterations,and the convergence behavior is not sensitive to the initial guess.

4.2 Elastic-net model

Example 4.2The kernelk(s,t)is given bythe exact solutionu†consists of two bumps and it is shown in Figure 2.The penalties areandto retrieve the groupwise sparsity structure,which is known as elastic-net in statistics(see[19]).The integration interval is[0,1].The size of the problem is 100.

It is observed from Table 2 that the hybrid principle gives slightly too small but otherwise reasonable estimate for the optimal choice.A close look at Figure 2 indicates that the solutionul2has almost no zero entries,and thus it fails to distinguish between relevant and irrelevant factors.Meanwhile,many entries of theℓ1solution are zero,and thus some relevant factors are correctly identified.However,it tends to select only a part instead of all relevant factors.The elastic-net combines the best of bothℓ1andℓ2models,and it achieves the desired goal of identifying the group structure.

Figure 2 Numerical results for Example 4.2 with ε=5%noise.

Table 2 Numerical results for Example 4.2

4.3 Image deblurring

Example 4.3The kernelk(s,t)performs standard Gaussian blur with standard deviation 1 and blurring width 5.The exact solutionu†is shown in Figure 3.The size of the image is 50×50.The penalties are

This example represents a more realistic problem ofimage deblurring.Here one half of the data points are retained,which renders the problem far more ill-posed.Theℓ1solution is very spiky,cf.Figure 3,and neighboring pixels act independently of each other.In particular,many pixels in the blocks and the cross are missing.In contrast,the solutionul2is smooth,but there are many small spurious oscillations in the background.The elastic-net model achieves the best of the two:Retaining the block structure with only few spurious nonzero coefficients.The numbers are also very telling:ebdp=2.96e-1,eo=2.44e-1,el1=9.21e-1,andel2=3.42e-1.Hence,the errorebdpagrees well with the optimal choice,and it is smaller than that with the optimal choice for eitherℓ1orℓ2models.

5 Conclusions

We have studied multi-parameter regularization from the viewpoint of augmented Tikhonov regularization,and shown a unified way to derive the balancing principle and balanced discrepancy principle.A priori and a posteriori error estimates for the principles were provided,and efficient numerical algorithms(Broyden’s method and fixed point algorithm)were presented and discussed.Numerical results were presented to illustrate the feasibility of the balanced discrepancy principle.

Figure 3 Numerical results for Example 4.3 with ε=1%noise.The selected regularization parameters are ηbdp=(4.70e-3,4.65e-3),ηopt=(1.26e-2,1.31e-3),ηl1=5.67e-1,and ηl2=3.51e-3.

AcknowledgementThis work was partially carried out during the visit of the first author at Institute for Applied Mathematics and Computational Science of Texas A&M University.He would like to thank the institute for the hospitality.

[1]Belge,M.,Kilmer,M.E.and Miller,E.L.,Efficient determination of multiple regularization parameters in a generalized L-curve framework,Inverse Problems,18(4),2002,1161–1183.

[2]Broyden,C.G.,A class of methods for solving nonlinear simultaneous equations,Math.Comp.,19(92),1965,577–593.

[3]Burger,M.and Osher,S.,Convergence rates of convex variational regularization,Inverse Problems,20(5),2004,1411–1420.

[4]Chen,Z.,Lu,Y.,Xu,Y.and Yang,H.,Multi-parameter Tikhonov regularization for linear ill-posed operator equations,J.Comput.Math.,26(1),2008,37–55.

[5]Engl,H.W.,Hanke,M.and Neubauer,A.,Regularization ofinverse Problems,Kluwer,Dordrecht,1996.

[6]Hofmann,B.,Kaltenbacher,B.,Poeschl,C.and Scherzer,O.,A convergence rates result for Tikhonov regularization in Banach spaces with non-smooth operators,Inverse Problems,23(3),2007,987–1010.

[7]Ito,K.,Jin,B.and Takeuchi,T.,Multi-parameter Tikhonov regularization,Methods Appl.Anal.,18(1),2011,31–46.

[8]Ito,K.,Jin,B.and Takeuchi,T.,A regularization parameter for nonsmooth Tikhonov regularization,SIAM J.Sci.Comput.,33(3),2011,1415–1438.

[9]Ito,K.and Kunisch,K.,BV-type regularization methods for convoluted objects with edge,flat and grey scales,Inverse Problems,16(4),2000,909–928.

[10]Ivanov,V.K.,Vasin,V.V.and Tanana,V.P.,Theory of Linear Ill-Posed Problems and Its Applications,VSP,Utrecht,2nd edition,2002.

[11]Jin,B.and Lorenz,D.A.,Heuristic parameter-choice rules for convex variational regularization based on error estimates,SIAM J.Numer.Anal.,48(3),2010,1208–1229.

[12]Jin,B.and Zou,J.,Augmented Tikhonov regularization,Inverse Problems,25(2),2009,025001,25 pages.

[13]Lu,S.and Pereverzev,S.V.,Multi-parameter regularization and its numerical regularization,Numer.Math.,118(1),2011,1–31.

[14]Lu,S.,Pereverzev,S.V.,Shao,Y.and Tautenhahn,U.,Discrepancy curves for multi-parameter regularization,J.Inv.Ill-Posed Probl.,18(6),2010,655–676.

[15]Lu,Y.,Shen,L.and Xu,Y.,Multi-parameter regularization methods for high-resolution image reconstruction with displacement errors,IEEE Trans.Circuits Syst.I.Regul.Pap.,54(8),2007,1788–1799.

[16]Math´e,P.,The Lepskii principle revisited,Inverse Problems,22(3),2006,L11–L15.

[17]Stephanakis,I.M.,Regularized image restoration in multiresolution spaces,Opt.Eng.,36(6),1997,1738–1744.

[18]Xu,P.,Fukuda,Y.and Liu,Y.,Multiple parameter regularization:numerical solutions and applications to the determination of geopotential from precise satellite orbits,J.Geod.,80(1),2006,17–27.

[19]Zou,H.and Hastie,T.,Regularization and variable selection via the elastic net,J.R.Stat.Soc.Ser.B,67(2),2005,301–320.