{'fraction_non_alphanumeric': 0.04200962390781309, 'fraction_numerical': 0.012330631885526149, 'mean_word_length': 4.135506421720046, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 29, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Learning Theory Approach to a Computationally Efficient Parameter Selection for the Elastic Net

Despite recent advances in regularisation theory, the issue of parameter selection still remains a challenge for most applications. In a recent work the framework of statistical learning was used to approximate the optimal Tikhonov regularisation parameter from noisy data. In this work, we improve their results and extend the analysis to the elastic net regularisation, providing explicit error bounds on the accuracy of the approximated parameter and the corresponding regularisation solution in a simplified case. Furthermore, in the general case we design a data-driven, automated algorithm for the computation of an approximate regularisation parameter. Our analysis combines statistical learning theory with insights from regularisation theory. We compare our approach with state-of-the-art parameter selection criteria and illustrate its superiority in terms of accuracy and computational time on simulated and real data sets.


Introduction
Many applications deal with the problem of the recovery of an unknown quantity of interest x ∈ R d from its corrupted observations y ∈ R m . In most cases the relationship between x and y can be approximately described by the following inverse problem y = Ax + σw, where A ∈ R m×d is a known matrix describing the measuring apparatus, w is a zero-mean isotropic random vector, modelling the noise, and σ > 0 is the noise level. Inverse problems of this type are ubiquitous in image processing tasks, such as: denoising, where A is an identity operator; deblurring, where A is a convolution operator; and inpainting, where A is a masking operator.
The recovery of the original signal x from the corrupted observation y is an ill-posed inverse problem. Thus, theory of inverse problems suggests the use of suitable regularization techniques [18]. Specifically, we are interested in approximating x with the minimizer of the following functional argmin z∈R d where · 2 is the Euclidean norm modelling data-fidelity, J is a penalty term encoding an a priori knowledge on the nature of the true solution, and λ is a regularization parameter determining a trade-off between these two terms. Having the penalty term fixed, a central issue concerns selection of λ. The optimal parameter λ is the one that minimizes the discrepancy between the minimizer z λ of (1) and the exact solution x λ opt = argmin Unfortunately, x is unknown and there is often no information regarding either the noise w or the noise level σ. Hence, λ opt needs to be approximated. Choosing a good approximation to λ opt is a non-trivial, problem-dependent task and has motivated significant research on parameter selection over the last decades. However, there is still no universal framework for fast and efficient data-driven parameter selection. In this paper, we aim at (partially) closing this gap and provide a novel concept for automated parameter selection by recasting the problem to the framework of statistical learning theory. Specifically, inspired by very recent and (to our knowledge) first results in this spirit for Tikhonov regularization [11], we propose a method for learning the optimal parameter for elastic net regularization [39] from a dataset of corrupted observations. Existing parameter selection methods. Parameter selection rules used in theory and in practice can be broadly classified as those based on the discrepancy principle ( [1,28] and references therein), generalized cross-validation (GCV) [22], balancing principle [26,36], quasi-optimality [24,33] and various estimations of the mean-squared error (MSE) (see [12,29] and references therein). GCV is a popular parameter selection rule for linear methods since it gives a closed form of the regularization parameter and does not require tuning of any additional parameters or the knowledge of the noise. Follow-up studies [38] extended GCV to nonlinear estimators, but those cover specialized cases and do not provide a closed form or implicit equations for computing λ opt . Balancing principle has received a lot of attention in the inverse problems community and has also been studied in the framework of learning theory. Quasi-optimality is one of the simplest available parameter choice methods. Though it is not as stable as the balancing principle, it does not require any information about the problem. Discrepancy and MSE-based principles still remain the preferred methods for parameter selection for nonlinear estimators due to their simplicity and accuracy. We refer to a recent rather comprehensive comparative study on the existing approaches [3].
In order to select the regularization parameter most existing methods require the regularized solution z λ to be computed for various λ, i.e., over a predefined grid of values. The 'optimal' regularization parameters are then chosen according to some criteria, e.g., loss over a validation set. To find regularization parameters by an exhaustive search is a computationally expensive task, especially in the high-dimensional data scenario, with often no guarantees on the quality of approximation. Moreover, selection criteria often presuppose that some a priori information is readily available, such as an accurate estimate of the noise level (in e.g. discrepancy principle) or bounds on the noise error (in e.g. balancing principle). Furthermore, most parameter selection methods require additional, method-specific parameters to be preselected. Considering these issues it is fair to say that state-of-the-art methods usually require tedious manual adjustments before they can be utilized. The situation becomes even more challenging when the penalty term J takes a more complex form, such as the elastic net or total variation (as commonly used for image processing tasks). The question of how to obtain an accurate regularized solution z λ , with a near-optimal parameter λ, while ensuring low computational complexity and minimizing the need for manual intervention is the main motivation of our work. In particular, we propose to recast the task of parameter selection to the framework of statistical learning, where one is interested in learning a function R : y → λ opt , from a training set of corrupted data, such that λ opt is a good approximation of the optimal parameter λ opt .
Elastic net regularization. Inspired by recent results in this spirit [11], in this paper we propose a method for learning the optimal parameter for elastic net regularization [39] which is efined by where α ≥ 0 is a hyperparameter controlling the trade-off between the 1 (Lasso) and the 2 (Tikhonov) penalty terms. Our main motivation for considering the elastic net is that they have been shown to produce sparse models comparable to 1 , while often achieving superior accuracy in real-world and simulated data. The elastic net overcomes main limitations of 1 minimization, namely, it encourages the grouping effect, which is relevant for many real-life applications such as microarray classification and face detection, see [10] and references therein.
In the following, we summarize the main algorithmic and theoretical developments related to the elastic net, which was first proposed by Zou and Hastie [39]. Therein the authors rewrite the elastic net functional as Lasso regularization with augmented datum. Then, the authors apply LARS [16] for the effective reconstruction of the entire solution path, and use cross-validation for the selection of the optimal regularization parameter. Later work [10] studies the theoretical properties of (3) in the context of learning theory, analyses the underlying functional and uses iterative soft-thresholding to compute the solution. Moreover, as the parameter choice rule authors provide an adaptive version of the balancing principle [26,36]. The rule aims to balance approximation and sample errors, which have opposite behaviour with respect to the tuning parameter. At the same time, the proposed rule requires multiple estimations of z λ for various values of λ, which makes it less applicable for high-dimensional tasks. We will rework some of the arguments from [10] for the computation of z λ , while keeping our focus on an efficient approach for parameter learning. In [25] the authors propose an active set algorithm for solving (3). Addressing the problem in the framework of classical regularization theory, the authors consider the discrepancy principle [6,28] for determining the parameter. As for the balancing principle, implementation of the discrepancy principle requires multiple estimations of the solution z λ for different parameter values, and a pre-tuning of other, method-specific parameters. Moreover, it is assumed that the noise level is known, which is often not the case in practice.
Paper [27] proposes a hybrid method for tuning parameters λ and α in a model fitting problem y i = a i x, i = 1, . . . , n, where y i ∈ R, a i ∈ R p and x ∈ R p . The authors propose an alternating method which updates the solution z λ using coordinate descent and then updates λ and α in one iteration. The main advantage is its efficiency, as one does not need to calculate z λ for multiple parameters at once, but rather calculate the solution on a much coarser parameter grid. The method is in spirit similar to the LARS method, but has better scalability. At the same time, it requires a non-convex problem to be solved, and hence has inherent limitations. Moreover, this method is not directly applicable for our inverse-problem setting, where we do not observe changes in the same solution x, but rather assume that the design matrix A is fixed and each response y is generated by a new original signal x. In summary, we are not aware of any parameter choice rule for elastic net that allows to select a parameter without any a priori assumptions and manual adjustments.
This paper leverages the work [11] where the parameter selection is considered in the context of nonparametric regression with random design. In particular, the authors propose a fully data-driven method for the determination of the optimal parameter for Tikhonov regularization under the assumption that a training set of independent observations y 1 , . . . , y N , is made available, each of them associated with an (unknown) signal x 1 , . . . , x N , through y i = Ax i + σw i . The goal is to find a prediction function R : y → λ opt based on the given training set, which provides a good approximation to the optimal parameter λ opt for a new datum y. The starting point of the method is to find an empirical proxy x to the real solution x by assuming that x 1 , . . . , x N are distributed over a lower-dimensional linear subspace. The authors show that x = A † Πy, where A † is the pseudo-inverse of A and Π is the projection onto a suitable eigenspace of the empirical covariance matrix, is a good proxy to x provided N is large enough. Then, using the proxy, one can select the regularization parameter as the minimizer of where z λ Tik is the minimizer of the Tikhonov functional min The analysis and the techniques related to the construction of x are independent of the choice of the optimisation scheme, whereas the selection of λ opt is defined by the regularization scheme. However, it is worthwhile to mention that if A † is not injective, A † Πy is not a good proxy of x. For Tikhonov regularization this is not an issue as, without loss of generality, we can always assume that A is injective. Specifically, one can replace x in (2) with x † = A † Ax and recall that z λ Tik belongs to ker A ⊥ for all λ. Therefore, for a wider applicability of the suggested framework, it is important to address the construction of x, and the selection of λ opt for a wider class of problems.
In this paper, we extend the framework of [11] by providing the analysis and algorithms to the elastic net, which is a non-linear method, and to non-injective operators. Moreover, we improve some of the theoretical results, and develop an efficient, fully automated algorithm. To do this, we analyse our problem in two settings: (i) simplified case: for the case A = Id (corresponding to image denoising) and w when are independent Bernoulli random vectors. We provide a bound on |λ opt − λ opt | with high probability and discuss the number of samples required for an optimal learning, see Proposition 3.4. We also consider the case of a bounded w, which motivates our algorithm. Though the model might be oversimplified, it captures the essence of the problem and our experiments confirm the results for more general settings.
(ii) general case: for a general matrix A, we provide an efficient and accurate algorithm based on gradientdescent for the computation of an approximate optimal parameter. We illustrate the performance of our algorithm comparing it to state-of-the-art parameter choice methods on a number of synthetic and imaging problems. The obtained results show that our approach has superior accuracy and often outperforms other methods from a computational perspective.

Outline
In Section 2, we describe the main ingredients for our approach. Therein we define and prove bounds on our empirical estimators, discuss minimizers of the elastic net (3), and define loss functions that will be used for approximately optimal parameter selection. Section 3 provides the main theoretical results of the paper, showing that in a simplified setting the empirical estimator of the optimal parameter is indeed close to the true one. In Section 4 we present an efficient and accurate algorithm for the computation of an approximate optimal parameter. We discuss the performance of the parameter learning method by means of several numerical experiments on synthetic and imaging data in Section 5. Therein, we compare our method with state-of-the-art parameter selection criteria in terms of accuracy of the solution recovery, closeness to the optimal parameter, and computational time. Our focus on imaging data is on wavelet denoising, where we work on synthetic data sets, as well as real-world brain MRI data set with heteroscedastic noise. Finally, we conclude with a brief discussion about future directions in Section 6. The Appendix contains proofs of some auxiliary, and technical results.
We denote by {e i } i=1,...,d the canonical basis of R d . For a matrix M, we denote its transpose by M , its Moore-Penrose pseudo inverse by M † , and its spectral norm by M 2 . Furthermore, range(M) and ker(M) are the range and the null space of M, respectively. For a square-matrix M, we use trace(M) to denote its trace. The identity matrix is denoted by Id and we use 1 D to denote the indicator function of a set D ⊂ R d . For The value ξ ψ 2 is called the sub-Gaussian norm of ξ, with which the space of sub-Gaussian vectors becomes a normed vector space [35]. The (non-centred) covariance of a random vector ξ is denoted as Finally, we write a b if there exists an absolute constant C > 0 such that a ≤ Cb.

Problem setting
We consider the following stochastic linear inverse problem: given a deterministic matrix A ∈ R m×d , we are interested in recovering a vector x ∈ R d from observations where (A1) the unknown datum x ∈ R d is a sub-Gaussian vector, such that x ψ 2 = 1; (A3) the noise w ∈ R m is an independent sub-Gaussian vector, such that σw ψ 2 ≤ 1, Σ(w) = Id and σ > 0 is the noise level.
Conditions (A1) and (A3) are standard assumptions on the distributions of the exact datum x and the noise σw ensuring that the tails have fast decay, and note that normalisation conditions on x and w can always be satisfied by rescaling. Furthermore, it follows from the definitions that V is the smallest subspace such that x ∈ V almost surely. Thus, by (A2), the exact datum x is h-sparse almost surely and, since ker (A) ∩ V = {0}, it is a unique vector with such a property. Define now W = range (Σ(Ax)) . The following simple result for an injective A was shown in [11]; here we extend it to the general case.

Proof. A direct computation gives
where P denotes the orthogonal projection onto V. Assumption (A2) says that A is injective on V, and thus Σ(Ax) and Σ(x) have the same rank h. Furthermore (AP) maps R d onto V, so that range(Σ(Ax)) = (APΣ(x))V = APV = W, where Σ(x)V = V, since Σ(x) is symmetric.

Empirical estimators
Lemma 2.1 suggests that V could be directly recovered if A were invertible and W were known. In most practical situations though, neither of those assumptions is satisfied: we only have access to noisy observations and A could be not only non-invertible, but also non-injective. We will address this issue by recasting the entire problem to a statistical learning framework, similar to [11]. Namely, suppose we are given observation samples y 1 , . . . , y N such that y i = Ax i + σw i for i = 1, . . . , N, where x i , w i and σ are unknown, and let be the empirical covariance of the observations. Standard statistical learning theory suggests that Σ (y) is a good approximation to Σ (y) provided N is large enough. As a consequence, we will show that the vector space W, which is spanned by the first h eigenvectors of Σ (y), is a good estimator of W. To justify the above claims, observe first that since Σ(w) = Id holds by (A3), we have Therefore, Σ (y) and Σ (Ax) have the same eigenvectors and the spectrum of Σ (y) is just a shift of the spectrum of Σ (Ax) by σ 2 . Let λ 1 ≥ . . . ≥ λ h be the non-zero eigenvalues of Σ(Ax), counting for multiplicity, and α 1 ≥ . . . ≥ α m and α 1 ≥ . . . ≥ α m be the eigenvalues of Σ(y) and Σ(y), respectively. From (6) it thus follows Let Π be the (orthogonal) projection onto W, which has rank h due to Lemma 2.1, and let Π be the (orthogonal) projection onto W. Using matrix concentration inequalities we will now show the fundamental tool of our study: that Π is an accurate and an unbiased approximation of Π. We distinguish between cases of bounded and unbounded y and improve upon the results in [11].
Assume now that y 2 ≤ √ L holds almost surely and consider a family of independent m × h matrices Applying the matrix Bernstein inequality for rectangular matrices (Theorem 6.1.1. in [34]) we have for all u > 0 A direction computation gives M ≤ L Σ(y) 2 . It follows that holds with probability greater than 1 − exp(−u) for every u > 0. Moreover, by analogous argumentation (12) holds provided N ≥ (log(2m) + u), see A.2 in the Appendix for details. Thus, (9) follows by applying (16).
The previous result comes with a certain caveat. Namely, it is of only theoretical value as the proof implicitly assumes that either h or the spectral gap are known (which informs the choice of the approximate projector Π). In practice, however a proper eigenspace can only be detected if there is a spectral gap, and if it occurs around the eigenvalue corresponding to the eigenspace we want to recover, i.e., if λ h > δ, where Indeed, under this assumption the empirical covariance matrix also has a spectral gap around the desired eigenvalue, as shown by the next result.
Proposition 2.3. Assume (17) holds. Then the empirical covariance matrix has a spectral gap at the h-th eigenvalue, with probability greater than 1 − 2 exp(−u), provided δ < λ h and N , by adding and subtracting α j and α j+1 inside the first term.
and thus, h will be wrongly estimated. In many practical situations though, the situation is not as pessimistic as this result would suggest. In other words, estimating h could be relatively easy when the singular values are properly visualized and interpreted, even though applying a rule as direct as the spectral gap would lead to an unwanted conclusion. We devote more attention to this question in Section 5.1, and suggest alternative heuristic means of estimating the intrinsic dimension h.
We are now ready to define our empirical estimator of x. Define Q = AA † . Note that Q is the orthogonal projector onto range(A) = ker ⊥ (A ), whereas P = A † A is the orthogonal projector onto range(A ) = ker ⊥ (A). The empirical estimator of x is defined as For η = y − Πy the empirical estimator x satisfies the (empirical) inverse problem A x + Q η = Qy. (19) In the following, we use x to learn the optimal regularization parameter for elastic net minimization, though this part of the analysis is independent of the choice of the optimization scheme.
Remark 2.4. One might think of using x as a sought solution and then completely avoiding regularization and, thus, the issue of the parameter choice. This has been partially discussed in [11] and the authors show that in some cases, esp. when the training set size N is small or noise level is small, x − x is larger than z λ − x . Moreover, as we also show in our numerical experiments, x is not a good proxy to x when A is not injective. In addition, x does not preserve the structure of the original signal, e.g., x is, in general, not sparse, whereas x is sparse. Therefore, we remind that we are interested in using an estimator x for which λ opt − λ opt is small with high probability, and we are not interested in directly controlling x − x , which is the goal in manifold learning [4].
To begin, we observe that minimizers of the empirical problem coincide with minimizers of the original problem.
Since the second term does not depend on z we get argmin as desired.

Elastic net minimisation
From now on we focus on the parameter choice for the elastic net minimization, where The first regularization term, z 1 , enforces the sparsity of the solution, whereas the z 2 2 terms gives smoothness and ensures that in case of highly correlated features and we can correctly retrieve the whole relevant group (which is one of the limitations of LASSO minimization). We first recall some basic facts about existence, uniqueness and sensitivity of solutions with respect to regularization parameters [25].
Lemma 2.6. The functional J is strictly convex and coercive (i.e. J(z) ≥ C z 2 2 , for an absolute constant C > 0). Moreover, for each λ > 0 the minimizer of (1) exists and is unique and the mapping λ → z λ is continuous with respect to λ > 0.
In the remainder of this paper, we redefine the elastic net regularized solution by recasting λ to a bounded interval as z t (y) = argmin where t ∈ [0, 1] and α > 0 is a fixed parameter. For t ∈ (0, 1) the solutions of (21) correspond to (20) for λ = 1−t t . On the other hand, for t = 0 we get z 0 (y) = 0, and for t = 1 we define z 1 (y) := x α , where This definition is driven by the following observations. First, the set N = A † y ⊕ ker (A) is non-empty (since A has finite rank). Furthermore, it was shown in [10] and [25] that in case of elastic nets minimization In other words, vector x α plays the role of the Moore-Penrose solution x † in linear regularization schemes [18]. By Lemma 2.6 the minimizer of (21) always exist and is unique, the map t → z t is continuous for t ∈ (0, 1). Equation (23) implies that t → z t is continuous at t = 1, and later in (32) we show that the continuity also holds at t = 0.

Quadratic loss functionals
To drive the selection of regularization parameters we go back to the first principles and consider the following quadratic loss functionals.
are called the true and the empirical quadratic loss, respectively. Furthermore, are the true and the empirical optimal regularization parameters.
In view of Lemma 2.6 and the subsequent discussion, the decision to cast the regularization parameter from the set of all positive reals to [0, 1], has a clear purpose: true and empirical loss functions are both continuous, defined on a bounded interval, and, hence, achieve a minimum. Since x is unknown, R(t) is only of theoretical value. Thus, our aim is to minimize R while ensuring t opt − t opt is small.
Before going into details, let us briefly discuss some difficulties associated with elastic net minimization which need to be addressed. The elastic net solution admits a closed form solution only when A A = Id and, in other cases, the solution can only be approximated. Furthermore, as we will see below, loss functionals R, R are not globally differentiable, but rather only piecewise. The unavailability of a closed form of z t (y) and the non-smoothness of loss functionals means that the parameter in general cannot be analysed in full detail. Therefore, in the following we split the analysis into two parts: a simplified case where for A A = Id we show that t opt − t opt is small, and the general case where for any A we provide a simple and efficient algorithm for parameter learning. Furthermore, we need to amend the empirical loss function (24) in the case when measurement matrix A is non injective since R(t) cannot be reliably estimated for non linear methods, see Figure 1 and [17]. We follow the idea of SURE-based methods [21], which provide an unbiased estimate of R(t) by projecting the regularized solution onto ker ⊥ (A). In the left panel A is injective and t opt is a good proxy for t opt . In the middle panel rank(A) = 40 and we see that t opt does not approximate t opt well. On the other hand, the right panel shows that in case of a non-injective matrix, modified and projected losses behave really well where P = A † A is the orthogonal projection onto ker ⊥ (A), is called the projected empirical loss. Furthermore, is the projected empirical optimal regularization parameters.
Moreover, using the fact that x = A † Πy, we define the modified projected risk.
is called the modified projected risk. Furthermore, is the modified empirical optimal regularization parameters.
By considering R M instead of R P , we avoid the computation of the Moore-Penrose inverse A † , which might be either costly to compute or indeed numerically unstable if A is poorly conditioned. As we will show in Section 5 and can be seen in the right-most panel in Figure 1, using R P and R M instead of R when A is non-injective leads to a drastically improved performance. Projecting onto ker ⊥ (A) does affect the fundamental shape of the functional by making it somewhat smoother (decreasing the gradients), which we will need to take into account in the numerical part of our study.

Solution via soft-thresholding
As mentioned earlier, elastic net minimization does not admit a closed form solution in case of a general forward matrix A. In the original paper by Zou and Hastie [39], the elastic net problem is recast as a Lasso problem with augmented data, which can then be solved by many different algorithms (e.g. the LARS method [15]). Alternative algorithms compute the elastic net minimizer directly, and are generally either of the active set [25] or the iterative soft-thresholding-type [10]. Here we adhere to iterative soft-thresholding, and rework the arguments in [10] to show that the solution to (21) can be obtained through fixed point iterations for all t ∈ [0, 1]. To begin, define the soft-thresholding function by and the corresponding soft-thresholding operator S τ (u), acting component-wise on vectors u ∈ R m . The following lemma shows that the elastic net solution is a fixed point solution of a contractive map. The proof is included in the Appendix for the sake of completeness.
with Lipschitz constant where σ m and σ M are the smallest and the largest singular values of the matrix A, respectively.
For t = 0, the solution is z 0 = 0, which is consistent with (31). Furthermore, by (30) and (31), we get Our definition of the solution at t = 1, in (22), also satisfies z = T t (z) since S 1−t is identity and, thus, z t=1 (y) = A † y, though the map T 1 is not a contraction. In summary, the solutions are consistent with the statement of Lemma 2.6, as expected.

Closed form solution
As shown in [39], in the case of orthogonal design, i.e. A A = Id, the solution of the elastic net minimization (1) is given by Plugging λ = 1−t t into (33) we have

Parameter Error
Since the elastic net solution is not given in closed form for a general A, a rigorous study of the parameter error is unfeasible in full generality. Therefore, we will restrict our attention to simplified cases, though we should emphasize that our approach in practice performs well on significantly broader model assumptions, which we will demonstrate in Section 5. Consider thus the true and the empirical loss functions (24) for the case of the orthogonal design (A A = Id). Without loss of generality, we can assume A = Id (otherwise redefine y to be A y). Let now y = x + σw and assume |y 1 |≥ |y 2 |≥ . . . ≥ |y m |. Plugging (34) into (24) we get 1], and differentiable on intervals Focusing on one interval at a time, a direction computation yields that for k = 1, . . . , m − 1 the minimizer of R| I k is given by i . An analogous statement holds for k = m, whereas R(t) is constant on I 0 , as argued in (32). Therefore, the minimizer of R(t) is t opt = argmin k=0,...,m R(t * ,k ). The empirical loss function R(t) is also continuous on [0, 1] and is piecewise differentiable on the same set of intervals since they depend only on y. Consequently, minimizers t ,k of R(t) are also of the form (35), where we only ought to replace x i by x i .
Notice that unless further assumptions are made, minimizers t opt and t opt are not given explicitly: we still need to evaluate R(t) and R(t) at m + 1 locations. To address and showcase this issue we analyse two cases, that of Bernoulli noise and of bounded sub-Gaussian noise. This analysis will also give a theoretical intuition that will drive our algorithm.

