A comparison of parameter choice rules for ℓ p - ℓ q minimization

Images that have been contaminated by various kinds of blur and noise can be restored by the minimization of an ℓ p - ℓ q functional. The quality of the reconstruction depends on the choice of a regularization parameter. Several approaches to determine this parameter have been described in the literature. This work presents a numerical comparison of known approaches as well as of a new one


Introduction
In many areas of Science and Engineering, ranging from medical diagnostics to natural sciences [17,29,34], we are faced with the solution of linear systems of equations of the form where A ∈ R m×n is a large matrix whose singular values decrease to zero with increasing index number without a significant gap.The matrix A then is severely ill-conditioned and may be rank-deficient.In many applications, the right-hand side b δ ∈ R m represents available data that is contaminated by error.The quantity η ∈ R m collects measurement and discretization errors; it is not explicitly known.The vector x ∈ R n represents the signal that we would like to determine.Problems of this kind often are referred to as linear discrete ill-posed problems; see, e.g., [19].We remark that, in imaging applications, as the ones considered here, the unknown x as well as the observed data b δ and the additive noise η are vectorized forms of n 1 × n 2 and m 1 × m 2 2D signals, respectively, with n = n 1 n 2 and m = m 1 m 2 .Let b = b δ − η denote the unknown error-free vector associated with b δ .We are interested in determining the solution x † of the least squares problem min x∈R n ∥Ax−b∥ 2 of minimal Euclidean norm.This solution can be expressed as x † = A † b, where A † denotes the Moore-Penrose pseudo-inverse of A. Since the vector b is not available, it is natural to try to determine the vector x = A † b δ .However, due to the ill-conditioning of A and the error η in b δ , the latter vector often is a meaningless approximation of the desired vector x † .
To compute a meaningful approximation of x † one may resort to regularization methods.These methods replace the ill-posed problem (1) by a nearby well-posed one that is less sensitive to the perturbation in b δ and whose solution is an accurate approximation of x † .Among the various regularization methods described in the literature, the ℓ p -ℓ q minimization method has attracted considerable attention in recent years; see, e.g., [4,12,21,23].This method computes a regularized solution of (1) by solving the minimization problem arg min where 0 < p, q ≤ 2, µ > 0 is a regularization parameter, and the matrix L ∈ R s×n is chosen so that N (A) ∩ N (L) = {0}; here N (M ) denotes the null space of the matrix M and Note that z → ∥z∥ s is not a norm for 0 < s < 1, but for convenience, we nevertheless will refer to this function as a norm also for these values of s.
The first term in (2) ensures that the reconstructed signal fits the measured data, while the second term enforces some a priori information on the reconstruction.The parameter µ balances the two terms and determines the sensitivity of the solution of (2) to the noise in b δ .An ill-suited choice of µ leads to a solution of (2) that is a poor approximation of x † .It therefore is of vital importance to determining a suitable value of µ.
It is the purpose of this paper to review and compare a few popular parameter choice rules that already have been proposed for ℓ p -ℓ q minimization, and to apply, for the first time, a whiteness-based criterion described by Lanza et al. [27].
We briefly comment on the choices of p and q, for which several automatic selection strategies have been developed; see, e.g., [28] and references therein.However, since the focus of this work is the selection of the regularization parameter µ, the parameters p and q are considered to be fixed a priori ; we assume that they have been set to suitable values.
The choice of p should be informed by the statistical properties of the noise.We consider two types of noise, namely noise that can be described by the Generalized Error Distribution (GED) and impulse noise.In the first case the entries of η are independent realization of a random variable with density function where c θ,σ is a constant such that R ξ θ,ν σ (t)dt = 1, ν ∈ R, and θ, σ > 0. For θ = 2, (3) reduces to the Gaussian density function, while for θ = 1 we obtain the Laplace distribution.It is shown in [6] that, for this kind of noise, the maximum a posteriori principle prescribes that p = θ.
The data b δ are said to be affected by impulse noise of level σ if where r j is a uniformly distributed random variable in the dynamic range of b.Numerical experience indicates that it is beneficial to let p < 1 when b δ is contaminated by impulse noise; see, e.g., [6,21,23] for illustrations.The choice of the parameter q is determined by a priori knowledge of the desired solution that we would like to impose on the computed solution.In particular, it is often known that Lx is sparse, i.e., Lx has few nonvanishing entries.This is true, for instance, if L is a discretization of a differential operator or a framelet/wavelet operator.In this case, one ideally would want to let q = 0, where the ℓ 0 -norm of a vector x measures the number of nonzero entries.
However, minimizing the ℓ 0 -norm is an NP-hard problem.Therefore, it is prudent to approximate the ℓ 0 -norm by an ℓ q -norm with 0 < q ≤ 1.For smaller values of q > 0, the ℓ q -norm approximates the ℓ 0 -norm better, however, the minimizing algorithm requires more iterations the smaller q > 0 is, and may suffer from numerical instability for "tiny" q > 0. Therefore, it is usually a good practice to set q small enough, but not too small.We remark that the ℓ q -norm does not satisfy all properties of a norm for 0 < q < 1.This paper is organized as follows: Section 2 outlines two algorithms for computing an approximate solution of (2).In Section 3 we describe the methods for determining µ that are compared in this paper.We report some numerical results in Section 4 and draw conclusions in Section 5.

