Detangling robustness in high dimensions: composite versus model-averaged estimation

Robust methods, though ubiquitous in practice, are yet to be fully understood in the context of regularized estimation and high dimensions. Even simple questions become challenging very quickly. For example, classical statistical theory identifies equivalence between model-averaged and composite quantile estimation. However, little to nothing is known about such equivalence between methods that encourage sparsity. This paper provides a toolbox to further study robustness in these settings and focuses on prediction. In particular, we study optimally weighted model-averaged as well as composite $l_1$-regularized estimation. Optimal weights are determined by minimizing the asymptotic mean squared error. This approach incorporates the effects of regularization, without the assumption of perfect selection, as is often used in practice. Such weights are then optimal for prediction quality. Through an extensive simulation study, we show that no single method systematically outperforms others. We find, however, that model-averaged and composite quantile estimators often outperform least-squares methods, even in the case of Gaussian model noise. Real data application witnesses the method's practical use through the reconstruction of compressed audio signals.


Introduction
We investigate the benefits of model-averaged as well as composite estimators in highdimensional problems where the underlying goal is superior prediction quality. Robustness in data analysis with potentially more parameters than samples is a critical practical question and is of particular interest in constructing recoveries of compressed images and signals which should have high precision.
Model averaging, often used as a first tool to improve estimation quality, forms a weighted average of estimators and is here utilized for regularized sparsity-encouraging estimation in a high-dimensional regression setting. Model averaging is also well-known in the Bayesian setting (Hoeting et al., 1999), though we focus on its frequentist version in which a user determines the weights assigned to the separate estimators (e.g., see Claeskens and Hjort, 2008;Hjort and Claeskens, 2003;Yuan and Yang, 2005). Model averaging enjoys a wide application, see, for example, the recent overview paper for model averaging in ecology by Dormann et al. (2018) and for application to hydrology by Hge et al. (2019). In econometrics, the terminology "forecast combinations" appears (e.g., in Cheng et al., 2015;Bates and Granger, 1969); whereas "multimodel inference" is another commonly used term for this procedure (Burnham and Anderson, 2002).
While the technique is quite thoroughly investigated for low-dimensional models, far fewer results have been obtained in high dimensions. Ando and Li (2014) consider highdimensional linear regression. By computing the marginal correlation between each covariate and the response and forming groups according to the obtained values, regularized estimation is avoided. The authors fit a fixed number of low-dimensional models by the least squares method and subsequently average them. Zhao et al. (2016) extend this method to dependent data, while Ando and Li (2017) extend this approach to generalized linear models, again by only fitting low-dimensional models, this time via maximum likelihood estimation. In these papers, the weights are obtained via cross-validation; see also Hansen (2007) and Hansen and Racine (2012) for similar weight finding approaches in low-dimensional models.
Our setting is different and is theoretically valid (see Theorem 1 below). We explicitly work with l 1 -regularized estimators that are averaged, and we do not rely on the correct low-dimensional representation of the model. When designing the optimal weights, we explicitly take variable selection effects (of regularization itself) into account. Is the dependence among regularized estimators an impediment or a hidden benefit in obtaining robust predictions, i.e., predictions that do not change much when the data is changed a little?
A second approach to robustness is through composite estimation. While model averaging combines estimators after optimization of their respective loss functions, composite estimation weights the loss functions directly (before optimization). For quantile regression in low dimensions, Koenker (2005, Theorem 5.2) stated the asymptotic equivalence of model-averaged and composite quantile regression estimators, provided each method uses its own, optimal set of weights that minimize the asymptotic variance. Hence, with optimal weights, there is no asymptotic preference between the two methods in low dimensions. For high-dimensional quantile regression, when one restricts the attention to inference regarding the true nonzero part of the regression coefficient and ignores the variable selection effect, Bloznelis et al. (2019) obtained the same equivalence for high-dimensional quantile regression using different types of regularizations (SCAD, lasso, adaptive lasso).
In practice, however, one works with an estimated coefficient vector for which one is not sure that the regularization has led to the correct selection. Therefore, incorporating imperfections of variable selection is especially important for achieving robustness. This is where our approach differs from Bloznelis et al. (2019) or Bradic et al. (2011), where an irrepresentable condition (needed for consistent model or asymptotically perfect selection) has been used to specify weights and analyze robustness. The approximate message passing (AMP) algorithm is crucial in our approach to take the variable selection into account when studying the estimators' asymptotic mean squared errors. The use of such algorithms has been investigated by Donoho et al. (2009) and Bayati and Montanari (2011a) for compressed sensing. Donoho and Montanari (2016) explain the use of AMP algorithms for obtaining the variance of high-dimensional M-estimators for which n/p → δ ∈ (1, ∞). Here, n denotes the sample size and p the number of regression coefficients. However, the robustness of sparsity encouraging AMP estimators is still largely unknown.
In this paper, we first extend the robust AMP (RAMP) of Bradic (2016) to regularized composite estimation. Second, we construct estimators and develop new theory for the asymptotic mean squared error (AMSE) both for model-averaged and for composite estimators. Note that model-averaged AMSE required an extension of AMP theory for a challenging case of dependent estimates. Besides, we establish new Stein-type risk estimates of the AMSE in both cases.
The new estimates of the AMSE of the model-averaged and composite estimators enable a theoretically justified and data-driven optimal weight choice by minimizing the estimated AMSE (without relying on perfect variable selection). The estimated AMSE provides more information regarding the estimators than merely considering which variables have been selected.
Organization of the paper. First, in Section 2, we detail the model-averaged and composite estimators in a high-dimensional setup. Next, we explain the model-averaged robust message passing algorithm in Section 3. The limiting behavior of the estimators in the algorithm is studied by state evolution parameters in Section 4. We obtain the estimators' asymptotic mean squared error as well as an estimator of that quantity in Section 5. We showcase the procedure for high-dimensional regularized quantile regression in Section 6 and present numerical results in Section 7. Section 8 concludes. All proofs, together with the assumptions and some technical lemmas, are collected in the Appendix.

Model-averaged and composite estimation
We consider a high-dimensional linear model Y = Xβ + ε with Y ∈ R n , the design matrix X ∈ R n×p and the parameter vector β ∈ R p . The ith row of X is denoted X i· , i = 1, . . . , n, the jth column of X is denoted by X ·j , j = 1, . . . , p. We assume the components of ε to be independent and identically distributed with mean zero, cumulative distribution function F ε and probability density function f ε . We allow for a sparse high-dimensional setup.
Denote by s the l 0 norm of the parameter vector, s = β 0 , which counts the number of nonzero components of the vector β. We assume that the ratios n/p → δ ∈ (0, 1) and n/s → a ∈ (1, ∞) when p, n, s tend to ∞.
We consider two types of weighted estimation methods. First, model-averaged estimation where estimators from different models or estimation methods are weighted and summed to arrive at a final estimator, see (2). Second, composite estimation where a weighted average of loss functions is minimized; see (3).
For model-averaged estimation of the parameter β, define for k = 1, . . . , K the regularized estimators β k (λ k ) = arg min where ρ 1 , . . . , ρ K are nonnegative convex loss functions and λ = (λ 1 , . . . , λ K ) is a vector of possibly different nonnegative regularization parameters. For a set of weights w = (w 1 , . . . , w K ) , the model-averaged estimator is defined as Often one assumes that the weights w 1 , . . . , w K are all nonnegative and sum to 1, although this is not necessary for the computation of the estimator. For composite estimation we consider again K loss functions, though only with a single nonnegative regularization parameter λ, such that the regularized composite estimator is defined as β C (λ) = arg min Computationally, composite estimation is harder than model-averaged estimation and requires that all weights are positive to ensure a nonnegative and convex weighted loss function, even when all ρ k are nonnegative and convex. Hence, for composite estimation it is required that the weight vector w ∈ [0, 1] K such that K k=1 w k = 1. As a worked-out scenario throughout the paper, we consider quantile loss functions ρ k (·), k = 1, . . . , K that are defined below. For more information about quantile regression with i.i.d. errors, see Koenker (2005, Sec. 3.2.2). In this paper, we assume that the design matrix X does not contain a column of ones; see assumption (A1) in the Appendix. With τ ∈ (0, 1), the τ -quantile of the response Y is obtained as Xβ + F −1 ε (τ ) = Xβ + u τ . Figure 1 presents first a single quantile loss function with τ = 0.3,