Bernoulli Noise
We first consider a simplified model which however simple still encodes the main features of the problem. In particular, let and assume x = (x 1 , . . . , x h , 0, . . . , 0) , and |x i | ≥ 2σ, for i = 1, . . . , h. Without loss of generality we can assume |y 1 |≥ |y 2 |≥ . . . ≥ |y m |. It then follows b i = 2 |x i ± σ| + 1, for 1 ≤ i ≤ h, and b i = 2σ + 1 otherwise. Moreover, b j > b h+1 , for all j = 1, . . . , h. In the following, we will for the sake of simplicity consider the case when α = 1. The details regarding the general case, α = 1, are in the Appendix. Let us first consider the true loss function. We have The function R(t) is differentiable for all t = b i , i = 1, . . . , h + 1, in which case a direct computation yields We will first show that t opt ≥ b −1 h+1 .
Proof. R is constant on I 0 by definition. For a fixed j = 1, . . . , h and t ∈ I j we have Thus, the claim would follow from We can now explicitly compute the minimizer of R(t).

Lemma 3.2. True loss functional R(t) is minimized for t opt
Proof. For t ∈ b −1 h+1 , 1 we get The root of the expression above is exactly (36). Arguing as in the proof of Lemma 3.1, we have t > 1 b h+1 . Restricting to [0, 1] the claim follows. In case of the empirical loss function, it is in general not true that x i = 0 for i > h, nor is y i − x i σ a Bernoulli random variable. However, y and b i 's remain the same, and an entirely analogous computation gives We can now show a bound on the error in approximation the optimal regularization parameter.
Proof. By assumption t opt = t . Recall that Therefore, using (39) we have with v defined by v j = 2 sgn(y j )b j for j = 1, . . . , h, and v j = sgn(w j ) for j > h. It follows that v 2 ≤ b 2 . Thus, we can write Using b 2 2 = m + 4 y 1 + 4 y 2 2 and Πw 2 = h, Lemma 2.2 provides that with probabilty greater than 1 − 2 exp(−u), N h + u + σ 2 m and C > 0 is a constant. By increasing N, we get that t < 1, so that t = t opt and the claim follows. The proof for a bounded y is bounded, the proof is entirely analogous.
Remark 3.5. The proof of Theorem 3.4 is valid not only for α = 1 but rather for all α Namely, in the appendix A.3 we show that, holds for any α > 0.