Majorization-minimization in generalized Krylov subspaces
In [21] the authors proposed an effective iterative method for the solution of (2).At each iteration, a smoothed version of the ℓ p -ℓ q functional, denoted by J ε , is majorized by a quadratic functional that is tangent to J ε at the current iterate x (k) .Then, the quadratic tangent majorant is minimized and the minimizer x (k+1) is the new iterate; see below.Two approaches to determine quadratic majorants are described in [21].We will outline both.
Majorization step.Consider the functional that is minimized in (2).It is shown in [8] that this functional has a global minimizer.When 0 < min{p, q} < 1, the functional ( 4) is neither convex nor differentiable.To construct a quadratic majorant, the functional has to be continuously differentiable.We, therefore, introduce a smoothed version for some ε > 0, where Since Φ s,ε is a differentiable function of t, J ε (x) is everywhere differentiable.We will comment on the choice of ε in Section 4.
We would like to compute When min{p, q} > 1, the functional J ε (x) is strictly convex and therefore the minimization problem (5) has a unique minimum.On the other hand, if min{p, q} < 1, which is the situation of interest to us, then the functional (5) generally is not convex.The methods we describe here therefore determine a stationary point.Let x (k) be the currently available approximation of x * in (5).We first describe the majorant referred to as "adaptive" in [21].This majorant is such that the quadratic approximation of each component of J ε (x) at x (k) has an as large as possible positive second order derivative.In general, each component is approximated by a different quadratic polynomial.We denote the adaptive quadratic tangent majorant of J ε (x) at x (k) by Q A (x, x (k) ).It is characterized by ) is a quadratic functional, (6) where ∇ x denotes the gradient with respect to x.
The functional Q A (x, x (k) ) can be constructed as follows: Evaluate the residual vectors and compute the weight vectors , where all the operations are meant element-wise.The weight vectors determine the diagonal matrices where c ∈ R is a constant that is independent of x; see [21] for details.The new approximation of x * is given by the minimizer x (k+1) of Q A (x, x (k) ).We discuss the computation of an approximation or this minimizer below.
The second type of majorant considered in [21] is referred to as "fixed."This majorant is constructed so that each component of the quadratic polynomial majorants have the same leading coefficient.As shown below, the evaluation of the fixed majorant at x requires less computational work than the computation of the adaptive majorant.However, in general, the method defined by fixed majorants requires more iterations than the method defined by adaptive majorants to reach numerical convergence.
For the fixed case the weight vectors are given by ω and determine the fixed quadratic tangent majorant k) .Here ⟨•, •⟩ denotes the standard inner product and the constant c ∈ R is independent of x.The functional Q F (x, x (k) ) satisfies the properties (6) with Q A (x, x (k) ) replaced by Q F (x, x (k) ); see [21] for details.
Minimization step.We describe how to compute minimizers of Q A and Q F when the matrix A ∈ R m×n is large.To reduce the computational effort, we compute an approximate solution in a subspace V k of fairly small dimension k ≪ min{m, n}.Let the columns of the matrix V k ∈ R n× k form an orthonormal basis for V k .We determine approximations of the minima of Q A and Q F of the form where y (k+1) ∈ R k .We first discuss the adaptive case.We would like to solve min and denote the solution by x (k+1) .This minimization problem is equivalent to Introduce the economic QR factorizations and compute where the superscript T denotes transposition.Assume that This condition typically holds in applications of interest to us.The solution y (k+1) of ( 12) then is unique, and the approximate minimizer of ) is given by (8).
We now enlarge the solution subspace V k by including the normalized residual of the normal equations associated with (9).The residual is given by and the columns of the matrix form an orthonormal basis for the new solution subspace V k+1 .We would like to point out that the vector r (k+1) is proportional to the gradient of k+1) .We refer to the solution subspace V k+1 = range(V k+1 ) as a generalized Krylov subspace.Note that the computation of r (k+1) requires only one matrix-vector product with A T and L T , since we can exploit the QR factorizations (11) and the relation ( 8) to avoid computing any other matrix-vector products with the matrices A and L.Moreover, we store and update the "skinny" matrices AV k and LV k at each iteration to reduce the computational cost.The initial space V 1 is usually chosen to contain a few selected vectors and to be of small dimension.A common choice is We will use this choice in the computed examples reported in Section 4. Summarizing, each iteration of the adaptive approach requires one matrixvector product evaluation with each one of the matrices A, L, A T , and L T , as well as the computation of economic QR factorizations of two tall and skinny matrices, whose column numbers increase by one with each iteration.The latter computations can be quite demanding if the matrices A and L are large and many iterations are required.The algorithm requires the storage of the three matrices V k , AV k , and LV k .In addition, storage of some representations of the matrices A and L is needed.
We turn to the fixed approach.The weight vectors are now given by ( 7), and we would like to solve the minimization problem for x (k+1) , where η = µε q−p .This problem can be expressed as The solution y (k+1) of ( 14) yields the solution x (k+1) = V k y (k+1) of (13).
Introduce the economic QR factorizations Substituting these factorizations into (14), we obtain .
Once we have computed y (k+1) and x (k+1) , we enlarge the solution subspace by including the residual reg of the normal equations associated with (13).Thus, let v new = r (k+1) / r (k+1) 2 .Then the columns of the matrix V k+1 = [V k , v new ] form an orthonormal basis for the solution subspace V k+1 .We remark that the residual is proportional to the gradient of Note that differently from (10), the least-squares problem (14) does not have a diagonal scaling matrix.We therefore may compute the QR factorizations of AV k+1 and LV k+1 by updating the QR factorizations of AV k and LV k , respectively.This reduces the computational work and leads to that each new iteration with the fixed approach is cheaper than with the adaptive approach.Updating formulas for the QR factorization can be found in [13,21].
Each iteration with the fixed approach requires one matrix-vector product evaluation with each one of the matrices A, L, A T , and L T , similarly as for the adaptive approach.Moreover the memory requirements of the fixed and adaptive approaches are essentially the same.
The memory requirement of both the adaptive and fixed approaches outlined grows linearly with the number of iterations.It follows that when the matrix A is large, the memory requirement may be substantial when many iterations are required to satisfy the stopping criterion.This could be a difficulty on computers with fairly little fast memory.Moreover, the arithmetic cost for computing QR factorizations in the adaptive approach and for updating QR factorizations in the fixed approach grows quadratically and linearly, respectively, with the number of iterations.

