Robust High-Dimensional Regression with Coefficient Thresholding and its Application to Imaging Data Analysis

In the high-dimensional linear regression ${\bf{y}}={\bf{X}}\boldsymbol{\beta}^{\star}+\boldsymbol{\varepsilon},$ ${\bf{y}}=(y_{1},y_{2},\dots,y_{n})^{T}$ is an $n$-dimensional response vector, ${\bf{X}}=({\bf{x}}_{1},\dots,{\bf{x}}_{p})$ is an $n\times p$ deterministic design matrix, $\boldsymbol{\beta}^{\star}$ is a $p$-dimensional true regression coefficient vector and $\boldsymbol{\varepsilon}$ is an error vector with mean zero. We wish to recover the true sparse vector $\boldsymbol{\beta}^{\star}$, whose support $S(\boldsymbol{\beta}^{\star})=\{j:\beta_{j}^{\star}\neq 0\}$ satisfies that its cardinality $|S(\boldsymbol{\beta}^{\star})|\ll p$. We use $S$ as the shorthand for the support set $S(\boldsymbol{\beta}^{\star})$.

The irrepresentable condition is known to be sufficient and almost necessary for the variable selection consistency of the $\ell_{1}$ penalization (Zhao & Yu, 2006; Zou, 2006; Wainwright, 2009). Specifically, the design matrix ${\bf{X}}$ should satisfy that

for some incoherence constant $\gamma\in[0,1)$, This condition requires that the variables in the true support are weakly correlated with other variables that are not in the true support.

There are several versions of relaxation of the irrepresentable condition in the current literature. Zou (2006) used the adaptive weights $\hat{w}_{j}$’s in the $\ell_{1}$ penalization and relaxed the irrepresentable condition as:

where $\circ$ denotes the Hadamard (or componentwise) product of two vectors. The folded concave penalization  (Fan & Li, 2001; Zhang, 2010) also relaxed the irrepresentable condition. Let $p_{\lambda}(\cdot)$ be the folded concave penalty function. Given the local linear approximation and the current solution $\tilde{\boldsymbol{\beta}}$, the folded concave penalization relaxed this condition as:

where $\dot{p}_{\lambda}(\cdot)$ is the subgradient of the function $p_{\lambda}(\cdot)$. Both the adaptive $\ell_{1}$ penalization and folded concave penalization utilized the differential penalty functions to relax the restrictive irrepresentable condition. The procedure of adaptive lasso and solving folded concave penalized problem using local linear approximation (LLA) both require a good initialization. Their dependence on the initial solution is inevitably affected by the irrepresentable condition. As a promising alternative, we consider the assignment of adaptive weights on the design matrix ${\bf{X}}$, which down-weights the unimportant variables. Consider the ideal adaptive weight function $g(\boldsymbol{\beta}):=(g(\beta_{1}),\dots,g(\beta_{p}))$, where $g(\beta_{j})$ goes to $1$ when $j\in S$ and goes to $0$, otherwise. We propose the following nonconvex estimation procedure:

Let ${\bf{z}}_{j}=g(\boldsymbol{\beta})\circ{\bf{x}}_{j}$. The irrepresentable condition is now relaxed in this paper as follows:

Given the ideal adaptive weight function $g(\boldsymbol{\beta})$, the proposed approach further relaxes the irrepresentable condition in comparison to those considered by existing methods. Unlike existing methods, we propose an integrated nonconvex estimation procedure without depending on initial solutions.

The weight function plays an important role in relaxing the irrepresentable condition. If we knew the oracle threshold $\eta^{\star}=\underset{j\in S}{\min}\{|\beta_{j}^{\star}|\}$, the hard thresholding function $I\{|\beta_{j}|\geq\eta^{\star}\}\}$ would be the ideal choice for the weight function $g(\boldsymbol{\beta})$. However, the discontinuity of $I\{|u|\geq\eta\}$ is challenging for associated optimization. To overcome, we consider a smooth approximation of $I\{|u|\geq\eta\}$ given by

where $h_{\tau}(w)=1/2+\arctan(w/\tau)/\pi$. Since $h_{\tau}(w)\to I\{w\geq 0\}$ as $\tau\to 0$ when $w\neq 0$, we have $g_{\tau,\eta}(u)\to I\{|u|\geq\eta\}$ as $\tau\to 0$ when $u\neq 0$. Figure 2 illustrates the smooth approximation of $\beta\circ g_{\tau,\eta}(\beta)$ to $\beta\circ I\{|\beta|\geq\eta\}$ when $\tau$ is small (e.g., $\tau=0.001,0.005,0.01).$