Bounded noise
A particularly helpful property of Bernoulli noise is that it reduces the number of distinct intervals I k since b i = b h+1 for all i > h. This made it possible to directly compare t opt and t opt , since all we had to do was to show that they belong to the same interval, and then estimate the error. In case of non-Bernoulli noise this will not be possible since b i = b h+1 does not hold in general. To illustrate this, we now consider the case when the noise is bounded and there is a gap between the noise and the signal. We show that there exists a unique minimizer, which will also justify our numerical algorithm for its estimation. For simplicity of computation, we keep α = 1. As before, let y = x + σw and assume x = (x 1 , . . . , x h , 0, . . . , 0) , and |x i | > 2σ w j , for all i = 1, . . . , h and j = 1, . . . , m. Without loss of generality, we assume that y is ordered so that Consequently, the loss functional R(t) = z t − x 2 2 is piecewise differentiable on intervals I k , k = 0, . . . , m, where now Same as in the case of Bernoulli noise, we will show that R(t) is decreasing 1 for all t ≤ b −1 h+1 . So, let t ∈ I j for j < h. The function R(t) is continuously differentiable in I j , so to show that it is monotonically decreasing it is sufficient to show that R (t) is positive. We have 1 We will show that R(t) can be monotonously increasing only for a large enough t. Thus, for all t smaller than that value (denoted as ϑ j ), it will be a monotonously decreasing function Hence, t opt > b −1 h+1 , as desired. We will now show that R(t) admits only one minimizer. Assume there exists t such that t ∈ I j for some j > h and R (t ) = 0. This means We will now prove by induction that R(t) is increasing on all intervals I j for j > j . For t ∈ I j with j > h, it follows The rest of the proof then follows by mathematical induction. Analogous computation yields the same type of a result for the empirical loss function R. Lemma 3.6. Let the above assumptions hold. Loss function R(t) is then either monotonically decreasing on the entire interval [0, 1], or it is decreasing until some interval I j , for j > h + 1 where it achieves a (unique) minimum, and it is monotonically increasing on all the subsequent intervals. The same holds for R. Lemma 3.6 says that R(t) and R(t) achieve a unique minimum in [0, 1], and that they are monotonically decreasing before, and monotonically increasing after this minimum. Furthermore, the minimizer is bigger than (1 + 2σ |w h+1 |) −1 , which means that for moderate noise levels, it will be close to 1. The issue though is that minimizers ϑ j and ϑ j do not need to be in the same interval, that is j = j . This complicates having a direct comparison, which we will thus avoid.