Parameter choice rules
This section discusses several approaches to determine a suitable value of the regularization parameter µ in (2).The strategies considered in the literature for the choice of the regularization parameter when variational models, including Tikhonov regularization, are employed can be divided into three main classes: (i) Methods relying on the noise level, that either may be known or accurately estimated.A very popular approach belonging to this class is the Discrepancy Principle (DP).(ii) Methods that rely on different statistical and non-deterministic properties of the noise, such as its whiteness.
(iii) Heuristic methods, which are typically only based on the knowledge of the data b δ , among which we mention Generalized Cross Validation (GCV) and the L-curve criterion.
A general discussion on heuristic methods is provided by Kindermann [22].Further references are provided below.In what follows, we are going to review some of the approaches mentioned and their modifications when applied to the solution of (2).Specifically, we are considering stationary and non-stationary scenarios: in the former, the methods are employed a posteriori, i.e., the minimization problem in ( 2) is solved for different µ-values, and the optimal value, µ * , is selected based on a chosen strategy.In the latter scenario, the chosen strategy is applied during the iterations of the generalized Krylov method; therefore the optimization problem ( 2) is solved only once.
We remark that the formulation of the stationary strategies discussed in Section 3.1 does not depend on the selected algorithmic scheme.When considering the non-stationary rules presented in Section 3.2, several issues should be considered, such as the existence of a fast parameter update and the convergence of the overall iterative procedure.In this review, we focus on the robustness of the considered approaches when embedded in the iterations of the generalized Krylov method, and leave a rigorous analysis of the mentioned issues to a future study.