Equal weights
Unequal weights For model averaging we specify K different quantile levels and use K different such quantile loss functions for estimation of β: For composite quantile estimation we assume that the K quantile levels τ 1 < · · · < τ K , then also the quantiles of ε are sorted u τ 1 < · · · < u τ K . Define u τ 0 = −∞ and u τ K+1 = ∞. The middle panel of Figure 1 depicts such a composite quantile loss function ρ C = K k=1 w k ρ k for K = 3 quantile levels 0.25, 0.5 and 0.75 with equal weights w = (1/3, 1/3, 1/3) . The panel on the right in Figure 1 uses the same quantile levels but depicts the quantile loss function ρ C with weights w = (0.15, 0.55, 0.3) .
In general, the composite quantile loss function can be rewritten in the following way, Note that a single quantile loss function can be seen as a particular case of a composite loss function: take K = 1 and the single weight w 1 = 1. Theoretical results regarding regularized estimation for a single quantile loss function can be found in Bradic (2016). Henceforth, we concentrate on the example of the composite case.
One aim of this paper it to investigate the weight choice w by minimizing the asymptotic mean squared error of the estimators β MA (λ) and β C (λ).

Robust approximate message passing
The idea behind approximate message passing algorithms is to provide an iterative procedure that has as its fixed point the estimator of interest; in this case the minimizer (2) of the regularized loss function in the case of model averaging, and the estimator (3) in the case of composite estimation. Due to a convergence in the mean square between the solution of the approximate message passing algorithm and the estimator (2), respectively (3), the asymptotic mean squared error that holds for the solution of the approximate message-passing algorithm, is also the asymptotic MSE of the other estimator. Studying effects of regularization while allowing n/p → δ ∈ (0, 1) is challenging. The AMP provides theoretical advantages in these cases as it enables a complete and tractable, albeit challenging, structure for obtaining AMSE. This paper is the first to obtain and use the asymptotic mean square error of the regularized estimators to optimize the weight choice of both the model-averaged estimator and the composite estimator. We extend the theory of the RAMP to apply to the model-averaged estimator; see Theorem 1. Challenges arise with incorporating dependence into the AMSE expression; see Theorem 2. Theorem 2, in turn, leads to a new Stein-type estimator of RAMPs asymptotic MSE. While we focus on the weight choice, the availability of an estimated AMSE may be used in other contexts, for instance, for the construction of confidence intervals.

Notation
When the composite loss function ρ C = K k=1 w k ρ k is used in the RAMP algorithm with tuning parameter α we denote the estimator at iteration number t by β C,(t) (α). When the value of the tuning parameter is clear from the context, we also denote the RAMP estimator by β C, (t) .
For constructing the model averaging estimator we denote the separate estimators from the RAMP algorithm using regularity parameters α k , k = 1, . . . , K by β k,(t) (α k ) and the model-averaged estimator is denoted by β MA,(t) (α) = K k=1 w k β k,(t) (α k ) with α = (α 1 , . . . , α K ) . When the value of the tuning parameters is clear from the context, we denote the model averaging RAMP estimator by β MA,(t) .
A generic estimator, without referring to a specific loss function or construction, is denoted by β (t) , using tuning parameter α; the subscript (t) refers to the iteration number.

The robust approximate message passing algorithm
We first revise the (robust) approximate message passing algorithm, which consists of three steps iterated until convergence. In comparison with the more straightforward AMP for the case with a differentiable convex loss function (Donoho et al., 2009), this procedure for robust high-dimensional parameter estimation (Donoho and Montanari, 2016;Bradic, 2016) adjusts the residuals to incorporate the valid score directly. While more details are given in Algorithm 1, which is applied to the different loss functions ρ 1 , . . . , ρ K and to their weighted sum ρ C = K k=1 w k ρ k , we here provide the main outline. The used notation does not explicitly indicate a dependence on the number of coefficients p to not overcomplicate the formulas. Donoho and Montanari (2016) proposed to use the following proximal mapping operator to adjust the residuals. With b > 0, which minimizes the square loss regularized by the non-differentiable loss, ρ. The parameter b controls how the proximal operator map points to the minimum of the non-differentiable loss, where small values correspond to a small movement towards the minimum of ρ.
The fixed point solution of the proximal operator coincides with the minimum of the loss function ρ. For more information, see Parikh and Boyd (2014). We continue with the worked out example on quantile regression, see (4). For = 0, . . . , K, define where we define a summation sign to be equal to zero in the case where the upper summation index is smaller than the lower one, that is, b i=a x i = 0 if b < a. The proximal operator for the composite quantile case, see (4), is See Section B.2.1 for the derivation of the algorithm. We now describe the three steps in more detail.
Step 1: Create adjusted residuals. We use the estimates β (t−1) and β (t) from iteration steps t−1 and t to compute the adjusted residuals where the soft-thresholding function η(x; θ) = sign(x) max(|x|−θ, 0) and the score function G is defined in (11). In Algorithm 1, see Section 4, we give details on how to set the soft-thresholding parameter θ, which might change in each iteration, and we explain that a proper choice of θ as a function of the regularity constant λ leads to an equivalence of the RAMP estimator and the regularized estimator.
The effective score function used in Donoho and Montanari (2016) is Throughout, we use ∂ 1 as the notation for the partial derivative or partial subgradient of a function with respect to its first argument. Functions (e.g. G) are applied componentwise to vectors. For the example on composite quantile regression the subgradient of ρ C is computed as, where h( ) is defined in (5). The effective score function for composite quantile regression, see Section B.2.2, is To incorporate the sparsity, Bradic (2016), see also Bayati and Montanari (2011a), used the rescaled, min regularized effective score function, where ω = E[ B 0 0 ], see condition (A2) in the Appendix, which corresponds to the limit of s/p, with s = B 0 0 , the true number of nonzero components, as p tends to infinity.
Step 2: Use the effective score function to set b. We choose the scalar b (t) such that the empirical average of the effective score function G(z; b) has slope 1, thus n −1 n i=1 ∂ 1 G(z i,(t) ; b (t) ) = 1. In the case of a non-differentiable loss function, Bradic (2016) proposed to solve ν(b (t) ) = 1 with Condition (A3) in the Appendix defines γ l , r l and the differentiable functions v 1 and v 2 (See also Condition (R) of Bradic, 2016), f P is the estimated density of Prox(z i,(t) ; b (t) ) for i = 1, . . . , n.
The derivation of the estimator ν(b (t) ), see also Section B.2.4, relies on the limiting behaviour of the system, see Section 4.
For the composite quantile loss, see (4), we clearly see the dependence on the quantiles. The estimator of ν in (12) uses v 1 (z) = 0 and v 2 (z) = z − u τ z ∈ [u τ + bh( − 1), u τ + bh( )], = 1, . . . , K, corresponding to the differentiable pieces in (10). The step functions v 3 (z) = bh( ) when z ∈ (u τ + bh( ), u τ +1 + bh( )), = 0, . . . , K. Solving for b in the equation ν(b) = 1 is equivalent to solving for b in the following equation, where F z is the cumulative distribution function and f z the density function of the adjusted residuals. In practice, a grid search is performed to approximate the solution b t . For each b in the grid, we use the empirical cumulative distribution, that is, )/h} with φ being the standard normal density function. The solution b t is taken to be the average of the smallest b in the grid that makes the righthand side of (13) smaller than s n and the next value in the grid.
Step 3: Update the estimator of β. Use the estimated b (t) from the previous step to update the estimate of β to The estimator β (t) , before applying the soft-thresholding function, is of interest too since it can be interpreted as a debiased estimator (Javanmard and Montanari, 2014a,b;Van de Geer et al., 2014;Javanmard et al., 2018); a thorough study of which, however, is beyond the current work.