OptEN Algorithm
Driven by insights in Section 3, we are ready to present an efficient heuristic algorithm for learning the Optimal regularization parameter for the Elastic Net (OptEN). The algorithm is based on the minimization of a given loss function (be it e.g. the empirical loss function (24) or the modified projected risk (28)). In Section 3 we showed that in a simplified, yet instructive, setting that the optimal parameter tends to be in the vicinity of t = 1 (depending on the noise level and the signal-to-noise gap). This is supported by experimental evidence in more general situations such as non-injective A, as we will see in Section 5. Moreover, the loss function is monotonically decreasing as we get away from t = 1, and experimental evidence suggests that the angle of the slope near t = 1 is steep and depends on σ. These observations drive our algorithm which assumes that the minimizer lies in a steep-sided valley not too far from t = 1, see Figure 1. Therefore, we will perform a line search on the graph of a given loss function, starting from t = 1.
Line search methods follow iterations t k+1 = t k + s k p k , where p k is the search direction and s k the step size: • Search direction. We estimate R (t), for t = t k , by central differences, where > 0. For t = 1 we compute the one-sided differences∆ R(1) := R(1)−R(1− ) . Then p k is given as p k = −∆ R(t).
• Step size. We estimate s k by using the Armijo rule with quadratic models and the backtracking line search (consult [2] for an overview of line search methods).
Our approach is presented in Algorithm 1, while the extensive numerical study on synthetic and imaging data is provided in the next section.