Stationary rules
We describe the stationary rules considered in this paper.

Discrepancy principle
Let the noise that corrupts the data be Gaussian.Then we set p = 2. Let x µ denote the solution of (2) with p = 2, i.e., Assume that a fairly accurate estimate δ > 0 of the norm of the noise is available ∥η∥ 2 ≤ δ.
The discrepancy principle (DP) prescribes that the parameter µ be chosen such that µ DP = sup µ : Ax µ − b δ 2 ≤ τ δ , where τ > 1 is a user-defined constant that is independent of δ.
A non-stationary strategy to estimate µ DP is described in [9], and we will discuss it in Section 3.2.1.Here we propose a stationary way to determine an estimate of the DP parameter.Extensive numerical experience shows that the computation of x µ is fairly stable with respect to the choice of µ and, therefore, a rough estimate of µ DP is usually enough to compute a satisfactory solution.
The analysis of the DP requires that b ∈ R(A); see, e.g., [14].Then µ DP is well defined, since Define the function r(µ) = Ax µ − b δ 2 − τ δ and assume that r(µ) is continuous.This assumption is satisfied for q = 2, see, e.g., [14], but to the best of our knowledge a proof for general q is not currently available.We may employ a root-finder to determine µ DP ; see, e.g., [7,31].

Residual whiteness principle
The whiteness property of the corrupting noise in linear inverse problems of the form (1) has been extensively explored in the context of variational methods.It is convenient to apply the whiteness property since it does not require knowledge of the standard deviation of the noise.Moreover, as it will be made clear in the following, it exploits more information of the data vector b δ than the DP.
The whiteness property has been incorporated in variational models for image denoising and deblurring problems in [3,[24][25][26]32].Despite the high-quality results achieved in these works, these approaches suffer from the strong nonconvexity of the variational models that have to be solved.This makes minimization a very hard task.Other approaches that exploit that the residual image is expected to model white Gaussian noise are described in [18,33] and [1], where the authors propose two statistically-motivated parameter choice procedures based on the maximization of the residual whiteness by the normalized cumulative periodogram and the normalized auto-correlation, respectively.The approach introduced in [1] and applied as an a posteriori parameter choice criterion for image deconvolution problems has been revisited in [27,30].There the authors propose to automatically update the regularization parameter during the iterations of the algorithm used for the minimization of a wide class of convex variational models for image restoration and super-resolution problems.
In what follows, we recall the main steps of the a posteriori criterion described in [1,27] and referred to as the Residual Whiteness Principle (RWP).We remark that, although the RWP has been originally designed for Gaussian noise corruption, it can been applied whenever the noise η in b δ has independent and identically distributed entries.
Consider the noise realization η ∈ R m in (1) represented in its original m 1 × m 2 form.Thus, The sample auto-correlation of η is defined as The scalar components a l,k (η) are given by where the integer pairs (l, k) are referred to as lags, ⋆ and * denote the 2D discrete correlation and convolution operators, respectively, and η ′ (i, j) = η(−i, −j).
Clearly, for ( 16) to be defined for all lags (l, k) ∈ Θ, the noise realization η must be padded with at least m 1 − 1 samples in the vertical direction and m 2 − 1 samples in the horizontal direction.We will assume periodic boundary conditions for η, so that ⋆ and * in ( 16) denote 2D circular correlation and convolution, respectively.Then the auto-correlation has some symmetries that allow us to only consider the lags If the error η in (1) is the realization of a white noise process, then it is well known that the sample auto-correlation a(η) satisfies the asymptotic property: We note that the discrepancy principle relies on exploiting only the lag (0, 0) -among the m asymptotic properties of the noise auto-correlation given in (17).Imposing whiteness of the residual image of the restoration by constraining the residual auto-correlation at non-zero lags to be small is a much stronger requirement.
The whiteness principle can be made independent of the noise level by considering the normalized sample auto-correlation of the noise realization η in (15), namely where ∥η∥ F denotes the Frobenius norm of the matrix η.It follows easily from (17) that We introduce the following σ-independent non-negative scalar measure of whiteness W : R m1×m2 → R + of the noise realization η: Clearly, the nearer the restored image x µ in (2) is to the target uncorrupted image x † , the closer the associated m 1 × m 2 residual image obtained from d µ = Ax µ − b δ is to the white noise realization η in (1) and, hence, the whiter is the residual image according to the scalar measure in (18).The RWP for automatically selecting the regularization parameter µ in variational models of the general form (2) therefore can be formulated as where the scalar cost function W is defined by We refer to W as the residual whiteness function.