Given the above smooth thresholding function $g_{\tau,\eta}(\cdot)$, we propose the robust regression with coefficient thresholding:

where following Mei et al. (2018) and Loh & Wainwright (2017) we assume that the regression coefficients $\boldsymbol{\beta}$ are bounded in the Euclidean ball ${\bf B}^{p}(r)\equiv\{\boldsymbol{\beta}\in\mathbb{R}^{p}:\ \|\boldsymbol{\beta}\|_{2}\leq r\}$. Here, to tackle possible outliers, we choose $L(\cdot)$ to be a differentiable and possibly nonconvex robust loss function such as the pseudo Huber loss (Charbonnier et al., 1997) or Tukey’s biweight loss. For the ease of presentation, throughout this paper, we use the pseudo Huber loss of the following form:

Note that $L(a)$ provides a smooth approximation of the Huber loss (Huber, 1964) and bridges over the $\ell_{1}$ loss and the $\ell_{2}$ loss. Specifically, $L(a)$ is approximately an $\ell_{2}$ loss $a^{2}/2$ when $a$ is small and an $\ell_{1}$ loss with slope $\omega$ when $a$ is large. It is worth pointing out that the utility of the pseudo Huber loss essentially adds a weight on each observation $x_{i}$ by $w(x_{i})$ in the following form:

where $\hat{y_{i}}$ is a fitted value. Hence, outliers are down-weighted to alleviate potential estimation bias.

To deal with a large number of predictors, we propose the following penalized estimation of high-dimensional robust regression with coefficient thresholding:

where $p_{\lambda}(|\boldsymbol{\beta}|)$ is a chosen penalty function such as the $\ell_{1}$ penalization or folded concave penalization. In the scalar-on-image analysis, often we need to incorporate some known group information into $p_{\lambda}(|\boldsymbol{\beta}|)$. Suppose that the coefficient vector $\boldsymbol{\beta}$ is divided into $B$ separate subvectors $\boldsymbol{\beta}_{1},.\cdots,\boldsymbol{\beta}_{B}$. We can use the group penalty $\lambda\sum_{b=1}^{B}\|\boldsymbol{\beta}_{b}\|_{2}$ (Yuan & Lin, 2006) or sparse group penalty $\lambda_{1}\|\boldsymbol{\beta}\|_{1}+\lambda_{2}\sum_{b=1}^{B}\|\boldsymbol{\beta}_{b}\|_{2}$ (Simon et al., 2013) in the penalized estimation (3). With $p_{\lambda}(|\boldsymbol{\beta}|)=\lambda\sum_{b=1}^{B}\|\boldsymbol{\beta}_{b}\|_{2}$, the sparse robust regression with coefficient thresholding and group selection can be written as follows:

which includes the $\ell_{1}$ penalization as a special example with $B=p$. The penalized robust regression with coefficient thresholding (3) and (4) can be efficiently solved by a customized composite gradient descent algorithm with provable guarantees, and the details will be presented in Subsection 3.3.

It is known that $\ell_{0}$-penalization enjoys the oracle risk inequality (Barron et al., 1999). Fan & Lv (2013) proved that the global solution of $\ell_{0}$-penalization, denoted by $\boldsymbol{\hat{\beta}}^{\mathcal{O}}=(\hat{\beta}^{o}_{1},\cdots,\hat{\beta}^{o}_{p})^{T}$, satisfies the hard thresholding property that $\hat{\beta}^{o}_{j}$ is either $0$ or larger than a positive threshold whose value depends on ${\bf{X}}$ and ${\bf{y}}$ and tuning parameter $\lambda$. Specifically, Fan & Lv (2013) studied the oracle properties of regularization methods over the thresholded parameter space $\mathcal{B_{\eta}}$ defined as

Given the thresholded parameter space $\mathcal{B_{\eta}}$, we introduce the high-dimensional regression with coefficient hard-thresholding as follows:

This thresholded parameter space guarantees that the hard thresholding property is satisfied with proper $\ell_{2}$ and $\ell_{0}$ norm upper bounds. Suppose that global solutions to both (3) and (8) exist. It is intuitive that if $\tau$ is very small, an optimal solution to (3) should be “close” to that of (8). Apparently, $g_{\tau,\eta}(\cdot)$ converges to $I(|\cdot|\geq\eta)$ pointwisely almost everywhere as $\tau\mapsto 0^{+}$. However it is known that almost surely pointwise convergence does not guarantee the convergence of minimizers (e.g., (Rockafellar & Wets, 2009, Section 7.A)). Here we provide the following proposition to show that the global solution of the proposed method (3) converges to the minimizer of coefficient hard-thresholding in (8) with additional assumptions.

Given Proposition 1, we can prove that the total number of falsely discovered signs of the cluster point $\boldsymbol{\hat{\beta}}$, i.e., the cardinality $\mathrm{card}(\{j:\mathrm{sgn}(\boldsymbol{\hat{\beta}})\neq\mathrm{sgn}(\boldsymbol{\beta}^{\star})\})$, converges to zero for the folded concave penalization under the squared loss and sub-Gaussian tail conditions. Recall that the folded concave penalty function $p_{\lambda}(t)$ satisfies that (i) it is increasing and concave in $t\in[0,\infty)$, (ii) $p_{\lambda}(0)=0$, and (iii) it is differentiable with $\dot{p}_{\lambda}(0_{+})=a\lambda$ for some $a>0$. Define $\kappa_{c}$ as the smallest possible positive integer such that there exists an $n\times\kappa_{c}$ submatrix of $n^{-1/2}{\bf{X}}$ having a singular value less than a given positive constant $c$. $\kappa_{c}$ is named robust spark in (Fan & Lv, 2013) to ensure model identifiability. Note that $\mathrm{card}(\{j:\mathrm{sgn}(\boldsymbol{\hat{\beta}})\neq\mathrm{sgn}(\boldsymbol{\beta}^{\star})\})$ is an upper bound on the total number of false positives and false negatives. The following proposition implies the variable selection consistency of the cluster point $\boldsymbol{\hat{\beta}}$.

Propositions 1 and 2 indicate that the cluster point $\boldsymbol{\hat{\beta}}$ enjoys the similar properties as those of the global solution given by the $\ell_{0}$-penalized methods. The proposed method mimics the $\ell_{0}$-penalization to remove irrelevant predictors and improve the estimation accuracy as $\tau\to 0$.

In the sequel, we will show that with high probability our proposed method has a unique global minimizer (See Theorem 2(b)). In this aspect, our proposed regression with coefficient thresholding estimator provides a practical way to approximate the global minimizer over the thresholded space.

After presenting technical conditions, we first provide the landscape analysis and establish the uniform convergence from empirical risk to population risk in Subsection 3.2.1 and then prove the oracle inequalities for the unique minimizer of the proposed method in Subsection 3.2.2.

Now we introduce the notation to simplify our analysis. Let $g(u)$ be a shorthand of $g_{\tau}(u,\eta)$, and let $G(\boldsymbol{\beta})=(\beta_{1}g(\beta_{1}),\beta_{2}g(\beta_{2}),\cdots,\beta_{p}g(\beta_{p}))^{T}$ on $\mathbb{R}^{p}$. Given fixed $\tau$, we claim that $G(\boldsymbol{\beta})$ is third continuously differentiable on its domain. After the direct calculations given in Lemma 1 of the supplementary file, we have the explicit upper and lower bounds for the derivatives of $G(\boldsymbol{\beta})$, i.e.,

where $D_{G}(\boldsymbol{\beta}):=\nabla G(\boldsymbol{\beta})$, $D^{2}_{G}(\boldsymbol{\beta}):=\nabla D_{G}(\boldsymbol{\beta})$, and $D^{3}_{G}(\boldsymbol{\beta}):=\nabla D^{2}_{G}(\boldsymbol{\beta})$. Here, $\underline{k_{0}},\overline{k_{0}},\overline{m_{0}}$ and $\overline{s_{0}}$ are uniform constants independent of $\boldsymbol{\beta}$. And $A\preccurlyeq B$ represents that $B-A$ is semi positive definite.