Algorithm 1 OptEN algorithm approximated the optimal elastic net regularization parameter using backtracking line search
Input: y 1 , . . . , y N , y ∈ R m , A ∈ R m×d ; Output: Approximate optimal regularization parameter t opt ; Compute x according to (18). Set a loss function → R or R P or R M . In the rest of the algorithm we will refer to it as R; until |p k | < tol or k < max_iter; Output: Approximate regularization parametert := t k .

Experimental Results
To highlight the performance of our approach and its adaptivity to different scenarios we conduct experiments on synthetic and real imaging data. In the first set of experiments, we perform a thorough comparison of our method with the state-of-the-art parameter selection rules by exploring their sensitivities with respect to noise level and other notions. The second set of experiments deals with image denoising where we use wavelet thresholding with elastic nets. We consider two data-sets: natural images and a real-world brain MRI data. What is important in the latter is that we cannot assume homoscedasticity of given data samples, but our method still delivers near-optimal performance. Note that we do not aim to compare our method with state-of-the art denoising methods, but rather only with state-of-the-art methods regarding the selection of the regularization parameter for the elastic net. Before presenting the experimental results, we start with a discussion of methods that can be used for automatic detection of the sparsity level h and show that when a sufficient amount of training points is given, we can reliably estimate h from the data.