Cross Validation
Two heuristic approaches to determine a suitable value of the regularization parameter µ in (2) for any positive values of p and q based on cross validation (CV) are described in [10].This and the following subsections reviews these methods.For details on CV, we refer to [35].Let 1 ≤ d ≪ m and choose d distinct random integers i 1 , . . ., i d in {1, . . ., m}.Remove rows i 1 , . . ., i d from A and b δ .This gives the "reduced" matrix and data vector A ∈ R (m−d)×n and b δ ∈ R m−d , respectively.Consider a set {µ 1 , . . ., µ l } of positive parameter values.For each µ j , we solve Our aim is to determine the parameter µ j that yields a solution that predicts the values b δ i1 , . . ., b δ i d well.Therefore, we compute for each j the residual error and let j * = arg min j r j .

Modified Cross Validation
Cross validation can be applied to predict entries of the computed solution instead of entries of the data vector b δ .This is described in [10] and there referred to as Modified Cross Validation (MCV).We outline this approach.The idea is to determine a value of µ such that the computed solution is stable with respect to loss of data.Let 1 ≤ d ≪ m and select two sets of d distinct indices between 1 and m referred to as 2 , . . ., i d } and d }.Let A 1 and b δ 1 , and A 2 and b δ 2 , denote the matrices and data vectors obtained by removing the rows with indexes I 1 and I 2 from A and b δ , respectively.Let {µ 1 , . . ., µ l } be a set of regularization parameters and solve We then compute the quantities and define µ * = µ j * .To reduce variability, the computations are repeated K times.This results in K values of µ * , which we denote by µ * 1 , . . ., µ * K .The computed approximation of x † is given by arg min

Non-stationary rules
We now describe a few non-stationary rules for determining the regularization parameter.

Discrepancy principle
An iterative, nonstationary, approach to determine the regularization parameter µ for the minimization problem (2) when p = 2 and the error η in b δ is white Gaussian is described in [9].We review this method for the fixed approach described at the end of Section 2. Our discussion easily can be extended to the adaptive case.