State evolution
Within each iteration step t of the approximate message passing algorithm, state evolution studies the limiting behaviour of the estimators when the sample size goes to infinity. We now define the state evolution parameterζ 2 (t) which is critical for Algorithm 1. We start by defining the empirical version as follows Algorithm 1: RAMP algorithm for a single loss function with tuning parameter α 1 Function singleRAMP(α): empirical state evolutionζ 2 (0) using (15), tuning parameter of the soft-thresholding function θ (0) = αζ (0) 2 while iteration t ≤ T and tolerance tol > ε tol do 1. Adjust residuals: adjust the residuals z (t) ∈ R n : 2. Effective score: (b) update the state evolution parameterζ 2 (t) using (15) (c) update the tuning parameter θ (t) ← αζ (t) .
This quantity is linked to the state evolution recursion which describes the limiting behaviour of large systems, see Theorem 2 in Bayati and Montanari (2011a) and Lemma 1 in Bradic (2016). It holds that, see Section B.2.3 for details, whereζ (t) is the state evolution parameter for the large system, Z is a random variable with standard normal distribution independent of everything else andσ(t) is defined in (18). Due to the symmetry of Z, the state evolution parameter is formally defined as This definition explicitly features the extra Gaussian componentσ (t) Z in the limiting version, with varianceσ with B 0 defined in (A2). To connect the theoretical expression ofσ 2 (t) to Algorithm 1, we apply Eq.(3.6) in Bayati and Montanari (2011a), and Eqs.(7.10) and (7.19) in Bradic (2016). This leads to Note that (19) features the debiased estimator from (14). We now explain the connection between the estimators that explicitly use an l 1 -regularization and the corresponding estimators from the RAMP algorithm.
By applying Theorem 2 of (Bradic, 2016), we get the immediate connection between the regularized estimators β k (λ k ) for k = 1, . . . , K and the corresponding estimators obtained by applying the RAMP algorithm with a suitable choice of its regularity parameter α. We explain this below. Since the regularized estimators β k (λ k ) for k = 1, . . . , K are used for β MA (λ), (2), the connection between the model-averaged estimators from regularization and from application of the RAMP algorithm, follows immediately from the connections between the K separate estimators. The composite estimator β C (λ), (3), is a special case of a model-averaged estimator with K = 1, weight equal to one, and loss function ρ C = K k=1 w k ρ k . Denote (ζ 2 , b) as the fixed point solution when the iteration number t → ∞ of the following equations, Note that (20) is the state evolution recursion for the large system while in (21) the first equality is the population version of the requirement in step 2 in Algorithm 1 which states that n −1 n i=1 ∂ 1 G(z i,(t) ; b (t) ) = 1. The second equalities of both (20) and (21) follow by using the definition of G in (11), with G being defined in (8).
Then, under assumptions (A1)-(A5) (see the Appendix), for the RAMP algorithm with θ = αζ, where the tuning parameter α > 0 (which motivates the definition of θ (t) = αζ (t) in Algorithm 1), and for the l 1 -optimization with it follows by Theorem 2 of Bradic (2016) that The convergence in (23) explicitly connects the two composite estimators: one estimator uses an explicit l 1 -regularization as in (3), the other estimator is obtained via the RAMP algorithm. Similar results can be found in Huang (2020, Theorem 2.2) for a generalized AMP algorithm with non-negative convex loss function, and in Bayati and Montanari (2011b, Theorem 1.8) for the AMP algorithm with least squares loss function.
For the model averaging estimator we use such an equivalence for estimation with each separate loss function ρ k , k = 1, . . . , K. When using explicit l 1 -regularization as in (1) with the regularization constants λ k matching as in (22) the values θ k = α kζ , for k = 1, . . . , K that are used in the RAMP algorithm, again Theorem 2 of Bradic (2016) applies. It hence follows that

Theoretical contributions
This section contains detailed theoretical developments for the composite as well as the model-averaged AMP estimators in high-dimensions.
In Section 4, we defined the empirical state evolution parameterζ 2 emp,(t) , and we described the connections between the empirical updates in Algorithm 1 and the theoretical state evolution recursion, which connects to the theoretical expression of the AMSE. While Algorithm 1 and the theoretical state evolution recursion involve only a single estimator, the model-averaged estimator, on the other hand, is the weighted sum of K such estimators β k , k = 1, . . . , K, each obtained by Algorithm 1. Consequently, the estimators β k , k = 1, . . . , K are correlated.
Lemma 1 extends Theorem 2 in Bayati and Montanari (2011a) and (3.16) in Lemma 1(b) in Bayati and Montanari (2011a) to the almost sure convergence of the product for any two recursions among K paralleled recursions. All proofs are contained in Appendix B.2. (17) and (18). These are the state evolution parameters for the kth estimation with initializationσ 2 k,(0) = lim n→∞ 1 n n i=1 q 2 (0),i /δ. Then Lemma 1 in Bayati and Montanari (2011a) holds individually for each of the K recursions; additionally, for all pseudo-Lipschitz functionsψ c : R t+2 → R of order κ c for some 1 ≤ κ c ≤ κ/2 with κ as in (A4) and t a natural number larger than or equal to 0,