Estimating the Sparsity Level
In real applications the sparsity level of a vector is either not available or is only an approximate notion, i.e. the desired vector is sparse only when we threshold its entries. Such situations require h to be estimated from available data, which in our case means looking at the spectrum of the corresponding covariance matrix. This question belongs to the class of low-rank matrix recovery problems since what we are trying to recover is the geometry (i.e. projection onto the range) of the noiseless Σ(Ax), which is of lower rank, by using only the the covariance matrix of noisy observations Σ(y), which is of full rank. Thus, the problem of estimating h boils down to thresholding singular values of the empirical covariance matrix according to some spectral criteria that exploits the underlying structure. ; (c) the cumulative spectral In the latter two cases one sets a value, say 0.95, and selectsh as the first k for which the corresponding spectral energy reaches that threshold.
We study the behaviour of these four metrics on three different types of forward matrices A: random Gaussian, random circulant Rademacher, and random Toeplitz Gaussian matrices. These matrices were chosen because they have different spectral behaviour and commonly appear in inverse problems. In each case A is a 100 × 100 real matrix, normalised so that A = 1, and we take N = 100 samples x i , sampled according to (D2) for h = 20, w i ∼ N (0, I 100 ) and σ = 0.3. We compute Σ(y) = 1 N ∑ N i=1 y i ⊗ y i for y i = Ax i + σw i . In Figure 2, we show the application of the aforementioned spectral criteria to the singular values of Σ(y). From the results it is clear (as it is generally understood) that all four methods are flawed and would fail if used without taking further information into account. For example, the spectral gap criterion (in the upper right panel) would in the strictest sense dictate the selection of h = 1, but a more careful look at the plot suggests that the behaviour of the spectral gap changes dramatically around h = 20, which would be a better choice. Such ad hoc solutions are sensible and often lead to great gains in performance, but can be hard to quantify, especially on real data.
The last spectral criteria, 1 − λ k λ k+1 , is perhaps the most promising, but is also subject to demands on N, as shown in the bottom row of Figure 2. Namely, if N is not large enough, then computing 1 − λ k λ k+1 on the spectrum of Σ(y) has a heavy tail, and thus, if the criterion is applied blindly we could severely overestimate h, as in the bottom left corner of the Figure. An ad hoc solution in this case would be to look for the relative gap within the first m/2 singular vectors, which is the approach used in the results in Section 5.2 to estimate h. Moreover, in the case of the first three spectral criteria the situation does not change as the number of samples increases. On the other hand, for the last criterion it does: heavy tails flatten back to zero for all three choices of random matrices, see the plot in the bottom right corner of Figure 2.
Recovery of low rank matrices by thresholding the singular values has also been studied by Gavish and Donoho in [20]. Let λ med be the median of singular values λ i , for i = 1, . . . , m, of Σ(y), and defineh as h = #{λ i : λ i ≥ 2.858λ med }, i.e. as the number of singular values greater than 2.858λ med . The authors show that this rule produces an asymptotically better mean squared error for the recovery of a low rank matrix than what we would get with using the true h. For the present testing scenario, for N = 100 this rule estimates h as 23, 24, and 23 for the three random matrices, and for N ≥ 150 this criterion always selects h = 20 for each of the matrices.
In Section 5.3 we will apply our algorithm to wavelet denoising, where the decay of corresponding eigenvalues is smooth and it is exceedingly difficult to decide where to make the eigenvalue cut using the rules mentioned above. Namely, in such cases there is no natural choice of h that could be inferred by looking at the singular values and none of the methods mentioned above gives a consistently effective rule for choosing h. We will therefore consider two testing scenarios; when we are given an oracle h (by which we mean the h giving the highest PSNR), and when h has to be estimated from the data.

Synthetic Examples
Experimental Setting. We consider the inverse problem of the type y = Ax + σw, where the data are generated according to Comparison. We compare our algorithm with the state-of-the-art parameter selection methods such as the discrepancy principle [28], monotone error rule [31], quasi optimality [32], L curve method [23], (Monte-Carlo) balancing principle [26] and its elastic net counterpart [10], (Monte-Carlo) generalised cross-validation [21] and nonlinear cross validation [19]. In the remainder of this paper we will refer to other methods by their acronyms, and to our method as OptEN -a method that uses Algorithm 1 on the empirical loss function R. The first five methods are commonly used in inverse problems (a detailed account and an experimental study can be found in [3]), whereas Monte-Carlo and nonlinear cross-validation are adaptations of generalised cross-validation for non-linear regularization methods.
Before presenting the comparison results, we provide a concise description and the rationale that underpins the considered methods. Most of these methods require some additional information about the problem, predominantly the noise level σ, to be either known or estimated. The performance of these methods is affected by the accuracy of that information. We will provide the true noise level whenever a given method requires it and furthermore, we will perform judicious testing and tuning of all other quantities, taking into account recommendations from relevant literature.

Discrepancy Principle [DP]
Discrepancy principle is one of the oldest parameter choice rules which selects a solution so that the norm of the residual is at the noise level. Thus, the regularization parameter is chosen by the first n ∈ N such that where we fix τ = 1.

Monotone Error Rule [ME]
This rule is based on the observation that the monotone decrease of the error z n − x 2 in the noisy case can only be guaranteed for large values of the regularization parameter. Therefore, the best parameter t n * is chosen as the first t-value for which one can prove that the error is monotonically decreasing. This means that we can choose the parameter by the smallest n such that We fix τ = 1 for our experiments. The left hand side of (43) is replaced with (42) whenever the denominator is 0.

Quasi-Optimality Criterion [QO]
Quasi-optimality is a parameter rule that does not need the noise level, and thus has enjoyed a reasonable success in practice, esp. for Tikhonov regularization and truncated singular value decomposition. The regularization parameter is chosen as

L-curve method [LC]
The criterion is based on the fact that the log − log plot of ( Az n − y 2 , z n 2 ) often has a distinct L-shape. As the points on the vertical part correspond to under-smoothed solutions and those on the horizontal part correspond to over-smoothed solutions, the optimal parameter is chosen at the elbow of that L-curve. There exist several versions of the method; here we use the following version n * = argmin

(Monte-Carlo) Balancing Principle [BP]
The principle aims to balance two error contributions, approximation and sampling errors, which have an opposite behaviour with respect to the tuning parameter. More precisely, we select the parameter by n * = argmin{t n | z n − z k 2 ≤ 4κσρ(k), k = n, . . . , N max }, here κ > 0 is a tuning parameter. Alternative versions [3] of the balancing principle were proposed that achieve the same accuracy but are more computationally efficient. As our main focus on the accuracy of the parameter choice, we will use the original and more computationally heavy version of the balancing principle.
The value of σρ(k) is, in general, unknown but can be estimated from data in case of white noise. Following [3], we calculate , and A −1 n is the map that assigns y to z n .

Elastic Nets Balancing Principle [ENBP]
In [10] the authors proposed the balancing principle for elastic net. In particular, The method stops the first time two solutions are sufficiently far apart. The constant C needs to be selected, and in our experience this task requires a delicate touch.

(Monte-Carlo) Generalised Cross-Validation [GCV]
The rule stems from the ordinary cross-validation, which considers all the leave-one-out regularized solutions and chooses the parameter that minimizes the average of the squared prediction errors. The appealing fact of GCV is that the calculations can be done without computing all the regularized solutions. Specifically, GCV selects n according to where A −1 n is the map such that z n = A −1 n (y). In the case of elastic nets the map A −1 n is not linear and, thus, we cannot assign a meaning to its trace. Instead, we follow the ideas of [21] and estimate the trace stochastically using only one data sample.