Residual whiteness principle
The iterated version of the RWP outlined in Section 3.1.2was proposed in [27] for convex variational models, which can be solved by the Alternating Direction Method of Multipliers (ADMM).The approach described in [27] for the situation when the noise η is white Gaussian relies on the following steps: 1. Find an explicit expression of the residual whiteness function W defined in (19) in terms of µ for a general ℓ 2 -ℓ 2 variational model.2. Exploit this expression of W as a function of µ to automatically update the regularization parameter µ in the ℓ 2 -ℓ 2 subproblem arising when employing ADMM using a suitable variable splitting.
The iterative procedure proposed in [27] can not be directly applied to the generalized Krylov method considered here.However, to explore the potential of the RWP when applied during the iterations of a generalized Krylov scheme, one can minimize the residual whiteness measure in (19) at each iteration in a similar fashion as for the DP.
For simplicity let us consider the adaptive approach, and note that the extension to the fixed case is straightforward.At each iteration, we solve Denote by .
The iterated RWP thus reduces to finding the µ-value by solving the minimization problem with The minimization problem ( 22) is solved by means of a Matlab optimization routine that relies on the existence of a closed form for y µ .Alternatively, one could compute y µ on a grid of different µ-values and compute the corresponding whiteness measure function.Nonetheless, we believe that the design of an efficient algorithmic procedure for tackling problem ( 22) is worth further investigation.
We remark that the outlined strategy defines a new non-stationary method, for which a convergence proof is not yet available.

Generalized Cross Validation
This subsection summarizes the approach presented in [11].Similarly as above, we consider a non-stationary method for determining µ.This method can be applied to any kind of noise.
Let us first assume that the noise is Gaussian and recall Generalized Cross Validation (GCV) for Tikhonov regularization, i.e., when p = q = 2 in (2): Define the functional The regularization parameter determined by the GCV method is given by When the GSVD of the matrix pair {A, L} is available, G(µ) can be evaluated inexpensively.If the matrices A and L are large, then they can be reduced to smaller matrices by a Krylov-type method and the GCV method can be applied to the reduced problem so obtained; see [11] for more details.This is how we apply the GCV method in the methods of the present paper.Specifically, in the adaptive solution method, we have to solve the minimization problem (12), which is of the same form as (23).Therefore, one can use GCV to compute an appropriate value of µ.Since the matrices R A and R L are fairly small, it is possible to compute the GSVD quite cheaply.It follows that the parameter µ GCV can be determined fairly inexpensively.Moreover, the projection into a generalized Krylov subspace accentuates the convexity of G, making minimization easier; see [15,16] discussions.
We iterate the above approach similarly as we iterated the discrepancy principle, and compute the parameter µ GCV in each iteration.This furnishes a non-stationary algorithm.In detail, consider the kth iteration with regularization parameter µ (k) , Compute the GSVD of the matrix pair {R A , R L }.It gives the factorizations where Σ A and Σ L are diagonal matrices, the matrix X is nonsingular, and the matrices U and V have orthonormal columns.We would like to compute the minimizer of where (24), we get Since Σ A and Σ L are diagonal matrices, the value of G(µ) can be computed cheaply for any value of µ; see, e.g., [5] for a derivation.We let at each iteration Note that, even though it does not appear explicitly in the formulas, the function G varies with k.
We now consider the case in which the noise is not Gaussian.Since, in the continuous setting, the value of Ax and the parameter µ is determined by minimizing G smooth (µ).The non-stationary ℓ p -ℓ q method in this case is obtained in a similar fashion as above.We therefore do not dwell on the details.In the following, we will refer to this method as the GCV method regardless of the type of noise, and use G for Gaussian noise and G smooth for non-Gaussian noise.