Lemma 1. Let the sequences of design matrices
is a bivariate standard normal vector with covariance not necessarily equal to zero.
Note that Algorithm 1 belongs to the general recursion in Bayati and Montanari (2011a), the initial condition takes q (0) = −β and the kth estimator calculated by Al- We obtain at iteration t, for k 1 , k 2 ∈ {1, . . . , K}, where Z k 1 and Z k 2 are possibly dependent standard normal random variables.
Since the estimators β kr , r = 1, 2 use the same design matrix, a correlation between Z k 1 and Z k 2 exists (see Corollary 2) and contributes to the correlation between β k 1 and β k 2 . Using Lemma 1, we obtain the theoretical AMSE for the regularized model-averaged estimator.
Theorem 1. Assume conditions (A1)-(A5) in the Appendix. At Algorithm 1's iteration step t for the estimator β k,(t) , for each k = 1, . . . , K, and for a weight vector w = (w 1 , . . . , w K ) , the model-averaged estimator β MA,(t) = K k=1 w k β k,(t) has asymptotic mean squared error Since the AMSE expression of the regularized model-averaged estimator is a quadratic function of the weight vector w, Corollary 1 readily provides the lower bound of the AMSE as well as the weight vector reaching this lower bound. The K-vector 1 K consists of ones only.

Corollary 1.
Constraining the weights to sum to one, the lower bound of the AMSE at iteration t for the model-averaged estimator as in (26)

Estimating optimal weights
The expression of the core matrix Σ (t) , which is the limit matrix for n, p → ∞, contains the random variable B 0 which satisfies assumption (A2) in the Appendix. Likewise, Σ 0,(t) which is the limit matrix for fixed p while n → ∞, contains the true coefficient β (see (27)). In practice, neither the true coefficient vector β nor the random variable B 0 is known. To make practical use of the expressions of the AMSE, we derive an estimator of the matrix Σ 0,(t) relying only on sequences generated in Algorithm 1.

