Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
112,158 characters · 23 sections · 67 citation commands
Detangling robustness in high dimensions: composite versus model-averaged estimation
We investigate the benefits of model-averaged as well as composite estimators in high-dimensional 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 HoetingMadiganRafteryVolinsky99, though we focus on its frequentist version in which a user determines the weights assigned to the separate estimators ClaeskensHjort08, HjortClaeskens03, YuanYang2005. Model averaging enjoys a wide application, see, for example, the recent overview paper for model averaging in ecology by Dormannetal2018 and for application to hydrology by HogeGuthkeNowak2019. In econometrics, the terminology “forecast combinations" appears ChengWangYang2015,BatesGranger1969; whereas “multimodel inference” is another commonly used term for this procedure BurnhamAnderson2002.
While the technique is quite thoroughly investigated for low-dimensional models, far fewer results have been obtained in high dimensions. ando2014model consider high-dimensional 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. ZhaoZhouLi2016 extend this method to dependent data, while ando2017weight 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 Hansen2007 and Hansen2012 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, Koenker2005quantile 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, BloznelisClaeskensZhou2019 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 BloznelisClaeskensZhou2019 or bradic2011penalized, 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 donoho2009message and bayati2011dynamics for compressed sensing. donoho2016high explain the use of AMP algorithms for obtaining the variance of high-dimensional M-estimators for which $n/p\to\delta\in(1,\infty)$. 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 bradic2016robustness 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 (ref), 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 (ref). The limiting behavior of the estimators in the algorithm is studied by state evolution parameters in Section (ref). We obtain the estimators' asymptotic mean squared error as well as an estimator of that quantity in Section (ref). We showcase the procedure for high-dimensional regularized quantile regression in Section (ref) and present numerical results in Section (ref). Section (ref) concludes. All proofs, together with the assumptions and some technical lemmas, are collected in the Appendix.
We consider a high-dimensional linear model $Y = X {\beta} + \bm{\varepsilon}$ with $Y \in \mathbbm{R}^n$, the design matrix $X \in \mathbbm{R}^{n \times p}$ and the parameter vector ${\beta} \in \mathbbm{R}^{p}$. The $i$th row of $X$ is denoted $X_{i\cdot}$, $i=1,\ldots,n$, the $j$th column of $X$ is denoted by $X_{\cdot j}$, $j=1,\ldots,p$. We assume the components of $\bm{\varepsilon}$ to be independent and identically distributed with mean zero, cumulative distribution function $F_\varepsilon$ and probability density function $f_\varepsilon$. We allow for a sparse high-dimensional setup. Denote by $s$ the $l_0$ norm of the parameter vector, $s=\|\beta \|_0$, which counts the number of nonzero components of the vector $\beta$. We assume that the ratios $n/p \to \delta \in (0,1)$ and $n/s \to a \in (1, \infty)$ when $p, n, s$ tend to $\infty$.
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 (ref). Second, composite estimation where a weighted average of loss functions is minimized; see (ref).
For model-averaged estimation of the parameter $\beta$, define for $k=1,\ldots,K$ the regularized estimators
where $\rho_1,\ldots,\rho_K$ are nonnegative convex loss functions and ${\bm\lambda}=(\lambda_1,\ldots,\lambda_K)^\top$ is a vector of possibly different nonnegative regularization parameters. For a set of weights $w=(w_1,\ldots,w_K)^\top$, the model-averaged estimator is defined as
Often one assumes that the weights $w_1,\ldots,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 $\lambda$, such that the regularized composite estimator is defined as
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 $\rho_k$ are nonnegative and convex. Hence, for composite estimation it is required that the weight vector $w \in [0,1]^K$ such that $\sum_{k=1}^K w_k=1$.
As a worked-out scenario throughout the paper, we consider quantile loss functions $\rho_k(\cdot), k = 1, \ldots, K$ that are defined below. For more information about quantile regression with i.i.d. errors, see Koenker2005quantile. In this paper, we assume that the design matrix $X$ does not contain a column of ones; see assumption (ref) in the Appendix. With $\tau \in (0,1)$, the $\tau$-quantile of the response $Y$ is obtained as $X\beta+F_\varepsilon^{-1}(\tau) = X\beta + u_\tau$.
Figure (ref) presents first a single quantile loss function with $\tau=0.3$, $$ \rho(x) = (x - u_{\tau}) (\tau - I\{x \leq u_{\tau}\}). $$ For model averaging we specify $K$ different quantile levels and use $K$ different such quantile loss functions for estimation of $\beta$: $$ \rho_k (x) = (x - u_{\tau_k}) (\tau_k - I\{x \leq u_{\tau_k}\}), \quad k \in \{1, \ldots, K\}. $$ For composite quantile estimation we assume that the $K$ quantile levels $\tau_1 < \cdots < \tau_K$, then also the quantiles of $\varepsilon$ are sorted $u_{\tau_1} < \cdots < u_{\tau_K}$. Define $u_{\tau_{0}}=-\infty$ and $u_{\tau_{K+1}}=\infty$.
The middle panel of Figure (ref) depicts such a composite quantile loss function $\rho_C=\sum_{k=1}^Kw_k\rho_k$ for $K=3$ quantile levels 0.25, 0.5 and 0.75 with equal weights $w=(1/3,1/3,1/3)^\top$. The panel on the right in Figure (ref) uses the same quantile levels but depicts the quantile loss function $\rho_C$ with weights $w=(0.15, 0.55, 0.3)^\top$.
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 bradic2016robustness. 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 $\widehat{\beta}_{{\rm MA}}(\bm{\lambda})$ and $\widehat{\beta}_{{\rm C}}(\lambda)$.
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 ((ref)) of the regularized loss function in the case of model averaging, and the estimator ((ref)) 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 ((ref)), respectively ((ref)), 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 \to \delta \in (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 (ref). Challenges arise with incorporating dependence into the AMSE expression; see Theorem (ref). Theorem (ref), 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.
When the composite loss function $\rho_C=\sum_{k=1}^K w_k\rho_k$ is used in the RAMP algorithm with tuning parameter $\alpha$ we denote the estimator at iteration number $t$ by $\widehat\beta_{\text{C},(t)}(\alpha)$. When the value of the tuning parameter is clear from the context, we also denote the RAMP estimator by $\widehat\beta_{\text{C},(t)}$.
For constructing the model averaging estimator we denote the separate estimators from the RAMP algorithm using regularity parameters $\alpha_k$, $k=1,\ldots,K$ by $\widehat\beta_{k,(t)}(\alpha_k)$ and the model-averaged estimator is denoted by $\widehat\beta_{\text{MA},(t)}(\bm{\alpha})=\sum_{k = 1}^K w_k \widehat\beta_{k, (t)}(\alpha_k)$ with $\bm{\alpha}=(\alpha_1,\ldots,\alpha_K)^\top$. When the value of the tuning parameters is clear from the context, we denote the model averaging RAMP estimator by $\widehat\beta_{\text{MA},(t)}$.
A generic estimator, without referring to a specific loss function or construction, is denoted by $\widehat{\beta}_{(t)}$, using tuning parameter $\alpha$; the subscript $(t)$ refers to the iteration number.
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 donoho2009message, this procedure for robust high-dimensional parameter estimation donoho2016high,bradic2016robustness adjusts the residuals to incorporate the valid score directly. While more details are given in Algorithm (ref), which is applied to the different loss functions $\rho_1,\ldots,\rho_K$ and to their weighted sum $\rho_C=\sum_{k=1}^K w_k\rho_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.
donoho2016high proposed to use the following proximal mapping operator to adjust the residuals. With $b>0$, $$\mbox{Prox}(z, b) = \arg \min_{x \in \mathbbm R} \{ b\rho(x) + \frac{1}{2}(x-z)^2 \}$$ which minimizes the square loss regularized by the non-differentiable loss, $\rho$. 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 $\rho$. The fixed point solution of the proximal operator coincides with the minimum of the loss function $\rho$. For more information, see ParikhBoyd2014.
We continue with the worked out example on quantile regression, see (ref). For $\ell = 0, \ldots, 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, $\sum_{i=a}^{b} x_i = 0$ if $b<a$. The proximal operator for the composite quantile case, see (ref), is
See Section (ref) for the derivation of the algorithm.
We now describe the three steps in more detail.
Step 1: Create adjusted residuals. \\ We use the estimates $\widehat\beta_{(t-1)}$ and $\widehat{\beta}_{(t)}$ from iteration steps $t-1$ and $t$ to compute the adjusted residuals
where the soft-thresholding function $\eta(x; \theta ) = \mbox{sign}(x)\max(|x|-\theta,0)$ and the score function $G$ is defined in (ref).
In Algorithm (ref), see Section (ref), we give details on how to set the soft-thresholding parameter $\theta$, which might change in each iteration, and we explain that a proper choice of $\theta$ as a function of the regularity constant $\lambda$ leads to an equivalence of the RAMP estimator and the regularized estimator.
The effective score function used in donoho2016high is
a subgradient is used in case of nondifferentiability. That is, for a value $x$ where $\rho$ is non-differentiable $$\partial\rho(x) = \{y: \rho(u)\ge \rho(x) + y(u-x), \forall u \}.$$ Throughout, we use $\partial_1$ as the notation for the partial derivative or partial subgradient of a function with respect to its first argument. Functions (e.g. $\widetilde G$) are applied componentwise to vectors.
For the example on composite quantile regression the subgradient of $\rho_{\rm{C}}$ is computed as,
where $h(\ell)$ is defined in (ref). The effective score function for composite quantile regression, see Section (ref), is
To incorporate the sparsity, bradic2016robustness, see also bayati2011dynamics, used the rescaled, min regularized effective score function,
where $\omega=E[\|B_0\|_0]$, see condition (ref) 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}\sum_{i = 1}^n\partial_1 G(z_{i,(t)}; b_{(t)})=1$. In the case of a non-differentiable loss function, bradic2016robustness proposed to solve $\widehat \nu(b_{(t)})=1$ with
Condition (ref) in the Appendix defines $\gamma_l, r_l$ and the differentiable functions $v_1$ and $v_2$ bradic2016robustness, $\widehat f_P$ is the estimated density of $\mbox{Prox}(z_{i,(t)};b_{(t)})$ for $i=1,\ldots,n$.
The derivation of the estimator $\widehat\nu(b_{(t)})$, see also Section (ref), relies on the limiting behaviour of the system, see Section (ref).
For the composite quantile loss, see (ref), we clearly see the dependence on the quantiles. The estimator of $\nu$ in (ref) uses $v_1(z) = 0$ and $v_2(z) = z - u_{\tau_\ell}$ $z\in [u_{\tau_\ell} + b h(\ell-1), u_{\tau_\ell} + b h(\ell)],\ \ell = 1, \ldots, K$, corresponding to the differentiable pieces in (ref). The step functions $v_3(z)=bh(\ell)$ when $z \in (u_{\tau_\ell} + bh(\ell), u_{\tau_{\ell} + 1} + bh(\ell)),\ \ell = 0, \ldots, K$. Solving for $b$ in the equation $\widehat\nu(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 $\widehat b_t$. For each $b$ in the grid, we use the empirical cumulative distribution, that is, $\widehat F_z(bh(K)) = n^{-1}\sum_{i=1}^n I\{ z_{i;(t)} \leq b h(K)\}$. A kernel density estimator of $f_z$ with the Gaussian kernel estimates defined as $\widehat f_z\{u_{\tau_k} +bh(k)\} = (nh)^{-1} \sum_{i =1}^n \phi\{ (z_{i;(t)} -u_{\tau_k} -bh(k))/h \}$ with $\phi$ being the standard normal density function. The solution $\widehat b_t$ is taken to be the average of the smallest $b$ in the grid that makes the righthand side of ((ref)) smaller than $\frac{s}{n}$ and the next value in the grid.
Step 3: Update the estimator of $\beta$. \\ Use the estimated $b_{(t)}$ from the previous step to update the estimate of $\beta$ to
The estimator $\widetilde\beta_{(t)}$, before applying the soft-thresholding function, is of interest too since it can be interpreted as a debiased estimator javanmard2014confidence, javanmard2014hypothesis, van2014asymptotically, javanmard2018debiasing; a thorough study of which, however, is beyond the current work.
\allowdisplaybreaks
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 $\bar\zeta_{(t)}^2$ which is critical for Algorithm (ref). We start by defining the empirical version as follows
This quantity is linked to the state evolution recursion which describes the limiting behaviour of large systems, see Theorem 2 in bayati2011dynamics and Lemma 1 in bradic2016robustness. It holds that, see Section (ref) for details,
where $\bar\zeta_{(t)}$ is the state evolution parameter for the large system, $Z$ is a random variable with standard normal distribution independent of everything else and $\bar \sigma(t)$ is defined in (ref).
Due to the symmetry of $Z$, the state evolution parameter is formally defined as
This definition explicitly features the extra Gaussian component $\bar\sigma_{(t)}Z$ in the limiting version, with variance
with $B_0$ defined in (ref). To connect the theoretical expression of $\bar\sigma_{(t)}^2$ to Algorithm (ref), we apply Eq.(3.6) in bayati2011dynamics, and Eqs.(7.10) and (7.19) in bradic2016robustness. This leads to
Note that (ref) features the debiased estimator from (ref).
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 bradic2016robustness, we get the immediate connection between the regularized estimators $\widehat{\beta}_k(\lambda_k)$ for $k=1,\ldots, K$ and the corresponding estimators obtained by applying the RAMP algorithm with a suitable choice of its regularity parameter $\alpha$. We explain this below. Since the regularized estimators $\widehat{\beta}_k(\lambda_k)$ for $k=1,\ldots,K$ are used for $\widehat{{\beta}}_{{\rm MA}}({\bm\lambda})$, (ref), 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 $\widehat{{\beta}}_{{\rm C}}(\lambda)$, (ref), is a special case of a model-averaged estimator with $K = 1$, weight equal to one, and loss function $\rho_{\rm C}=\sum_{k=1}^Kw_k\rho_k$.
Denote $(\bar\zeta^2, b)$ as the fixed point solution when the iteration number $t \to \infty$ of the following equations,
Note that (ref) is the state evolution recursion for the large system while in (ref) the first equality is the population version of the requirement in step 2 in Algorithm (ref) which states that $n^{-1}\sum_{i = 1}^n\partial_1 G(z_{i,(t)}; b_{(t)})=1$. The second equalities of both (ref) and (ref) follow by using the definition of $G$ in (ref), with $\widetilde G$ being defined in (ref).
Then, under assumptions (ref)--(ref) (see the Appendix), for the RAMP algorithm with $\theta=\alpha \overline\zeta$, where the tuning parameter $\alpha >0$ (which motivates the definition of $\theta_{(t)}=\alpha \bar\zeta_{(t)}$ in Algorithm (ref)), and for the $l_1$-optimization with
it follows by Theorem 2 of bradic2016robustness that
The convergence in (ref) explicitly connects the two composite estimators: one estimator uses an explicit $l_1$-regularization as in (ref), the other estimator is obtained via the RAMP algorithm. Similar results can be found in huang2020robust for a generalized AMP algorithm with non-negative convex loss function, and in bayati2011lasso 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 $\rho_k$, $k=1,\ldots,K$. When using explicit $l_1$-regularization as in (ref) with the regularization constants $\lambda_k$ matching as in (ref) the values $\theta_k=\alpha_k\bar\zeta$, for $k=1,\ldots,K$ that are used in the RAMP algorithm, again Theorem 2 of bradic2016robustness applies. It hence follows that $$ \lim_{t\to\infty}\lim_{p\to\infty}\frac{1}{p}\sum_{j=1}^p \{\widehat{\beta}_{{\rm MA},j}(\bm\lambda)- \widehat{\beta}_{{\rm MA},(t),j}(\bm\alpha)\}^2 = 0, \mbox{a.s.} $$
This section contains detailed theoretical developments for the composite as well as the model-averaged AMP estimators in high-dimensions.
We first define the asymptotic mean squared error as
Combining ((ref)) and (ref), we obtain
which corresponds to Eq.(3.4) in bradic2016robustness with $\widetilde\beta_{(t),j}$ the debiased estimator in ((ref)).
In Section (ref), we defined the empirical state evolution parameter $\bar\zeta_{\rm{emp}, (t)}^2$, and we described the connections between the empirical updates in Algorithm (ref) and the theoretical state evolution recursion, which connects to the theoretical expression of the AMSE. While Algorithm (ref) 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 $\widehat{\beta}_k$, $k=1,\ldots,K$, each obtained by Algorithm (ref). Consequently, the estimators $\widehat{\beta}_k$, $k=1,\ldots,K$ are correlated.
Lemma (ref) extends Theorem 2 in bayati2011dynamics and (3.16) in Lemma 1(b) in bayati2011dynamics to the almost sure convergence of the product for any two recursions among $K$ paralleled recursions. All proofs are contained in Appendix (ref).
Note that Algorithm (ref) belongs to the general recursion in bayati2011dynamics, the initial condition takes $q_{(0)} = - \beta$ and the $k$th estimator calculated by Algorithm (ref) takes $h_{k, (t+1)} = \beta - X^\top G(z_{k, (t)}; b_{k, (t)}) - \beta_{k, (t)}$.
We obtain at iteration $t$, for $k_1, k_2\in\{1,\ldots,K\}$,
where $Z_{k_1}$ and $Z_{k_2}$ are possibly dependent standard normal random variables.
Since the estimators $\widehat\beta_{k_r}$, $r=1,2$ use the same design matrix, a correlation between $Z_{k_1}$ and $Z_{k_2}$ exists (see Corollary (ref)) and contributes to the correlation between $\widehat \beta_{k_1}$ and $\widehat \beta_{k_2}$. Using Lemma (ref), we obtain the theoretical AMSE for the regularized model-averaged estimator.
Since the AMSE expression of the regularized model-averaged estimator is a quadratic function of the weight vector $w$, Corollary (ref) readily provides the lower bound of the AMSE as well as the weight vector reaching this lower bound. The $K$-vector $\mathbf 1_{K}$ consists of ones only.
The expression of the core matrix $\Sigma_{(t)}$, which is the limit matrix for $n, p \to \infty$, contains the random variable $B_0$ which satisfies assumption (ref) in the Appendix. Likewise, $\Sigma_{0, (t)}$ which is the limit matrix for fixed $p$ while $n \to \infty$, contains the true coefficient $\beta$ (see ((ref))). In practice, neither the true coefficient vector $\beta$ 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 $\Sigma_{0,(t)}$ relying only on sequences generated in Algorithm (ref).
Before deriving the estimator of the AMSE for the model-averaged estimator, we first define $\bar\zeta_{\textrm{emp}, (k_1, k_2), (t)}$ which is an estimator of the parameter $\bar\zeta_{(k_1, k_2), (t)}$, a quantity similar to the state evolution parameter $\bar\zeta_{k, (t)}^2$, which records the covariance between the unbiased sequences $\widetilde \beta_{k_1, (t)}$ and $\widetilde \beta_{k_2, (t)}$ generated in ((ref)) in Algorithm (ref) when $p \to \infty$. Since model-averaged estimators combine estimators constructed from the same data into one weighted average, the correlation between $\widehat\beta_{k_1}$ and $\widehat\beta_{k_2}$ is needed to understand the AMSE of the model-averaged estimator.
Notice that the unbiasedness of the sequence $\widetilde\beta_{k, (t)}$ follows from the argument that $\widetilde\beta_{k, j, (t)}$ converges weakly to $B_0 + \bar\zeta_{k, (t)} Z_k$ when $p \to \infty$, while assigning $1/p$ point mass to each entry of the vector. Then, $\widetilde \beta_{k, j, (t)} | (B_0 = \beta_j) \sim N(\beta_j, \bar\zeta^2_{k, (t)})$ for large $p$, indicating that $\widetilde \beta_{k, j, (t)}$ centers at $\beta_j$ ensuring the unbiasedness. Moreover, the vector $\widetilde \beta_{k, (t)}$ has Gaussian distribution. By applying the soft-thresholding function $\eta$ on $\widetilde \beta_{k, j, (t)}$ in Lemma (ref), we avoid the usage of the true coefficient vector $\beta$ in $\Sigma_{0, (t)}$ resulting in a Stein-type risk estimator requiring only observables from Algorithm (ref). A Gaussianity argument has also been used in bayati2011lasso, bayati2013estimating, mousavi2013parameterless, mousavi2018consistent to derive a similar Stein-type risk estimator for the Lasso. Details can be found in Section (ref). The bias of the estimator $\widehat\beta_{k, (t)}$ is introduced in Algorithm (ref) by applying the soft-thresholding function componentwise to the unbiased sequence $\widetilde\beta_{k, (t)}$.
Corollary (ref) indicates both the existence and a feasible estimation of the covariance between $Z_{k_1}$ and $Z_{k_2}$. As an estimator for $$\bar\zeta_{(k_1, k_2), (t)} = \bar\zeta_{k_1, (t)}\bar\zeta_{k_2, (t)} \mathrm{Cov}(Z_{k_1}, Z_{k_2})$$ we define
We now state an unbiased estimator for the matrix $\Sigma_{0, (t)}$, and a consistent estimator for the matrix $\Sigma_{(t)}$ upon convergence of Algorithm (ref).
This new estimator can be compared to the estimator used in bayati2013estimating and mousavi2018consistent for the case of a single estimator ($K=1$). The proof of Theorem (ref), see Section (ref) uses Stein's lemma (see Lemma (ref)) to estimate the covariances that appear in the matrix $\Sigma_{0,(t)}$. The soft-thresholding function $\eta(\cdot; \theta)$ that appears in the estimator $\widehat\Sigma_{0,(t)}$ links the estimator $\widehat{\beta}_k$ to the estimator $\widetilde\beta_k$. The proof also uses the joint asymptotic normality of the $j$th components of the vectors $\widetilde\beta_{k_1}$ and $\widetilde\beta_{k_2}$. The obtained estimator for $\Sigma_{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 $\widehat\Sigma_{0,(t)}$ at the final iteration in Theorem (ref). In combination with the sum-to-one constrained weights this gives the estimated weights that minimize the estimated AMSE for the model-averaged estimator $$\widehat{w}_{\rm MA} = \big(\widehat\Sigma_{(t)}\big)^{-1} \mathbf 1_{K} \big(\mathbf 1_{K}^\top (\widehat\Sigma_{(t)})^{-1} \mathbf 1_{K} \big)^{-1}.$$ When additional constraints such as positivity are needed, the optimal weights no longer have an explicit formula, but they are straightforward to compute, see (ref).
The AMSE of a composite estimator can be obtained from Theorem (ref) as a special case, treating the composite loss function as a single loss function with weight one, thus $\rho_{\rm C}=\sum_{k=1}^Kw_k\rho_k$ as in (ref). At iteration $t$, $$ \Sigma_{(t)} = E[\{\eta(B_0 + \bar\zeta_{(t-1)} Z; \theta_{(t-1)}) - B_0 \}]^2, \mbox{ and } \Sigma_{0, (t)} = p^{-1}\sum_{j = 1}^p (\widehat\beta_{(t), j} - \beta_j)^2. $$ The matrices $\Sigma_{(t)}, \Sigma_{0, (t)}$ are now real numbers and coincide with the AMSE of the estimator in ((ref)). We obtain the corresponding estimator for the AMSE
For the single loss function, $\rho_{\rm C}$, the estimator of AMSE in (ref) can be compared to the Stein-type estimator that has been obtained in Definition 2 in bayati2013estimating for the AMP algorithm using the least squares loss, which is a particular case of Algorithm (ref).
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 (ref), the composite estimator and its AMSE depend on the weights in a highly nonlinear fashion; e.g., observe that the soft-thresholding function in (ref) depends on $w$.
Therefore, optimization of the estimated AMSE with respect to the weights proceeds numerically; $$ w_{\rm C, 1} = \arg\min_{w} \widehat{\mbox{AMSE}}_{{\rm C}}(w).$$ See Section (ref) for more details.
donoho2016high and el2013robust showed that the asymptotic variance of the M-estimators in the case where $p, n \to \infty$ and $n/p \to \delta \in [1, \infty)$ contains an extra Gaussian component. Recently, lei2018asymptotics 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 $\delta \in (0, 1)$, it was shown that the sequence $\widetilde \beta_{(t)}$ in ((ref)) follows for the Lasso estimator bayati2013estimating 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 \to \delta \geq 1 $, following donoho2016high.
Under the dense linear model with $n \geq p$, the soft-thresholding function $\eta(\cdot; \theta)$ is replaced by the identity function and the ratio $\omega = E[\| B_0\|_0]=1$. Consequently, Algorithm (ref) is adjusted to estimate $$\widehat\beta_k = \arg\min_{\beta \in \mathbbm R^p}\Big\{ \sum_{i = 1}^n \rho_k(Y_i - X_i\beta) \Big\},$$ where $\beta$ is dense. It is trivial to show that Algorithm (ref) still belongs to the general recursion in bayati2011dynamics. For a single estimator at iteration $t$ denoted as $\widehat\beta_{k, (t)}$, the two state evolution parameters $\bar\zeta_{k, (t)}^2$ and $\bar\sigma_{k, (t)}^2$ coincide and Theorem 4.1 in donoho2016high holds.
Recall that the componentwise empirical distribution of $\widehat\beta_k(p)$, when $p\to \infty$, converges weakly to $B_0 + \bar\zeta_{k} Z_k$ following bayati2011dynamics and donoho2016high. Then for large $p$, while the iteration $t \to \infty$, $\widehat\beta_k(p) \sim N(\beta, \bar\zeta_{k}^2 I_p)$ donoho2016high, mousavi2013parameterless with $I_p$ the $p\times p$ identity matrix. The $(k_1, k_2)$th component of the empirical variance matrix is denoted by $\big(\Sigma_{\rm emp}(p)\big)_{(k_1, k_2)} = p^{-1} \sum_{j = 1}^p (\widehat\beta_{k_1, j} - \beta_j)(\widehat\beta_{k_2, j} - \beta_j)$, which is unbiasedly estimated by $$ \big(\widehat\Sigma_{\textrm{emp}}(p)\big)_{(k_1, k_2)} = \sum_{j = 1}^p (\widehat\beta_{k_1, j} - \frac{1}{p} \sum_{j = 1}^p \widehat\beta_{k_1, j})(\widehat \beta_{k_2, j} - \frac{1}{p}\sum_{j = 1}^p \widehat\beta_{k_2, j}) / (p - 1). $$ Note that this estimator coincides with (ref) for the special case that $n \geq p$ and the soft-thresholding function is replaced by the identity function.
The estimation of the quantile $u_{\tau_k}= F_\varepsilon^{-1}(\tau_k)$ follows a two-step procedure.
The regularized model-averaged estimator is obtained by averaging over $K$ paralleled estimators. See Algorithm (ref) for the pseudo-code, of which the core is Algorithm (ref); there the effective score function $G$ is that of a single quantile loss function with $K=1$, see also Example 2 in bradic2016robustness. 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.
The tuning parameter $\alpha$ of Algorithm (ref) controls the sparsity of the estimators and requires a tuning procedure to choose it in practice. In Section (ref), we consider the one dimensional Golden-section search algorithm kiefer1953sequential for tuning the value $\alpha$ in the range $[\alpha_{\rm min}, \alpha_{\rm max}]$ that minimize the estimated MSE of $\widehat\beta$ using the estimator derived in Section (ref). The upper bound $\alpha_{\rm max}$ is chosen to be 2.3 for the simulations and data analysis. The lower bound $\alpha_{\rm min}$ in the data analysis follows the lower bound in Proposition 9.2 in eldar2012compressed and is chosen to be the unique non-negative solution to the equation $(1 + \alpha^2)\Phi(-\alpha) - \alpha \phi(\alpha) = \delta/2$, where $\phi(x)$ and $\Phi(x)$ denote the p.d.f and c.d.f of the standard normal distribution respectively. In the simulation study, the lower bound $\alpha_{\rm 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.
To obtain the regularized model-averaged quantile estimations with the AMSE-type weight derived in Corollary (ref), we follow the following procedure:
It is worth mentioning that $\widehat w_{\rm MA, 1}$ is a constrained version of $w_{\rm MA}$ attaining the lower bound of the AMSE in Corollary (ref). $\widehat w_{\rm MA, 1}$ focuses on approximating the lower bound of the AMSE of the sparse coefficient vector $\beta$ without assuming that the nonzero entries are selected perfectly; whereas another type of weight choice derived in bradic2011penalized, BloznelisClaeskensZhou2019 aims at the lower bound of the variance of the nonzero part of $\beta$ by imposing the perfect selection assumption. A numerical comparison of these two types of weight choices is presented in Section (ref).
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 $\widehat\beta_{\rm C}(w_{\rm 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 $\widehat\beta_{\rm C}(w_{\rm C})$ is updated in each step by the one having the lowest estimated AMSE, i.e., $$w_{\rm C, 1} = \arg\min_{w_{\rm cand}} \widehat{\mbox{AMSE}}(\widehat\beta_{\rm C}(\alpha_{\rm opt}; w_{\rm cand}); \beta). $$
A more detailed search procedure is as follows.
The pseudocode of the search procedure is stated in Algorithm (ref).
In this section, we consider the following setup under the high-dimensional linear model setting.
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_{\rm MA, 1}$ for the model-averaged quantile estimator and $w_{\rm 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_{\rm MA, 2}$ BloznelisClaeskensZhou2019 and $w_{\rm C, 2}$ bradic2011penalized where, with the $(k_1, k_2)$th component of $A$ equal to $A_{k_1, k_2} = \min(\tau_{k_1},\tau_{k_2})\{1 - \max(\tau_{k_1},\tau_{k_2})\}$, $A_\varepsilon= \textrm{diag}( f_\varepsilon(u_{\tau_1}), \ldots, f_\varepsilon(u_{\tau_K}))$, and $a_\varepsilon = ( f_\varepsilon(u_{\tau_1}), \ldots, f_\varepsilon(u_{\tau_K}))^\top$
and
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 “non-zero 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 $$ \mbox{MSE}(\widehat\beta_{\rm part}) = \sum_{j_{\rm part} = 1}^{p_{\rm part}} (\widehat\beta_{j_{\rm part}} - \beta_{j_{\rm part}})^2 / p_{\rm part} $$ for the appropriate part of the full vectors. Results for the regularized model-averaged quantile estimator with different weights are presented in Table (ref). We observe that the model-averaged quantile estimator using the weight in ((ref)) has lower MSEs for estimating the non-zero part of $\beta$ and for the full vector $\beta$, 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 $\beta$. The Lasso estimator is considered as the baseline comparison, which from Table (ref) 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 (ref) summarizing the averaged true positive (TP) and true negative (TN) recovery rates which are defined as $$\mbox{TP (TN)} = \frac{\mbox{number of correctly identified as non-zeros (zeros)}}{\mbox{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 ((ref)) 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_{\rm 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_{\mathcal V}$ to be five and only randomly select 4 points in the neighborhood $\mathcal{V}_{w_{\rm cand}}$; the tuning parameter $\alpha$ of the soft-thresholding function is tuned once for the regularized composite quantile estimator with the weight $w_{\rm C, 2}$, then fixed after that.
Table (ref) summarizes the empirical MSEs of the regularized composite quantile estimator with different weights. Since the tuning parameter, $\alpha$ is selected for $w_{\rm 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_{\rm C, 2}$ leads to lower MSEs in most cases. However, it is worth noticing that using equal weights, while $\alpha$ 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 (ref). The regularized composite estimator with locally optimized $w_{\rm C, 1}$ consistently has the highest TP rate, and second-highest TN rate among all competitors, except the TN rate for $t_3$ distributed errors and TP rate for $0.5N(0,1) + 0.5 N(5, 9)$ distributed errors. At the same time, the non-zero parts of $\beta$ are generated from Dirac distribution at -1 and 1.
Tables (ref) and (ref) 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 (ref), column 7 and 12, and Table (ref), column 7, with the single regularized quantile estimator at the median $\tau = 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 $\beta$ in settings where errors are generated from $N(0, 1)$ and $0.5N(0,1) + 0.5 N(5, 9)$ distribution and the true non-zero subvector of $\beta$ 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 $\beta$ is generated from $N(0, 1)$.
The percentage of converged cases for the model-averaged and composite estimator, while setting the tolerance $\varepsilon_{{\rm tol}}$ to be $10^{-6}$ for different error distributions, are included in Table (ref), where we define a estimator to have converged when the needed number of iterations was less than 50.
Condition (ref) restricts Algorithm (ref) to a special design matrix that does not allow correlations between the $X_{\cdot j}$'s. However, since such correlation might be present in reality, it is of interest to see if Algorithm (ref) is still numerically robust while Condition (ref) is relaxed in practice. We consider a similar simulation setup as used before with $p=500$, the sample size $n = 250$, and $\delta = 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, \Sigma_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 $\Sigma_{X}$ of which its $(i, j)$th component $(\Sigma_{X})_{i,j} = \sigma_{X}^{|i - j|}, i, j = 1, \ldots, p$. We consider $\sigma_{X} = 0, 0.1, 0.3$. \\ 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 (ref) records the performance of Algorithm (ref) with tolerance $\varepsilon_{\rm 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 (ref) that parameter estimation using Algorithm (ref) remains accurate and stable when weak correlations such as with $\sigma_X = 0.1$ exist between the $X_{\cdot j}$'s; the accuracy drops when we further increase the correlations as with $\sigma_X= 0.3$; it is worth mentioning that the convergence percentages decrease when the correlation increases. Further research concerning correlated data is worth considering.
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.
Originated from the compressed sensing problem, the sparse linear model $Y = X\beta + \bm{\varepsilon}$ describes the image or signal compression. The $s$-sparse $p$-dimensional input signal $\beta$ is first compressed by a known sensing matrix $X \in \mathbbm R^{n\times p}$ with $n < p$; the compressed signal vector $X \beta \in \mathbbm R^{n}$ can be corrupted by the noise $\varepsilon$ with $\varepsilon_i$'s i.i.d. via transmission. Notice that the $p$-dimensional input signal vector $\beta$ is assumed to be $s$-sparse which is usually unsatisfied by signals expressed in the standard basis. To obtain the sparse representation of $\beta$ in practice, an intermediate stage of expressing the natural non-sparse vector $\beta^*$ in a proper orthonormal basis $ \Psi^* = (\psi_1^*, \ldots, \psi_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.
To imitate the compressed sensing process, we process the audio wave signal vector as follows:
In practice, the artificial vector $Y$ and the sensing matrix $X$ are observed. The accurate recovery of the original wavelet coefficient vector $\beta$ 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 $\bm\varepsilon$ under various seed numbers. However, we only present the reconstructions under one seed for each setting in Section (ref) due to limited space.
To reconstruct the signal vector $\beta$ expressed in the wavelet basis from the sensing matrix $X$ and the observed compressed signal vector $Y$ corrupted by potentially non-Gaussian distributed error $\bm \varepsilon$, 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 $\varepsilon_{\rm tol}=10^{-8}$. The Lasso estimator is considered as the baseline comparison. Notice that the regularized estimates $\widehat\beta_{\rm MA}$ and $\widehat\beta_{\rm C}$ after reconstruction are the representations in the wavelet domain. To compare the accuracy of the reconstruction, we perform a back-transform on the estimates and obtain the corresponding signal vectors $\widehat{\beta}_{\rm MA}^*$ and $\widehat{\beta}_{\rm C}^*$ with representations in the natural basis.
Example reconstructions of the audio signal for $K = 3$ using the regularized model-averaged estimator equipped with different weights, with the baseline recovery from the Lasso represented in the natural basis are presented in Figure (ref) for the mixture of normals distributed error, and in Figure (ref) 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 model-averaged estimators using $\widehat{w}_{{\rm MA, 1}}$ and equal weights outperform the counterpart with $\widehat{w}_{{\rm MA, 2}}$ for $0.5 N(0,1) + 0.5 N(5, 9)$ distributed errors; recovery differences for the weak signals of 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_{\rm 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.
BatesGranger1969 provide an alternative weight choice for the model-averaged estimator obtained by considering only the variances of $\widehat\beta_k$'s and ignoring the covariances. This leads to
where $\mbox{diag}(\widehat\Sigma_{0,(t)})$ denotes the diagonal matrix obtained from $\widehat\Sigma_{0,(t)}$ which keeps the diagonal and has zeros in all off-diagonal entries. Figure (ref) contains the recovery of the audio signal using the model-averaged estimator using this weight.
For the composite quantile estimator $\widehat\beta_{\rm C}$, we performed the same weight searching method as for the simulation study. This is, $S_{\mathcal V}=5$ and randomly select 4 candidate weights in the neighbourhood of the previous value. We select the tuning parameter $\alpha$ once for the starting weight $\widehat w_{\rm 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 (ref) where the MAPE is defined as
Table (ref) 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 $\widehat w_{\rm MA, 1}$ 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 $\widehat w_{\rm 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 $\alpha$ tuned only once for the starting weight.
Additionally, we present the estimated weights for both regularized model-averaged and composite estimators in Table (ref). An interesting observation is made by comparing the estimated weights $\widehat w_{\rm MA, 1}$ and $\widehat w_{\rm MA, 2}$ for the mixed normal distributed error. The weight $\widehat w_{\rm 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, $\widehat w_{\rm 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.
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 $\lambda$ (for composite estimation) and $\lambda_1,\ldots,\lambda_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 $\beta$ and $w$.
To avoid overly complicated mathematical expressions, we followed earlier literature in the use of a design matrix where $X_{ij}\sim N(0, 1/n)$. Other applications might require studying, for example, fixed designs, which are beyond the scope of the current paper.