For the group lasso penalty, denote by $\boldsymbol{\beta}^{a}=(\boldsymbol{\beta}_{1}^{a^{T}},...,\boldsymbol{\beta}_{B_{1}}^{a^{T}})^{T}$ the non-zero groups and by $\boldsymbol{\beta}^{c}=(\boldsymbol{\beta}_{1}^{c^{T}},...,\boldsymbol{\beta}_{B_{2}}^{c^{T}})^{T}$ the zero groups, respectively. Following the structure of classical group lasso, let $d_{i}^{a}$ and $d_{i}^{c}$ be the corresponding length of vector $\boldsymbol{\beta}_{i}^{a}$ and $\boldsymbol{\beta}_{i}^{c}$, respectively. Let $d^{a}_{max}=\underset{i\in\{1,\dots,B_{1}\}}{\max}\ d_{i}^{a}$, $d^{c}_{max}=\underset{i\in\{1,\dots,B_{2}\}}{\max}\ d_{i}^{c}$, and $d_{\max}=d^{a}_{\max}\vee d^{c}_{\max}$, where $d_{\max}$ is finite, not diverging with $n$ and $p$.

We make the following assumptions on the distribution of predictor vector ${\bf{x}}$, the true parameter vector $\boldsymbol{\beta}^{\star}$ and the distribution of random error $\boldsymbol{\varepsilon}$.

Assumption 1(a) presents the technical conditions on the design matrix. The same conditions were considered in Mei et al. (2018) for binary linear classification and robust regression. The sub-Gaussian assumption is a commonly used mild condition in high-dimensional regression. Assumption 1(b) imposes the sparsity on the true parameter vector $\boldsymbol{\beta}^{\star}$. Assumption 1(c) permits the possibly heavy tails of the error. We allow the size of the true support set to diverge at rate $o(n)$. Given the sparsity, it is reasonable to limit our theoretical analysis in the Euclidean ball $\Theta_{n,p}={\bf B}^{p}(r)\equiv\{\boldsymbol{\beta}\in\mathbb{R}^{p},\|\boldsymbol{\beta}\|_{2}\leq r\}$, which can avoid unnecessary technical complications. Also, $\eta$ does not need to be a constant bounded away from $0$ in the asymptotic analysis. Since $G(\boldsymbol{\beta})$ is continuous and injective, there exists a unique $\boldsymbol{\tilde{\beta}}^{\star}$, such that $\|\boldsymbol{\tilde{\beta}}^{\star}\|_{2}\leq r/3$ and $G(\boldsymbol{\tilde{\beta}^{\star}})=\boldsymbol{\beta}^{\star}$, i.e. $\boldsymbol{\beta}^{\star}$ is the thresholded version of $\tilde{\boldsymbol{\beta}}^{\star}$. In the sequel, we will study the statistical convergence to the surrogate $\tilde{\boldsymbol{\beta}}^{\star}$, which shares the same non-zero support with $\boldsymbol{\beta}^{\star}$. In Theorem 2, we show that $\|\boldsymbol{\hat{\beta}}-\boldsymbol{\tilde{\beta}}^{\star}\|_{2}=O_{p}(\sqrt{\frac{s_{0}\log p}{n}})$. Then $\|G(\boldsymbol{\hat{\beta}})-G(\boldsymbol{\tilde{\beta}}^{\star})\|_{2}=\|G(\boldsymbol{\hat{\beta}})-\boldsymbol{\beta}^{\star}\|_{2}=O_{p}(\sqrt{\frac{s_{0}\log p}{n}})$ given $G(\cdot)$ is Lipschitz with fixed $\eta$. In this way, we can use $G(\boldsymbol{\hat{\beta}})$ as an estimator for $\boldsymbol{\beta}^{\star}$. Assumption 1(d) allows random error with heavier tails than the standard Gaussian distribution, and it suits for many applications with outliers. For example, in our simulation, our noise is a mixture distribution with a small variance Gaussian distribution and a large variance Gaussian distribution. Define $\varphi(\cdot)=L^{\prime}(\cdot)$ and $h(\cdot)=\mathbb{E}_{\boldsymbol{\varepsilon}}[\varphi(\cdot+\boldsymbol{\varepsilon})]$. Assumption 1(d) can be relaxed as Assumption 1(d’) $h(0)=\mathbb{E}_{\boldsymbol{\varepsilon}}[\varphi(\boldsymbol{\varepsilon})]\geq 0$, which holds for the right skewed error or discrete error. Our theoretical results given in this paper still hold under Assumption 1(a)–(c) and 1(d’). See Remark 6 for more details.