Model-averaged estimator
Before deriving the estimator of the AMSE for the model-averaged estimator, we first defineζ emp,(k 1 ,k 2 ),(t) which is an estimator of the parameterζ (k 1 ,k 2 ),(t) , a quantity similar to the state evolution parameterζ 2 k,(t) , which records the covariance between the unbiased sequences β k 1 ,(t) and β k 2 ,(t) generated in (14) in Algorithm 1 when p → ∞. Since modelaveraged estimators combine estimators constructed from the same data into one weighted average, the correlation between β k 1 and β k 2 is needed to understand the AMSE of the model-averaged estimator.
Notice that the unbiasedness of the sequence β k,(t) follows from the argument that β k,j,(t) converges weakly to B 0 +ζ k,(t) Z k when p → ∞, while assigning 1/p point mass to each entry of the vector. Then, k,(t) ) for large p, indicating that β k,j,(t) centers at β j ensuring the unbiasedness. Moreover, the vector β k,(t) has Gaussian distribution. By applying the soft-thresholding function η on β k,j,(t) in Lemma 4, we avoid the usage of the true coefficient vector β in Σ 0,(t) resulting in a Stein-type risk estimator requiring only observables from Algorithm 1. A Gaussianity argument has also been used in Bayati and Montanari (2011b); Bayati et al. (2013); Mousavi et al. (2013Mousavi et al. ( , 2018 to derive a similar Stein-type risk estimator for the Lasso. Details can be found in Section B.2.8. The bias of the estimator β k,(t) is introduced in Algorithm 1 by applying the soft-thresholding function componentwise to the unbiased sequence β k,(t) .
Corollary 2. Assume conditions (A1)-(A5) in the Appendix. For any k 1 , k 2 = 1, . . . , K, at iteration t, whereζ k,(t) , k = k 1 , k 2 is the state evolution parameter corresponding to β k . Corollary 2 indicates both the existence and a feasible estimation of the covariance between Z k 1 and Z k 2 . As an estimator for We now state an unbiased estimator for the matrix Σ 0,(t) , and a consistent estimator for the matrix Σ (t) upon convergence of Algorithm 1.
Theorem 2. Assume conditions (A1)-(A5) in the Appendix, and that the state evolution parameter in (15) This new estimator can be compared to the estimator used in Bayati et al. (2013, Def. 2) and Mousavi et al. (2018, Eq.(9)) for the case of a single estimator (K = 1). The proof of Theorem 2, see Section B.2.8 uses Stein's lemma (see Lemma 4) to estimate the covariances that appear in the matrix Σ 0,(t) . The soft-thresholding function η(·; θ) that appears in the estimator Σ 0,(t) links the estimator β k to the estimator β k . The proof also uses the joint asymptotic normality of the jth components of the vectors β k 1 and β k 2 . The obtained estimator for Σ 0,(t) in the case K > 1 is nontrivial and new to the literature.
Estimated AMSE-optimal weights for the model-averaged estimator are obtained by using the estimator Σ 0,(t) at the final iteration in Theorem 2. In combination with the sumto-one constrained weights this gives the estimated weights that minimize the estimated AMSE for the model-averaged estimator When additional constraints such as positivity are needed, the optimal weights no longer have an explicit formula, but they are straightforward to compute, see (31).

Composite estimator
The AMSE of a composite estimator can be obtained from Theorem 1 as a special case, treating the composite loss function as a single loss function with weight one, thus ρ C = K k=1 w k ρ k as in (3). At iteration t, The matrices Σ (t) , Σ 0,(t) are now real numbers and coincide with the AMSE of the estimator in (25). We obtain the corresponding estimator for the AMSE For the single loss function, ρ C , the estimator of AMSE in (29) can be compared to the Stein-type estimator that has been obtained in Definition 2 in Bayati et al. (2013) for the AMP algorithm using the least squares loss, which is a particular case of Algorithm 1.
Finding optimal weights for the composite estimator is complicated. Indeed, while the model-averaged estimator has an AMSE, which is a quadratic function in the weights, see (26), the composite estimator and its AMSE depend on the weights in a highly nonlinear fashion; e.g., observe that the soft-thresholding function in (29) depends on w.
Therefore, optimization of the estimated AMSE with respect to the weights proceeds numerically; See Section 6.2 for more details.

5.3
The case of dense (non-sparse) linear models with n/p → δ ≥ 1 : asymptotic variance optimality Donoho and Montanari (2016) and El Karoui et al. (2013) showed that the asymptotic variance of the M-estimators in the case where p, n → ∞ and n/p → δ ∈ [1, ∞) contains an extra Gaussian component. Recently, Lei et al. (2018) obtained the coordinate-wise asymptotic normality of regression M-estimators in the moderate p/n regime for a fixed design matrix. In the sparse high-dimensional linear model setting where δ ∈ (0, 1), it was shown that the sequence β (t) in (14) follows for the Lasso estimator (Bayati et al., 2013) a similar normal distribution with the variance containing an extra Gaussian component. The above-mentioned literature focuses on the asymptotics for a single M-estimator; we extend the asymptotic result to the model-averaged estimator. In this section, we only characterize the asymptotic variance of the model-averaged estimator for dense linear models with n/p → δ ≥ 1, following Donoho and Montanari (2016).
Under the dense linear model with n ≥ p, the soft-thresholding function η(·; θ) is replaced by the identity function and the ratio ω = E[ B 0 0 ] = 1. Consequently, Algorithm 1 is adjusted to estimate where β is dense. It is trivial to show that Algorithm 1 still belongs to the general recursion in Bayati and Montanari (2011a). For a single estimator at iteration t denoted as β k,(t) , the two state evolution parametersζ 2 k,(t) andσ 2 k,(t) coincide and Theorem 4.1 in Donoho and Montanari (2016) holds.
Theorem 3. Assume conditions (A1)-(A5) in the Appendix. Let n/p → δ ≥ 1 when n, p → ∞. For the asymptotic variance of the model-averaged estimator β MA holds that The extra Gaussian component is identified in the convolution of the regression noise distribution and a Gaussian distribution: Recall that the componentwise empirical distribution of β k (p), when p → ∞, converges weakly to B 0 +ζ k Z k following Bayati and Montanari (2011a) and Donoho and Montanari (2016). Then for large p, while the iteration t → ∞, β k (p) ∼ N (β,ζ 2 k I p ) (Donoho and Montanari, 2016; Mousavi et al., 2013) Note that this estimator coincides with (28) for the special case that n ≥ p and the softthresholding function is replaced by the identity function.
6 Computational details 6.1 Regularized model-averaged quantile estimation The estimation of the quantile u τ k = F −1 ε (τ k ) follows a two-step procedure.
1. Obtain an initial slope estimate β init and calculate the residuals. Example initial slope estimates are the Lasso or regularized quantile estimation with a single quantile level.
2. For k = 1, . . . , K, estimate the quantile intercepts u τ k by taking the corresponding τ k × 100% quantile of the residuals from the previous step.
The regularized model-averaged estimator is obtained by averaging over K paralleled estimators. See Algorithm 2 for the pseudo-code, of which the core is Algorithm 1; there the effective score function G is that of a single quantile loss function with K = 1, see also Example 2 in Bradic (2016). In our numerical work, the upper bound for the number of iteration steps T is set to be 50 in both the simulation and the data analysis sections. With K = 1, this algorithm applies to the regularized composite estimator too.
AMSE refers to the estimated version. The β k s are recorded for calculating the weights in Corollary 1.
The tuning parameter α of Algorithm 2 controls the sparsity of the estimators and requires a tuning procedure to choose it in practice. In Section 7, we consider the one dimensional Golden-section search algorithm (Kiefer, 1953) for tuning the value α in the range [α min , α max ] that minimize the estimated MSE of β using the estimator derived in Section 5.2. The upper bound α max is chosen to be 2.3 for the simulations and data analysis. The lower bound α min in the data analysis follows the lower bound in Proposition 9. 2 in Eldar and Kutyniok (2012) and is chosen to be the unique non-negative solution to the equation (1 + α 2 )Φ(−α) − αφ(α) = δ/2, where φ(x) and Φ(x) denote the p.d.f and c.d.f of the standard normal distribution respectively. In the simulation study, the lower bound α min is chosen to be 1.3 for computational efficiency purposes, since the optimal tuning parameter for those settings was rarely less than 1.3.
2. Estimate the AMSE-type optimal weight w MA,1 with constraints by where the K × K matrix Σ 0 is the consistent estimator of Theorem 2.
It is worth mentioning that w MA,1 is a constrained version of w MA attaining the lower bound of the AMSE in Corollary 1. w MA,1 focuses on approximating the lower bound of the AMSE of the sparse coefficient vector β without assuming that the nonzero entries are selected perfectly; whereas another type of weight choice derived in Bradic et al. (2011); Bloznelis et al. (2019) aims at the lower bound of the variance of the nonzero part of β by imposing the perfect selection assumption. A numerical comparison of these two types of weight choices is presented in Section 7.
To equip the regularized composite quantile estimator with the weight minimizing the estimated AMSE, we cannot make use of an analytical solution to the weight minimization problem. Instead, a numerical search for a better weight choice in the neighbourhood of an initial weight proposal is employed. The basic idea is that the estimator β C (w C ) is treated as a function of the weights. We propose a collection of candidate weight vectors in the neighborhood of the weight chosen in the previous step. The weight for β C (w C ) is updated in each step by the one having the lowest estimated AMSE, i.e., w C,1 = arg min w cand AMSE( β C (α opt ; w cand ); β).
A more detailed search procedure is as follows.
1. Propose a reasonable initial weight vector w C,init , e.g. the vector of equal weights; estimate β C at the initial weight w C,init and obtain the estimate of AMSE β C (α opt ; w C,init ); β .
2. Initiate the searching step calculator s D = 0, the candidate optimal weight w C,1 = w C,init , and the corresponding candidate minimum MSE AMSE(w C,1 ) = AMSE β C (α opt ; w C,init ); β estimated by the AMSE estimator in Theorem 2 for K = 1, the collection of the used weight vectors V w = {w C,1 }.
3. Propose a set of candidate weight vectors V w cand . This is to exclude those recorded in the collection of the used weight vectors V w . In addition, V w cand should be in the neighborhood of the current optimal weight w C,1 . Rules of proposing candidate weight vectors are user-decided; here, we consider a (K − 1)-dimensional grid search centering at w C,1 .
4. Obtain the regularized composite quantile estimates at all candidate weight vectors in V w cand with Algorithm 2. Update the used weight vector collection V w , increase the counter s V = s V + 1, update the candidate optimal weight w C,1 by the weight with the lowest estimated AMSE in V w = {w C,1 }, and update the candidate minimum AMSE value AMSE(w C,1 ).
5. Stop the iteration if the searching step calculator s V > S V or the candidate weight vector collection V w cand = ∅; otherwise repeat steps 3 and 4.
The pseudocode of the search procedure is stated in Algorithm 3.
2 while searching step s V ≤ S V or candidate weight collection V w cand = ∅ do 1. Propose a new V w cand in the neighbourhood of w C,1 . Rules of proposing candidate weight vectors are user-decided; here, we consider a (K − 1)-dimensional grid search centering at w C,1 .
7 Numerical results

Simulation study
In this section, we consider the following setup under the high-dimensional linear model setting.
1. Fix the dimension p = 500, the sample size n = 250, the ratio δ = 0.5. The number of non-zero components s is taken to be 5 for the high-sparsity setting and 50 for the medium-sparsity setting; the non-zero part is generated from the Dirac distribution with a point mass equally distributed on -1 and 1, or a standard normal distribution.
2. In each repetition, we generate a new dataset by randomly generating a sensing matrix X, a coefficient vector β, and an error vector ε. The components of the sensing matrix X are independent and generated from N (0, 1/250).
The objective is to compare the performance of the regularized model-averaged estimator and the composite estimator with different weights, with emphasis on the weights where the selection uncertainty is taken into account. The simulation is repeated to get 500 estimates for each setup. For both the regularized model-averaged and composite quantile estimator, the weights considered are (1) the estimated AMSE-type weights (i.e. w MA,1 for the model-averaged quantile estimator and w C,1 for the composite quantile estimator), (2) the estimated weights based on minimising the asymptotic variance of the estimators of only the active set of coefficients, denoted by w MA,2 (Bloznelis et al., 2019) and w C,2 (Bradic et al., 2011) where, with the (k 1 , k 2 )th component of A equal and w C,2 = arg min w,a ε w=1,w k ≥0 w Aw .
Only considering the variance has been the standard practice so far. (3) Equal weights 1/K for each component. The number of quantiles K for both estimators is taken to be 3, with quantile levels 25%, 50%, 75%.
We present the empirical MSEs of the abovementioned estimators for estimation of three vectors of coefficients. First, we consider the estimator of the subvector of the full coefficient that consists of only the non-zero true coefficients, we refer to this as the "nonzero part". Second, we consider the estimator of the subvector of the coefficients that are truly zero. This is referred to as the "zero part". Third, we consider the full vector of estimated coefficients. Note that some truly zero coefficients might have a non-zero estimate, while some truly non-zero coefficients might be estimated as zero. For each of these three vectors, "parts", we compare the estimated values with the true values to get MSE( β part ) = ppart jpart=1 ( β jpart − β jpart ) 2 /p part for the appropriate part of the full vectors. Results for the regularized model-averaged quantile estimator with different weights are presented in Table 1. We observe that the model-averaged quantile estimator using the weight in (31) has lower MSEs for estimating the non-zero part of β and for the full vector β, and this for t 3 and the mixture of normally distributed errors in the high-sparse case where the number of non-zero components s = 5. Using equal weights leads to a fair performance of the model-averaged quantile estimator, especially for estimating the all-zero part of β. The Lasso estimator is considered as the baseline comparison, which from Table 1 seems to have a competitive performance, especially in the medium sparsity settings. However, the Lasso mostly gives over-sparse estimations, which can be observed in the top half of Table 2 summarizing the averaged true positive (TP) and true negative (TN) recovery rates which are defined as TP (TN) = number of correctly identified as non-zeros (zeros) number of true non-zeros (zeros) The Lasso has the highest TN rate consistently and mostly the lowest TP rate. Further, while increasing the standard deviation of the errors, the Lasso's overly-sparse estimation becomes clearer, i.e., Lasso gives sparser estimations and becomes all-zeros eventually. The regularized model-averaged estimator with equal weights mostly has the highest TP rate, except for the medium sparsity settings where the non-zero part of the true regression coefficient is sampled from a Dirac distribution at -1 and 1, and the errors are sampled from N (0, 1) or 0.5N (0, 1) + 0. 5N (5, 9). The model-averaged estimator with the weight in (31) has the second-highest TN rate consistently.
Since there is no analytical expression for the selection incorporated weight of the regularized composite quantile estimator w C,1 , the choice of weights can only be determined numerically by an exhaustive search. To reduce the searching time of the composite quantile estimator, we set the stopping criterion S V to be five and only randomly select 4 points in the neighborhood V w cand ; the tuning parameter α of the soft-thresholding function is tuned once for the regularized composite quantile estimator with the weight w C,2 , then fixed after that. Table 3 summarizes the empirical MSEs of the regularized composite quantile estimator with different weights. Since the tuning parameter, α is selected for w C,2 and a fixed tuning   parameter is used for obtaining the regularized composite quantile estimates with other weights, it is not surprising that using w C,2 leads to lower MSEs in most cases. However, it is worth noticing that using equal weights, while α is not optimally tuned, leads to the regularized composite quantile estimator's fair performances. The Lasso estimator consistently has the lowest empirical MSEs recovering the all-zero parts, through the largest empirical MSEs recovering the non-zero parts. This is caused by overly sparse estimations of the Lasso, which is indicated in the bottom half of Table 2. The regularized composite estimator with locally optimized w C,1 consistently has the highest TP rate, and secondhighest TN rate among all competitors, except the TN rate for t 3 distributed errors and TP rate for 0.5N (0, 1)+0.5N (5, 9) distributed errors. At the same time, the non-zero parts of β are generated from Dirac distribution at -1 and 1. Tables 3 and 2 illustrate that the regularized composite quantile estimator mostly improves the performance of regularized single quantile estimator. For the same simulations settings, we compare the averaged empirical MSEs, true positive and true negative rates of the regularized composite quantile estimator, see Table 2, column 7 and 12, and Table  3, column 7, with the single regularized quantile estimator at the median τ = 0.5. For settings where s = 5, the composite quantile estimator clearly dominates the single quantile estimator for all three error distributions. For settings where s = 50, the composite estimator still mostly outperforms the single quantile estimator, except for the following cases: (1) the MSE for the non-zero and zero estimated subvector of β in settings where errors are generated from N (0, 1) and 0.5N (0, 1) + 0. 5N (5, 9) distribution and the true non-zero subvector of β is generated from a Dirac distribution; (2) TN rates in settings where errors are generated from N (0, 1) and t 3 distribution and the true non-zero subvector of β is generated from N (0, 1).
The percentage of converged cases for the model-averaged and composite estimator, while setting the tolerance ε tol to be 10 −6 for different error distributions, are included in Table 4, where we define a estimator to have converged when the needed number of iterations was less than 50.
Condition (A1) restricts Algorithm 1 to a special design matrix that does not allow correlations between the X ·j 's. However, since such correlation might be present in reality, it is of interest to see if Algorithm 1 is still numerically robust while Condition (A1) is relaxed in practice. We consider a similar simulation setup as used before with p = 500, the sample size n = 250, and δ = 0.5. The number of non-zero components s is taken to be 5 or 50; the non-zero components are generated from the Dirac distribution with point mass equally distributed on -1 or 1, or a standard normal distribution. In each simulation replication, a design matrix is first generated from a multivariate Gaussian distribution N (0, Σ X ), then the components X i,j are centered and scaled such that the components of the rescaled matrix X have sample variance 1/n. Here, we allow for a Toeplitz covariance matrix Σ X of which its (i, j)th component (Σ X ) i,j = σ |i−j| X , i, j = 1, . . . , p. We consider σ X = 0, 0.1, 0.3.   5N (5, 9) 77 86 71 85 Table 4: Percentage of converged cases of both regularized model-averaged and composite quantile estimators with the convergence tolerance ε tol = 10 −6 . The convergence percentage of the regularized model-averaged estimator is calculated by including only those cases of which all single quantile component estimates converge in less than 50 iterations.
To investigate the effect of the correlation on the RAMP algorithm we consider the regularized single quantile estimator at quantile level 0.5. The error distribution considered is t 3 . Table 5 records the performance of Algorithm 1 with tolerance ε tol = 10 −6 for such a correlated design matrix; the performance is evaluated by the empirical MSEs, the TP and TN rates, and the percentage of convergence. We see from Table 5 that parameter estimation using Algorithm 1 remains accurate and stable when weak correlations such as with σ X = 0.1 exist between the X ·j 's; the accuracy drops when we further increase the correlations as with σ X = 0.3; it is worth mentioning that the convergence percentages decrease when the correlation increases. Further research concerning correlated data is worth considering.

Data analysis
We consider the audio wave file of a waveshape from Octave in the R package signal. The dataset is a list of 3 elements; the audio wave sample is a vector of 17380 entries stored in the element "sound", the sample rate is 22050 Hz stored in the element "rate", and the resolution of the wave file is 16 bits recorded in the element "bits". To alleviate the computational burden of the signal compression and reconstruction, we only consider the signal from the 6145th entry to the 8192th entry of the original sound wave signal.

The preprocessing -discrete wavelet transform
Originated from the compressed sensing problem, the sparse linear model Y = Xβ + ε describes the image or signal compression. The s-sparse p-dimensional input signal β is first compressed by a known sensing matrix X ∈ R n×p with n < p; the compressed signal vector Xβ ∈ R n can be corrupted by the noise ε with ε i 's i.i.d. via transmission. Notice that the p-dimensional input signal vector β is assumed to be s-sparse which is usually unsatisfied by signals expressed in the standard basis. To obtain the sparse representation of β in practice, an intermediate stage of expressing the natural non-sparse vector β * in a  Table 5: The mean, over 500 simulation repetitions, of the empirical MSEs, the true positive (TP), the true negative (TN), and convergence percentages of the regularized modelaveraged quantile estimators for t 3 distributed errors. Empirical MSEs are calculated for the non-zero parts, all-zero parts, and the full vector of the true coefficient β. The nonzero part of the true coefficient vector is generated from Dirac distribution with point mass equally distributed on -1 and 1 (top), or standard normal distribution (bottom).
proper orthonormal basis Ψ * = (ψ * 1 , . . . , ψ * p ) is required. Examples of such an orthonormal basis include the orthonormal wavelet basis, the Fourier basis, and so forth. To perform the discrete wavelet transform, we use the R package wavethresh. The collection of the coefficients at all resolution levels is used for further compression.

The artificially corrupted compression
To imitate the compressed sensing process, we process the audio wave signal vector as follows: 1. Perform the Daubechies' least asymmetric wavelet transform with 8 vanishing moments using the wd function in the R package wavethresh on the original signal β * ∈ R 2048 and obtain the corresponding wavelet coefficient vector β ∈ R 2047 with p = 2047.
2. Randomly generate the sensing matrix X with i.i.d components X ij ∼ N (0, 1/n), where n = δ p and δ is the undersampling ratio chosen to be 0.5 here; compress the corresponding wavelet coefficients β by computing Xβ.
3. Corrupt the compressed wavelet coefficients by the error vector ε with i.i.d. components ε i having p.d.f f ε ; obtain the artificial observed signal vector Y = Xβ + ε. Additionally, the standard normal N (0, 1), student-t with 3 degrees of freedom, and the bimodal mixed normal 0.5N (0, 1)+0. 5N (5, 9) are used as the corruption error distributions; the errors are sampled according to the distributions first, then centered and rescaled to have standard deviation 0.03.
In practice, the artificial vector Y and the sensing matrix X are observed. The accurate recovery of the original wavelet coefficient vector β is of practical interest. To obtain an impression on the performance of the AMSE-type optimal weight, we generate the sensing matrix X under a fixed seed number, which is set to be 1 in our case, then generate the error vector ε under various seed numbers. However, we only present the reconstructions under one seed for each setting in Section 7.2.3 due to limited space.

Signal recovery
To reconstruct the signal vector β expressed in the wavelet basis from the sensing matrix X and the observed compressed signal vector Y corrupted by potentially non-Gaussian distributed error ε, we consider the regularized model-averaged and the composite quantile estimator weighting over three equally-spaced quantiles (25%, 50%, 75%) using equal weights, the oracle-type weights and the new AMSE-type weights. The tolerance in the RAMP algorithm is set as ε tol = 10 −8 . The Lasso estimator is considered as the baseline comparison. Notice that the regularized estimates β MA and β C after reconstruction are the representations in the wavelet domain. To compare the accuracy of the reconstruction, Reconstructed audio signal fractions using regularized MAQR with K = 3 and the Lasso Figure 2: Reconstructed audio signal from using the regularized model-averaged estimator with the estimated AMSE-type weights in (31), oracle-type optimal weights in (32), and equal weights. The original audio curve is depicted in black. The Lasso reconstruction is presented at bottom-right. The error used for corruption follows the mixture of normals distribution 0.5N (0, 1) + 0.5N (5, 9).
we perform a back-transform on the estimates and obtain the corresponding signal vectors β * MA and β * C with representations in the natural basis. Example reconstructions of the audio signal for K = 3 using the regularized modelaveraged estimator equipped with different weights, with the baseline recovery from the Lasso represented in the natural basis are presented in Figure 2 for the mixture of normals distributed error, and in Figure 3 for the t 3 distributed error. We observe that the strong signals corresponding to large values located at the end of the sound signal are well captured by the model-averaged quantile estimator using different weights for both error distributions. For the weak signals clustering at the front of the signal, the modelaveraged estimators using w MA,1 and equal weights outperform the counterpart with w MA,2 for 0.5N (0, 1) + 0.5N (5, 9) distributed errors; recovery differences for the weak signals of Reconstructed audio signal fractions using regularized MAQR with K = 3 and the Lasso Figure 3: Reconstructed audio signal using the regularized model-averaged estimator with the estimated AMSE-type weights in (31), oracle-type optimal weights in (32), and equal weights. The original audio curve is depicted in black. The Lasso reconstruction is presented at bottom-right. The error used for corruption is t 3 distributed.
the model-averaged estimator using different weights are hardly observable for the t 3 distributed errors. Recovery using the Lasso is competitive to the model-averaged estimator using w MA,1 for strong signals. However, the Lasso estimates the signals in an over-sparse way with too many zeros entries; one can observe the almost flat recovery for the weak signals for both error distributions. Bates and Granger (1969) provide an alternative weight choice for the model-averaged estimator obtained by considering only the variances of β k 's and ignoring the covariances. This leads to w MA,3 = arg min where diag( Σ 0,(t) ) denotes the diagonal matrix obtained from Σ 0,(t) which keeps the diagonal and has zeros in all off-diagonal entries. Figure 4 contains the recovery of the audio signal using the model-averaged estimator using this weight. For the composite quantile estimator β C , we performed the same weight searching method as for the simulation study. This is, S V = 5 and randomly select 4 candidate weights in the neighbourhood of the previous value. We select the tuning parameter α once for the starting weight w C,2 , it remains unchanged thereafter. The recovered signals by the composite estimator with different weights are very similar in all cases.
To compare the recovery of the regularized model-averaged and composite estimator combined with different weights, as well as the Lasso estimator, we present the mean absolute percentage error (MAPE) in Table 6 where the MAPE is defined as Table 7 reports the MSE. We see that the Lasso has the lowest MAPE for both t 3 and mixed normal distributed errors; at the same time, it estimates the weak signals in an over-sparse way and is not capable of capturing the weak signals. Comparing the effect of different weight choices on the regularized model-averaged quantile estimator with its composite quantile counterpart, we see that the MAPEs of the composite quantile estimators are relatively stable using different weights. The model-averaged estimator with the AMSE-type weight w MA,1 f ε t 3 0.5N(0,1) + 0.5 N(5, 9) est: MA / C w est,1 w est,2 w eq w est,3 w est,1 w est,2 w eq w est , 5N (5, 9) w est,2 (0.147, 0.843, 0.010) (0.369, 0.345, 0.286) Table 8: The estimated weights w MA,1 and w MA,2 for the model-averaged estimator, and w C,1 and w C,2 for the composite estimator. The seed number used to generated the errors for corrupting the compressed signal vector is 37 for both t 3 and mixed normal distributed errors.
has excellent performance compared to the composite estimator, especially for the mixed normal distributed error. The Bates-Granger weighting provides good results regarding MAPE for the t 3 error case, but not for the mixed normal. Regarding MSE, it performs well for the mixed normal case but is worst for the t 3 errors, wherein this example the equal weights perform best, although all results are close. Searching for the selection incorporated weight w C,1 for the regularized composite quantile estimator is computationally infeasible for large p (2047 in our case). Estimating the regularized model-averaged quantile estimator averaging three quantiles here takes approximately 4 -5 hours whereas estimating the regularized composite quantile estimator takes more than 16 hours with only five steps in a nearby search with four surrounding candidate weights, and the tuning parameter α tuned only once for the starting weight. Additionally, we present the estimated weights for both regularized model-averaged and composite estimators in Table 8. An interesting observation is made by comparing the estimated weights w MA,1 and w MA,2 for the mixed normal distributed error. The weight w MA,1 presented here is quite representative; it assigns weight 0 to the quantile estimate at 50% quantile level suggesting the final model-averaged estimate is obtained by averaging estimates at 25% and 75% quantile levels. On the contrary, w MA,2 assigns the largest weight to the estimate at a 50% quantile level indicating the most significant contribution to the final model-averaged estimate.

Discussion
This paper is the first to take the selection uncertainty due to regularization into account when computing the weights used in model-averaged and composite estimation. While we have studied both composite estimation and model-averaged estimation, the flexibility of allowing for parallel computation and a component-specific choice of regularization, combined with an explicit expression of the optimal weights for model averaging, places this method in a preferred position from a computational point of view.
It would be interesting to investigate whether AMSE expressions for other types of regularization may be obtained similarly. Going yet one step further would be incorporating the effect of data-driven values of the regularization parameters λ (for composite estimation) and λ 1 , . . . , λ K (for model-averaged estimation) on the choice of the weights. To further study the weight selection and the effect of using data-driven weights, one should study the joint distribution of the estimated weights and the estimators of interest. To simplify such matters, sample splitting could be used such that the weights are computed on a hold-out sample and the estimation using those weights proceeds on the rest of the sample. In this paper, we used the same dataset for estimating both β and w.
To avoid overly complicated mathematical expressions, we followed earlier literature in the use of a design matrix where X ij ∼ N (0, 1/n). Other applications might require studying, for example, fixed designs, which are beyond the scope of the current paper.

A Assumptions
(A1) Design: The elements of the design matrix X, that is X ij for i = 1, . . . , p and j = 1, . . . , n, are independent and identically distributed according to a N (0, 1/n) which is also called a standard Gaussian design.
(A2) Coefficients: The p-vector β is such that the sequence of uniform distributions that is placed on its components converges, for p tending to infinity, to a distribution with a bounded (2k − 2)th moment for k ≥ 2. Denote by B 0 a random variable with this limiting distribution function F B 0 .
(A3) Loss function: (i) The subgradient ∂ρ(u) = 3 j=1 v j (u) where v 1 has an absolutely continuous derivative, v 2 is continuous and consists of piecewise linear parts and is constant outside a bounded interval, and v 3 is a non-decreasing step function.

B.1 Auxiliary definitions and lemmas
It follows that if φ is a pseudo-Lipschitz function of order κ, then there exists a constant L such that ∀x ∈ R m : |φ(x)| ≤ L (1 + x κ ).
Lemma 2 (Theorem 1 in Jameson (2014)). If x i ≥ 0 where i = 1, . . . , n and p ≥ 1, then The reversed inequality holds for p ∈ (0, 1) Lemma 3 (Extrema of quadratic forms in Rao (1973)). Let A be a m × m matrix, B be a m × k matrix, and U be a k-vector. Denote by S − any generalized inverse of B A −1 B, then where X is a column vector and the infimum is attained at A −1 BS − U .
Lemma 5 (Lemma 4 in Bayati and Montanari (2011a)). Let κ ≥ 2 and a sequence of vectors {β(p)} p≥0 whose empirical distribution converges weakly to probability measure f B 0 on R with bounded κth moment; additionally, assume that Then for any pseudo-Lipschitz function ψ : R → R of order κ:

B.2.4 Estimation of ν(b)
The effective score step in Section 3.2, in cases where G(·; b (t) ) is non-differentiable, requires a solution b (t) to the equation 1 = ν(b (t) ) where ν(b (t) ) is a consistent estimator of a population parameter ν(b (t) ) defined as with C (t) = ε −σ (t) Z the random variable characterizing the limit distribution of the adjusted residuals z (t) when p → ∞.
Using Condition (A3) and Lemma 3 of Bradic (2016), ∂ρ can be written as a sum of three functions of which v 1 and v 2 are differentiable. For the step function v 3 , we use Condition (A3) on ρ, where γ l is the step height on the interval (r l , denote the density of the variable C (t) − G(C (t) ; b) which is equivalent to Prox(C (t) ; b). The equivalence is obtained by setting the derivative of the bρ(x) + 1 2 (x − C (t) ) 2 w.r.t. x to zero and evaluate at Prox(C (t) ; b), due to the fact that the proximal operator is the minimizer of the function bρ(x) + 1 2 (x − C (t) ) 2 . Then we arrive at The consistent estimator in (12) is obtained by replacing the expectation above with the empirical mean and replacing the density of the proximal operator Prox(C (t) ; b) with its kernel density estimator.

B.2.5 Proof of Lemma 1
Since this proof is based on the general recursion and Lemma 1 in Bayati and Montanari (2011a), we first restate the general recursion to which Algorithm 1 belongs with slight changes in the notations. Given the noise ε ∈ R n and the coefficient vector β ∈ R p , the general recursion is defined where ξ 1,(t) = n −1 n i=1 ∂ 1 g 1,(t) (d (t),i , ε i ), ξ 2,(t) = (δp) −1 p j=1 ∂ 1 g 2,(t) (h (t),j , β j ). Further, to connect the general recursion to Algorithm 1, we also state the exact form of h (t+1) , m (t) , d (t) , q (t) taken in Algorithm 1. Lemma 1 in Bradic (2016) states that Algorithm 1 takes To proceed with the proof of Lemma 1, we first recall the technique used for proving Lemma 1 in Bayati and Montanari (2011a), which uses induction on the iteration t. To not fully repeat the long proof and all notations we only give details about where our proof differs from theirs. .
For the first term in the last inequality above, we see that the expectation E[|L (1 + |[X k 2 m k 2 ,(0) ] j | κc +|β j | κc )| 2+ρ ] is bounded by some constant, since the expectation is with respect to the matrices X k 1 , X k 2 ,X k 1 ,X k 2 of which the components are Gaussian distributed with mean 0 and variance 1/n; the rest terms are bounded by a constant; the moments of Gaussian distributed r.v. are all finite. Let us denote the upper bound of this expectation by L , then the first term of the inequality above is bounded by , which can be shown to be bounded by cp ρ/2 following a similar argument as in Lemma 1(b) in Bayati and Montanari (2011a). The second term similarly can be shown to be bounded by c p ρ/2 . Hence the variable defined in (40) satisfies the condition in Theorem 3 in Bayati and Montanari (2011a); thus the a.s. convergence holds.
In the special case where k 1 = k 2 , we show that the square of ψ c is still pseudo-Lipschitz of order 2κ c ≤ κ, then the almost sure convergence hold by directly applying the result in Lemma 1 in Bayati and Montanari (2011a).