Numerical experiments
We compare the parameter choice rules described above when applied to some image deblurring problems.These problems can be modeled by a Fredholm integral equation of the first kind where the function g represents the blurred image, k is a possibly smooth integral kernel with compact support, and f is the sharp image that we would like to recover.Since k has compact support and is smooth, the solution of ( 25) is an ill-posed problem.When the blur is spatially invariant, the integral equation ( 25) reduces to a convolution The kernel k is often referred to as the Point Spread Function (PSF).After discretization we obtain a linear system of equation of the form (1). Since we only have access to a limited field of view (FOV), it is customary to make an assumption about the behavior of f outside the FOV, i.e., we impose boundary conditions to the problem; see [20] for more details on image deblurring.We consider three different images and kinds of blur, and three different types of noise.Specifically, we regard the cameraman image (242 × 242) and blur it with a motion PSF, the boat image (222 × 222), which we blur with a PSF that simulates the effect of taking a picture while one's hands are shaking, and the clock image (226×226) with a Gaussian PSF; see Figure 1.We consider Gaussian noise, Laplace noise, and a mixture of impulse and Gaussian noise.In the first case, we scale the noise so that ∥η∥ 2 = 0.02 ∥b∥ 2 , the second case is obtained by setting θ = 1 and σ = 5 in (3), and in the third case we first modify randomly 20% of the pixels of b and then add white Gaussian noise so that ∥η∥ 2 = 0.01 ∥b∥ 2 .Following [6], we set p = 2 when the data is corrupted by Gaussian noise, p = 1 when the data is contaminated by Laplace noise, and p = 0.8 in the mixed noise case.We let L be a discretization of the gradient operator.Assume, for simplicity, that n = n 2 1 and define the matrix where I n1 denotes the n 1 × n 1 identity matrix and ⊗ is the Kronecker product.Natural images typically have a sparse gradient.Therefore, we set q = 0.1.We now briefly discuss the computational effort required by each parameter choice rule.The cost of the stationary methods depends on how many values of µ and, for CV and MCV, on how many training and testing sets are considered.For DP and RWP we sample 15 values of µ, while for CV and MCV we consider 10 values of µ and 10 different training sets.Therefore, DP and RWP require the MM algorithm to be run 15 time, while CV and MCV run the MM algorithm 101 and 201 times, respectively.The non-stationary methods require a single run, however, the regularization parameter µ k has to be tuned at each iteration.Therefore, the cost of a single run of a non-stationary method is, in general, more expensive that the cost of single run of a stationary method.Nevertheless, since the computation of the parameter µ k can be performed cheaply, thanks to the projection into the generalized Krylov subspace, the non-stationary methods are overall more computationally efficient than their stationary counterparts.We compare the performances of the considered parameter choice rules in term of accuracy using the Relative Restoration Error (RRE), defined as Table 1 displays values of the RRE.We report the optimal RRE as well, i.e., the RRE obtained by hand-tuning µ to minimize the RRE.For each test image and noise type, the RRE value that is closest to the optimal one is reported in bold.
Since the DP only can be used for Gaussian noise, we cannot determine the DP parameter for Laplace and mixed noise, while the RWP can be used only for white noise.Table 1 shows the DP and RWP to usually provide very accurate reconstructions, that is the achieved RRE is very close to the optimal one.The limitation of the DP is that it requires a fairly accurate estimate of ∥η∥ 2 to be known.However, due to the Bakushinskii veto [2], the DP method is the only one for which a complete theoretical analysis is possible.
The CV, MCV, and GCV methods can be applied to any type of noise; however, being so-called heuristic methods, they may fail in certain situations.Table 1 indicates that the MCV tends to select a µ-value that in some cases leads to significantly large RREs.The same behavior can be observed for the other CV-based strategies in the mixed noise case.Nonetheless, in general the CV algorithm provides the most consistent results.
For the strategies that admit both stationary and non-stationary formulations, the DP for Gaussian noise and the RWP for Gaussian and Laplace noises, we observe that for the two scenarios the corresponding RRE are very close, suggesting that the computed restorations are very similar.This behavior, which can be rigorously justified when the DP is used, provides empirical evidences of the robustness of the RWP.
We report in Figure 2 the reconstructions obtained with the optimal choice of the parameter µ in all considered cases.Note that, if the parameter µ is chosen properly, then ℓ p -ℓ q minimization is able to determine very accurate reconstructions.

Conclusions
This paper compares a few parameter choice rules for the ℓ p -ℓ q minimization method.Their pros and cons are discussed and their performances are illustrated.We have shown that, if the regularization parameter is tuned carefully, then the ℓ p -ℓ q model, solved by means of the generalized Krylov method, can provide very accurate reconstructions.The RWP, which here has been applied with the ℓ p -ℓ q model for the first time, can be seen to be particularly robust and is able to determine restorations of high quality.
µ − b δ 2 is infinite when b δ is corrupted by impulse noise, a smoothed version of GCV was proposed in [11].Let b δ smooth denote a smoothed version of b δ obtained by convolving b δ with a Gaussian kernel; see [11] for details.Consider the smoothed function

Table 1 :
RRE obtained with the different parameter choice rules in the computed examples.