Nonlinear Generalised Cross-Validation [NGCV]
Another way to reconfigure GCV for non-linear methods, namely shrinkage models, was described in [19]. The NGCV selects n as  Table 1: Comparison of errors for regularization parameter selection methods, with an injective matrix A ∈ R 500×100 and h = 10, α = 10 −3 , σ = 0.3. The values are averages over 100 independent runs Comparison/Error. For each method we computed: the (normalised) error in approximating the optimal regularization parameter, the (normalised) error, false discovery proportion (FDP), true positive proportion (TPP), and the computational time. TPP and FDP are measures that quantify the notions of true and false discovery of relevant features in sparsity based regression tasks [30]. FDP is the ratio between false discoveries and the total number of discoveries, TPP on the other hand is the ratio between the number of true (i.e. correct) discoveries and the number of true discoveries in the original signal, The closer FDP is to 0 and TPP is to 1, the better is the parameter choice in recovering the structure of the original sparse data. It is known that there is often an explicit (and sometimes even quantifiable) trade-off between FDP and TPP, in the sense that the support overestimation is an (undesirable) side-effect of the full recovery the support of the original (sparse) signal. In other words, a consequence of true support discovery is often a non-trivial false support discovery (see [30] for the case of the Lasso path). When computing FDP and TPP we will rather than demand for an entry to be exactly zero, instead threshold the values (with 0.5 being the threshold).
The results in Table 1 are averaged over 100 independent runs for α = 10 −3 , m = 500, d = 100, h = 10 (where V = span {e 1 , . . . , e h }), and σ = 0.3, which corresponds to SparseSNR ≈ 0.17. The first row in the table, x t opt , describes the behaviour of the elastic net minimizer for which the true optimal regularization parameter is provided.
Discussion. OptEN always returns the value which is the closest to the optimal regularization parameter, and its results are in general comparable to the ones provided by the minimizer with the optimal parameter. However, one can observe that other methods, e.g., discrepancy principle, provide a better balance between FDP and TPP (returning solutions that are more sparse), though at a cost of a larger approximation error. Balancing principle also provides very good results, but it is slow unless an effort is made to improve its computational time. Moreover, we have observed that the performance of all methods that require the noise level σ to be known deteriorates if we do not provide the exact value of the noise level, but only its rough estimate.
The overall results are mostly consistent over all experimental scenarios we looked at, with a couple of exceptions. As expected, FDP and estimation errors deteriorate not only for larger σ but also for larger α, though the ranking of the methods and the patterns of behaviour remain the same. This is because as α increases elastic nets sacrifice sparsity for smoothness. The empirical estimator x is a very accurate estimator of the original signal, and it sometimes outperforms both our method and indeed the elastic net solution that uses the optimal parameter. However, as has been observed in [11] for Tikhonov regularization and confirmed in Figure 3, the performance of the empirical estimator worsens in the small noise regime.

Comparison with empirical estimator: effects of σ and N train
We study the behaviour of the relative estimation error with respect to σ and the number of training samples. We compare OptEN with the empirical estimator x, DP, NGCV, and BP. We use m = 500, d = 100, h = 10 and (D1)-(D3) with σ ranging from 0.1 to 0.5 in the first experiment, whereas we vary the number of training samples N from 20 to 60 in the second experiments as depicted in Figure 3. Our method again outperforms other considered parameter selection rules. On the other hand, the empirical estimator has a slightly better performance than OptEN for larger noise levels (it is also better than the elastic nets solution with the optimal parameter), and it is worse for lower noise levels. This is essentially due to the fact that x is never truly sparse, but has a lower (thresholded) FDP. Namely, as a projection of a rank h the non-zero entries of x are very small, whereas the non-zero entries obtained with elastic net are larger, specifically with their size depending on the noise level. In other words, OptEN and elastic net, in general, are more affected by the increasing noise level. The parameter α plays a similar role; for small α OptEN beats x, and for larger α the situation is reversed.

Non-injective matrices
We also conduct experiments with non-injective matrices. The setting is as in Table 1, with the only difference that now A ∈ R 500×100 with rank (A) = 40. As mentioned in Section 2.2.1, we test our method by minimizing the projected loss functional R P and the modified error functional R M .
The results can be found in Table 2. Our method (using both the projected and modified functionals) again outperforms standard parameter selection rules in terms of the precision accuracy, and loses out to some methods when it comes to FDP and TPP. One can also observe that the performance of the empirical estimator deteriorates and x indeed should not be used as the solution itself but some additional regularization is required.    Table 2: Comparison of errors for regularization parameter selection methods, with a matrix A ∈ R 500×100 , rank(A) = 40, and h = 10, α = 10 −3 , σ = 0.3. The values are averages over 100 independent runs

Image Denoising
In this subsection, we demonstrate the practical value of the proposed algorithm by showing that it accurately and efficiently selects the regularisation parameter in case of image denoising with elastic nets. Image denoising is a classical problem in signal and image processing, which aims at finding an estimate Z of an unknown image X from its noisy measurements Y, where Y = X + σΞ and Ξ denotes (isotropic white) noise with covariance Id.
The goal is to improve the image quality by removing noise with preserving important image features such as edges and homogeneous regions. A vast amount of approaches have been proposed to deal with the image denoising problem, starting from simple and 'classical' wavelet thresholding [13,14] and non-linear filters to stochastic and variational methods [8,9]. Since the primary goal of this paper is to evaluate the performance of the proposed approach as an accurate and computationally efficient method for parameter selection in elastic net regularisation, we perform a thorough comparison of our method with other state-of-the-art parameter selection methods for elastic nets. In particular, we compare to the elastic net denoising schemes where parameters are selected using discrepancy principle, balancing principle, non-linear cross validation as well as an oracle (an optimal parameter). In all cases, we show that OptEN selects a parameter that leads to a better denoising results compared to other methods and even denoising with an oracle parameter, see Table 3.
Wavelet-based denoising. We will denoise the observed noisy image Y by minimising the elastic net func- where W is the wavelet transform corresponding to the family of db4 wavelets, and α = 10 −3 .

Denoising with an Oracle h
In case of wavelet denoising and other real-world scenarios, there is often no natural notion of h. By natural we here refer to any notion that would rely solely on the spectrum of the covariance matrix of noisy samples to represent its low rank structure or the sparsity of the original signal. Instead, we say that the true h is the one which minimises the MSE (i.e. maximises the PSNR) of the reconstructed image. In our first set of experiments, we assume that the optimal h is known, which means that we first select the best possible h and then use it to denoise the image.
We should emphasise the fact that by this procedure each of the training images has a random noise level, which is different from the noise level of the test image.
To learn the regularisation parameter, we first rescale all images, i.e. the test image and the training images, to 32 × 32 pixels and then compute the regularisation parameter in (47). This rescaling is a trade-off that sacrifices accuracy for the sake of faster computations. We should also note that rescaled images have a different noise level, which is 1.5 to 2 times smaller than the noise level of the original image. For BP and DP the regularisation parameter is selected from a sequence of parameter values t n = 1 1+µ 0 q n with µ 0 = 1, q = 0.95, and N max = 100, same as before. Moreover, we fix τ = 1, κ = 1/4 for BP and DP, and provide them with the true noise level. Our algorithm is learned using N = 1800 training images 3 for each given test image, using the method described above, and the empirical estimator is computed with an oracle h, i.e. the one returning the lowest MSE. Once the regularisation parameter is selected, we use its value to denoise the full-sized image.
Comparison/Error. We measure the recovery success of the methods by two performance metrics: peak signal-to-noise ratio (PSNR) and the similarity index (SSIM) between the original image X and the recovered version Z. PSNR is a standard pixel-based performance metric, defined through the MSE by PSNR(X, Z) = 10 log 10 max 1≤i≤p X i − min 1≤i≤p X i MSE(X, Z) , MSE(X, Z) = 1 p X − Z 2 .
Both MSE and PSNR are ubiquitous in image and signal analysis due to their simplicity, suitability for optimisation tasks, and other properties, but are also infamous for their inability to capture features salient for human perception of image quality and fidelity [37]. SSIM on the other hand, is a structure-based performance metric that tries to address this issue by using easy-to-compute structural statistics to estimate image similarity. It is defined through where X, Z are the means, and std(X), std(Z) are the standard deviations of pixels of corresponding images X and Z, and C 1 , C 2 are small positive constants 4 . 3 A significantly smaller number of samples produces comparable results. In this study we chose to use 1800 samples since the wavelet decomposition of a 32 × 32 image has 2025 coefficients 4 We take C 1 = 0.01, C 2 = 0.03 by the convention of python's skimage package  Table 3: Results on wavelet denoising of noisy images with different noise levels using elastic nets minimisation. Each column defines the method used to select the regularisation parameter. In bold is the method that achieved the best result. Columns titled noisy correspond to PSNR and SSIM values for the initial noisy image. Columns titled using λ opt correspond to the best values achievable for the selected elastic nets functional (where we find the optimal parameter by a grid search on the true loss functional Table 3 provides the PSNR and SSIM values generated by all algorithms on the considered images, while Figure 4 demonstrates a denoising result on a 128 × 128 detail of each image for σ = 0.075. We can see that our method achieves the highest PSNR on all images and that this effect is more pronounced for larger noise values.

Denoising with a heuristically chosen h
In this set of experiments, we challenge the performance of our method in a scenario where h cannot be chosen optimally and the training samples do not correspond to the same image, but only to the images with a (presumed) similar structure. Furthermore, the number of training samples is limited (significantly smaller than in the Section 5.3.1), and once again, the training images all have different noise levels.
In order to test this scenario, we run our experiments on a real-world dataset of brain images 5 . This data set consists of in vivo MRI images of 13 patients with brain tumour, taken pre-surgery. For each patient we took 25 consecutive MRI slices (resulting in a grand total of 325 images), and then isolated the area around the brain. All images where then resized to a uniform size of 128 × 128 pixels. From the resulting data set, we then create a test image by randomly selecting one of the images and adding to it isotropic white noise with σ ∈ {0.05, 0.075, 0.1}. The training set is constructed as in the previous section, by adding additive isotropic white noise to each of the remaining images, where in each instance the noise level is drawn uniformly at random from [σ − 0.05, σ + 0.05]. The sparsity level is then estimated by computing the PSNR of the regularised solution for h = 20 and subsequent h + 10 values. We stop when either PSNR starts to decrease, or h becomes larger than 324 (i.e. the larger than the size of the data set). The results for σ = 0.1 can be seen in Figure 5. The effects of denoising are visually not as striking as the results in Section 5.3.1. We attribute this to the fact that PSNR gains with the best possible choice of parameter using elastic net are quite small, namely, PSNR of the noisy image improves only by around 5 − 7%, when taking the optimal parameter (see Table 4). Other parameter selection rules (discrepancy principle and balancing principle) did not improve the PSNR and, thus, are not presented in the study. noiseless noisy OptEN Figure 5: Denoising results on a brain image data set. The PSNR values are in Table 4 6 Conclusion and Future work In this paper, we presented an approach for the estimation of the optimal regularisation parameter for elastic net. The theoretical guarantees are possible only in simplified scenarios but we used insights gained therein to steer and create an efficient algorithm. This algorithm exhibits good prediction accuracy, including cases when there  Table 4: PSNR values for the results in Figure 5 are no theoretical guarantees. Comparison with state-of-the-art methods show a clear superiority of our method, under the studied testing scenarios. Moreover, whereas other studied methods have a number of additional, intrinsic, parameters that ought to be adjusted to achieve satisfactory results, our method is autonomous given a sufficient number of training samples. We aim to use these as guiding principles for other situations. Namely, we will study the behaviour of the solution with respect to the other hyperparameter; α. Also, we consider and other optimisation schemes, predominantly focusing on imaging applications. Our initial work on TV minimisation is promising. We also aim to consider other geometries, e.g. when V is a union of subspaces, or a lower dimensional manifold.

A.2 Proofs for Theorem 2.2
We will here add an analogue of equation (12) for the case of bounded y. Assume y 2 ≤ √ L holds almost surely and consider a random matrix R = y y. Then ER = Σ(y), and R ≤ L. Furthermore, R = R and Let now y i ∼ y and define a family of independent m × m matrices R i = y i y i , i = 1, . . . , N, so that R i ∼ R. The empirical covariance Σ(y) = 1 N ∑ N i=1 R i is then the matrix sampling estimator and by Corollary 6.2.1 from [34] we have that for all s ≥ 0 P Σ(y) − Σ(y) ≥ s ≤ 2m exp − Ns 2 /2 m 2 (R) + 2Ls/3 .

A.3 Computations for α = 1
We now go back to the discussion in Section 3.1. So, let y = x + σw, where P (w i = ±1) = 1 2 , and assume x = (x 1 , . . . , x h , 0, . . . , 0) , and |x i | ≥ 2σ, for i = 1, . . . , h. In the following we will use v 1:k to denote a vector in R k that consists of the first k entries of a vector v ∈ R m , and denote η = σw. In Section 3 we showed that the minimum of R(t) = z t (y) − x Writing down each of the terms in the numerator we get d i c i − d i c i = −2αerr i (y i − x i ) − err i d i = −err i a i , where we use the fact that c i = d i + 2(y i − x i ). We also get d i c j + d j c i − d i c j − d j c i = err j ((α − 1)d i − αc i ) + err i (α − 1)d j − αc j = − err j a i + err i a j . Thus, a j err j , due to c i c j + c i c j = 2c i c j + (α − 1) err j c i + err i c j , and c i c i = c 2 i + (α − 1)c i err i . Putting it all together and rewriting we have