Robust Regression via Online Feature Selection under Adversarial Data Corruption

TL;DR

This paper introduces RoOFS, a robust online feature selection algorithm that effectively learns reliable regression coefficients from streaming, partially accessible, and corrupted data, addressing scalability and adversarial challenges.

Contribution

The paper presents a novel algorithm that simultaneously handles corrupted data estimation, online feature selection, and scalability in streaming environments.

Findings

Outperforms existing methods in feature recovery accuracy

Provides theoretical error bounds for the algorithm

Demonstrates effectiveness on synthetic and real-world datasets

Abstract

The presence of data corruption in user-generated streaming data, such as social media, motivates a new fundamental problem that learns reliable regression coefficient when features are not accessible entirely at one time. Until now, several important challenges still cannot be handled concurrently: 1) corrupted data estimation when only partial features are accessible; 2) online feature selection when data contains adversarial corruption; and 3) scaling to a massive dataset. This paper proposes a novel RObust regression algorithm via Online Feature Selection (\textit{RoOFS}) that concurrently addresses all the above challenges. Specifically, the algorithm iteratively updates the regression coefficients and the uncorrupted set via a robust online feature substitution method. We also prove that our algorithm has a restricted error bound compared to the optimal solution. Extensive empirical experiments in both synthetic and real-world datasets demonstrated that the effectiveness of our new method is superior to that of existing methods in the recovery of both feature selection and regression coefficients, with very competitive efficiency.

Figures

Tables

Table 1. TABLE I: Math Notations

Notations	Explanations
$p,n\in ℝ$	number of entire features and data samples
$p_t\in ℝ$	number of features in $t\text{<sup>th</sup>}$ time interval
$\mu \in ℝ$	ratio of feature sparsity, where ${\parallel 𝜷\parallel }_0=\mu$
$X_t\in {ℝ}^{p_t\times n}$	data samples containing features in the $t\text{<sup>th</sup>}$ time interval
$X\in {ℝ}^{p\times n}$	data samples containing the entire features
$𝜷,{𝜷}^{\ast }\in {ℝ}^{p\times 1}$	estimated and ground truth regression coefficient
$𝒖\in {ℝ}^{n\times 1}$	corruption vector with adversarial values
$𝜺\in {ℝ}^{n\times 1}$	dense noise vector, where ${\epsilon }_i\sim 𝒩(0,{\sigma }^2)$
$𝒚\in {ℝ}^{n\times 1}$	response vector, where $𝒚=X^T{𝜷}^{\ast }+𝒖+𝜺$
$𝒓\in {ℝ}^{n\times 1}$	residual vector, where $𝒓=\left|𝒚-X^T𝜷\right|$
$S\subseteq \left[n\right]$	estimated uncorrupted set
$S_{\ast }\subseteq \left[n\right]$	ground truth uncorrupted set, where $S_{\ast }=\overline{supp\left(𝒖\right)}$
$\Psi ,{\Psi }_{\ast }\subseteq \left[\mu \right]$	estimated and ground truth feature set

Table 2. TABLE II: F1 Scores for the Performance on Uncorrupted Set Recovery.

	p=2K, n=1K, $\mu$/$p$=20%				p=2K, n=2K, $\mu$/$p$=20%				p=4K, n=2K, $\mu$/$p$=20%
	10%	20%	30%	40%	10%	20%	30%	40%	10%	20%	30%	40%
Homotopy	0.980	0.912	0.849	0.682	0.977	0.923	0.854	0.834	0.970	0.923	0.845	0.775
DALM	0.976	0.915	0.865	0.825	0.973	0.921	0.885	0.924	0.962	0.946	0.926	0.898
TORR*	0.983	0.950	0.927	0.893	0.983	0.960	0.919	0.934	0.978	0.954	0.934	0.916
TORR25	0.965	0.899	0.842	0.762	0.961	0.909	0.828	0.770	0.958	0.905	0.848	0.752
RLHH	0.979	0.945	0.933	0.901	0.978	0.966	0.936	0.914	0.980	0.959	0.940	0.896
RoOFS	0.991	0.986	0.974	0.933	0.993	0.991	0.976	0.946	0.993	0.988	0.975	0.923

Table 3. TABLE III: F1 Score on Performance of Feature Selection (cr=30%).

	p=10K, n=10K, dense						p=20K, n=10K, dense
	10%	20%	30%	40%	50%	60%	10%	20%	30%	40%	50%	60%
Grafting	0.130	0.201	0.302	0.543	0.759	0.789	0.111	0.399	0.642	0.773	0.844	0.895
OS	0.706	0.649	0.626	0.643	0.810	0.678	0.611	0.611	0.606	0.689	0.836	0.975
Homotopy	0.116	0.217	0.304	0.406	0.498	0.667	0.109	0.202	0.417	0.625	0.750	0.833
DALM	0.130	0.219	0.309	0.418	0.516	0.597	0.114	0.215	0.390	0.404	0.502	0.603
TORR	0.275	0.297	0.368	0.446	0.525	0.647	0.320	0.338	0.395	0.458	0.535	0.623
RLHH	0.193	0.261	0.338	0.416	0.510	0.647	0.322	0.336	0.390	0.461	0.536	0.628
RoOFS	0.911	0.910	0.876	0.870	0.895	0.891	0.876	0.842	0.832	0.848	0.881	0.906

Table 4. TABLE IV: Mean Absolute Error of Sentiment Prediction.

	p=10K, n=10K
	5%	10%	20%	30%	40%	Avg
Grafting	3.162	3.162	3.162	3.162	3.162	3.162
OS	3.111	3.078	3.078	3.113	3.108	3.098
Homotopy	4.157	3.925	3.894	3.828	3.540	3.869
DALM	3.573	3.311	3.523	3.459	3.252	3.424
TORR	5.090	3.928	4.596	5.147	4.218	4.596
RLHH	5.448	4.198	3.716	4.269	4.626	4.451
RoOFS	3.074	3.073	3.072	3.070	3.066	3.071

Equations

$$\mathbf{y}=X^T{\mathbf{\beta }}^{\ast }+\mathbf{u}+\mathbf{\epsilon }$$

$$\begin{array}{c}\hat{\mathbf{\beta }},\hat{S}=\underset{\mathbf{\beta },S}{\mathrm{argmin}}\Vert {\mathbf{y}}_S-X_S^T\mathbf{\beta }{\Vert }_2^2\\ s.t.S\subset [n],|S|\ge \mathcal{G}(\mathbf{\beta }),\Vert \mathbf{\beta }{\Vert }_0\le \mu \end{array}$$

$${\mathcal{H}}_{\tau }(\mathbf{v})=\{i\in [n]:{\phi }_{\mathbf{v}}^{-1}(i)\le \tau \}$$

$$\mathbf{r}=|\mathbf{y}-X^T\mathbf{\beta }|$$

$$\hat{\tau }:=$$

$$f_S({\mathbf{\beta }}_1)-f_S({\mathbf{\beta }}_2)\ge {\nabla }^T{f}_S({\mathbf{\beta }}_2)({\mathbf{\beta }}_1-{\mathbf{\beta }}_2)+\frac{{\phi }_{\mu }}{2}\Vert {\mathbf{\beta }}_1-{\mathbf{\beta }}_2{\Vert }_2^2$$

$$\displaystyle r_{\tau}^{2}\leq\bigg{(}2\cdot\frac{\tau r_{\tau_{o}}}{\tau_{o}}\bigg{)}^{2}\ \stackrel{{\scriptstyle(a)}}{{\leq}}$$

$$|S_t\setminus S_{\ast }|r_{\tau }^2\stackrel{(b)}{\le }$$

$$f_{S_t\setminus S_{\ast }}(\mathbf{\beta })\le$$

$$+$$

$$f_{S_t\setminus S_{\ast }}(\mathbf{\beta })+f_{S_{\ast }\cap S_t}(\mathbf{\beta })\stackrel{(c)}{\le }$$

$$f_{S_t}(\mathbf{\beta })\stackrel{(d)}{\le }$$

$$\displaystyle f_{\hat{S}}(\hat{\bm{\beta}})-f_{S^{*}}(\bm{\beta}^{*})\leq\frac{\alpha\lambda}{1+\alpha}f_{S^{*}}(\bm{0})+\bigg{(}\frac{\lambda}{1+\alpha}-1\bigg{)}f_{S^{*}}(\bm{\beta}^{*})$$

$$f_S({\mathbf{\beta }}_1)-f_S({\mathbf{\beta }}_2)\ge ⟨\nabla f_S({\mathbf{\beta }}_2),{\mathbf{\beta }}_1-{\mathbf{\beta }}_2⟩+\frac{{\phi }_{\mu }}{2}\Vert {\mathbf{\beta }}_1-{\mathbf{\beta }}_2{\Vert }_2^2$$

$$f_S({\mathbf{\beta }}_1)$$

$$\ge$$

$$\stackrel{(a)}{=}$$

$$\ge$$

$$f_{S^{\ast }}(\hat{\mathbf{\beta }})-f_{S^{\ast }}({\mathbf{\beta }}^{\ast })\le$$

$$f_{S^{\ast }}(0)-f_{S^{\ast }}(\hat{\mathbf{\beta }})\ge ⟨\nabla f_{S^{\ast }}(\hat{\mathbf{\beta }}),-\hat{\mathbf{\beta }}⟩+\frac{{\phi }_{\mu }}{2}\Vert \hat{\mathbf{\beta }}{\Vert }_2^2$$

$$\mu {\hat{\mathbf{\beta }}}_{\min }^2\le \Vert \hat{\mathbf{\beta }}{\Vert }_2^2\le$$

$$\displaystyle 2\varphi_{\mu}\cdot\eta^{2}\Big{[}f_{S^{*}}(\hat{\bm{\beta}})-f_{S^{*}}(\bm{\beta}^{*})\Big{]}\stackrel{{\scriptstyle(c)}}{{\leq}}$$

$$\displaystyle f_{S^{*}}(\hat{\bm{\beta}})-f_{S^{*}}(\bm{\beta}^{*})\leq\alpha\Big{[}f_{S^{*}}(\bm{0})-f_{S^{*}}(\hat{\bm{\beta}})\Big{]}$$

$$\displaystyle f_{S^{*}}(\hat{\bm{\beta}})\leq\frac{\alpha}{1+\alpha}\Big{[}f_{S^{*}}(\bm{0})-f_{S^{*}}(\hat{\bm{\beta}})\Big{]}+f_{S^{*}}(\bm{\beta}^{*})$$

$$\displaystyle\frac{1}{\lambda}f_{\hat{S}}(\hat{\bm{\beta}})\leq f_{S^{*}}(\bm{\beta}^{*})\leq\frac{\alpha}{1+\alpha}\Big{[}f_{S^{*}}(\bm{0})-f_{S^{*}}(\hat{\bm{\beta}})\Big{]}+f_{S^{*}}(\bm{\beta}^{*})$$

$$\displaystyle f_{\hat{S}}(\hat{\bm{\beta}})-f_{S^{*}}(\bm{\beta}^{*})\leq\frac{\alpha\lambda}{1+\alpha}f_{S^{*}}(\bm{0})+\bigg{(}\frac{\lambda}{1+\alpha}-1\bigg{)}f_{S^{*}}(\bm{\beta}^{*})$$

Full text

I Introduction

The presence of noise and data corruption in real-world data can be inevitably caused by various reasons such as experimental errors, accidental outliers, or even adversarial data attacks. In traditional robust regression problem, reliable regression coefficients are learned in the presence of adversarial data corruptions in its response vector. A commonly adopted model from existing methods assumes that the observed response is obtained from the generative model $\bm{y}=X^{T}\bm{\beta}^{*}+\bm{u}$, where $\bm{\beta}^{*}$ is the true regression coefficients we wish to recover and $\bm{u}$ is the corruption vector with adversarial values. In the problem setting, the data matrix $X$ is assumed to contain all the features that can be accessed at any time and by arbitrarily many times.

Existing robust learning methods typically focus on modeling the entire dataset with all the features at once; however, they may meet the bottleneck in terms of computation and memory as more and more data sets are becoming However, the assumption is no longer suitable to the following scenarios in the applications that contain exponentially increasing user-generated contents:

1. features are too many to be loaded entirely. Features grows dramatically fast and becomes extremely large in. For instance, over 4.7 million movies and televisions with 8.3 million reviews in IMDb111https://www.imdb.com/ online movie and television review website, which makes it hard to load all the features entirely for any machine learning models using the movies as features.
2. features are generated dynamically. For example, people create and use new terms and hashtags all the time in Twitter, and the ”Likes” [1] on newly-generated articles in Facebook can be considered as new features describing the interestingness of the user.

Therefore, it is necessary to address online features in traditional robust regression as a new fundamental problem; however, current methods either focus on robust regression or online feature learning separately.

To the best of our knowledge, our proposed approach is the first robust regression algorithm that can handle the online features with adversarial data corruptions. It is nontrivial to consider online features and adversarial data corruption simultaneously in robust regression because 1) robust methods usually estimate data corruption based on the entire data, but online features make the data can only be partially accessible at one time; and 2) online feature selection methods can only select features based on uncorrupted data. Simply using robust regression and online feature selection methods sequentially makes the recovery result of coefficients worse, which is presented in our experiments in Section VI. To address the above challenges, we proposed a new robust regression algorithm via online feature selection (RoOFS). The main contributions of our study are summarized as follows:

• •

design of an efficient algorithm to simultaneously address the problem of data corruption and online feature. The algorithm RoOFS is proposed to recover the regression coefficients and uncorrupted set efficiently. Unlike using entire features, our approach alternately estimates the data corruption and selects the feature set via a robust online feature substitution method.

• •

theoretical analysis of the algorithm. We prove that our method yields a solution with a restricted error bound compared to ground truth coefficients under the Subset Restricted Strong Convexity (SRSC) property.

• •

demonstration of empirical effectiveness and efficiency. Our proposed algorithm was evaluated with 6 competing methods in both robust regression and online feature selection literatures. The results showed that our approach consistently outperforms existing methods in coefficients recovery and uncorrupted set estimation, delivering a competitive running time.

The reminder of this paper is organized as follows. Section II reviews the related work in robust regression model and online feature selection categories. Section III gives a formal problem formulation. The proposed RoOFS algorithm is presented in Section IV. Section V presents the theoretical analysis of proposed algorithm. In Section VI, the experimental results are analyzed and the paper concludes with a summary of our work in Section VII.

II Related Work

The work related to this paper is summarized in the categories of robust regression model and online feature selection as below.

II-A Robust Regression Model

A large body of literature on robust regression problem has been established over the last few decades. Most of studies focus on handling stochastic noise in small amounts [2]; however, these methods cannot be applied to data that may exhibit malicious corruption [3]. To recover regression coefficients with adversarial data corruption, Chen et al. [3] proposed a robust algorithm based on trimmed inner product. McWilliams et al. [4] proposed a sub-sampling algorithm for large-scale corrupted linear regression, but their theoretical recovery boundaries are not close to the ground truth [5]. Some $L_{1}$ penalty based methods [6, 7] pursue strong recovery results for robust regression problem, but these methods depend on severe restrictions of the data distribution such as row-sampling from an incoherent orthogonal matrix [7]. Zhang et al. [8] proposed a distributed robust algorithm to handle the large-scale data set under adversarial data corruption.

Most research in this area requires the corruption ratio parameter, which is difficult to estimate under the assumption that the dataset can be adversarially attacked. For instance, She and Owen [9] rely on a regularization parameter to determine the size of the uncorrupted set based on soft-thresholding. Chen et al. [3] require the upper bound of the outliers number, which is also difficult to estimate when the data contain the adversarial data corruption. Bhatia et al. [5] proposed a hard-thresholding algorithm with a strong guarantee of coefficient recovery under mild assumption on input data. However, the corruption ratio parameter is required by the algorithm and its recovery error can be more than doubled in size if the parameter is far from the true value. Recently, Zhang et al. [10] proposed a heuristic hard-thresholding based methods that learns the optimal uncorrupted set. However, all these approaches are based on batch feature selection under the assumption that all features can be accessed entirely at any time, which is infeasible to apply in massive and fast growing feature set.

II-B Online Feature Selection

Online feature selection methods [11, 12, 13] relaxes the requirement of batch selection and fit the scenarios that feature cannot be accessed entirely at one time. Statistical online feature selection algorithms [14, 15, 16] select features via certain statistical quantity such as mutual information, but these methods lack of specific objectives and usually have sub-optimal solutions for some certain tasks. Optimization based approaches [17, 18] use target oriented objective functions solved by some specific optimization techniques. These methods usually require the regression coefficient $\bm{\beta}$ be sparse, i.e., $\lVert\bm{\beta}\rVert_{0}\leq\mu$. Grafting [19] and its variation [18] relax the hard constraint of feature set into $L_{1}$ penalty, which makes it a convex problem. However, the parameter of $L_{1}$ norm [20] is difficult to determine because the usual cross validation strategy is unavailable for the online feature selection scenario [21]. Yang et al. [22] proposed a limited-memory substitution algorithm based on the $L_{0}$ norm constraint. Although the hard constraint leads to an NP-hard problem, a theoretical guarantee for the error bound of their local optimal solution is provided. However, none of these online feature methods can handle the adversarial data corruption.

III Problem Formulation

In this study, we consider the problem of robust regression with adversarial data corruption in the feature selection scenario in which only a few features are accessible at each time. Given data matrix $X_{t}\in\mathbbm{R}^{{p_{t}}\times n}$ where $p_{t}$ is the number of features available in the $t\textsuperscript{th}$ time interval, and $n$ are the number of data samples. The data matrix for all the time intervals is represented as $X=\{X_{t}\}_{t=1}^{\mathcal{T}}$. We assume the corresponding response vector $\bm{y}\in\mathbbm{R}^{n\times 1}$ is generated using the following model:

[TABLE]

where $\bm{\beta}^{*}$ represents the $\mu$-sparse ground truth coefficients of the regression model i.e., $\lVert\bm{\beta}^{*}\rVert_{0}\leq\mu$ and $\bm{u}$ is the unbounded corruption vector introduced by adversarial data attacks. $\bm{\varepsilon}\in\mathbbm{R}^{n\times 1}$ represents the additive dense noise, where $\varepsilon_{i}\sim\mathcal{N}(0,\sigma^{2})$. Different from the corruption vector $\bm{u}$ that can be arbitrarily distributed, the dense noise $\varepsilon_{i}$ follows normal distribution with zero mean and a relatively small variance $\sigma$. The notations used in this paper is summarized in Table I.

The goal of our problem is to learn a new robust regression problem with online feature selection, which is to recover the regression coefficients $\bm{\beta}^{*}$ and simultaneously determine the uncorrupted point set $\hat{S}$ with sequentially accessible features. The problem is formally defined as follows:

[TABLE]

Given a subset $S\subset[n]$, $\bm{y}_{S}$ restricts the row of $\bm{y}$ to indices in $S$ and $X_{S}$ signifies that the columns of $X$ are restricted to indices in $S$. Therefore, we have $\bm{y}_{S}\in\mathbbm{R}^{|S|\times 1}$ and $X_{S}\in\mathbbm{R}^{p\times|S|}$. We use the notation $S_{*}=\overline{supp(\bm{u})}$ to denote the ground truth set of uncorrupted points. Also, for any vector $\bm{v}\in\mathbbm{R}^{n}$, the notation $\bm{v}_{S}$ represents the $|S|$-dimensional vector containing the components in $S$. The notation $\Psi=supp(\bm{\beta})$ is used to represent the set of selected features, resulting in $\lvert\Psi\rvert\leq\mu$. Similarly, we use $X_{\Psi}$ to signify the rows of $X$ are restricted to indices in $\Psi$ and $X_{\Psi,S}$ to restrict both the rows and columns in set $\Psi$ and $S$. The function $\mathcal{G}(\cdot)$ determines the size of uncorrupted data according to the regression coefficients $\bm{\beta}$, which is explained in Section IV. It is worth mentioning that the features of data matrix $X$ in Equation (2) cannot be loaded entirely, but they can be accessed partially for each time interval. Therefore, the joint optimization of $\bm{\beta}$ and $S$ in our problem are very challenging because it amounts to a non-convex discrete optimization problem under the assumption that data matrix $X$ cannot be access entirely at one time.

IV The Proposed Methodology

To solve the problem in Equation (2) efficiently with the guarantee on the strong recovery of regression coefficients, we propose a novel robust regression algorithm with online feature selection, RoOFS. The algorithm is only allowed to access part of features at each time, which are defined as the newly incoming features in our problem. One naive solution to handle the sequentially incoming features is to retain all the features in the memory and then apply traditional robust feature selection methods. However, the solution has two major drawbacks: 1) the feature set can be too large to be retained in the memory, and 2) the algorithm becomes slower and slower when the feature set increases. Therefore, we proposed a new “robust online substitution” method to decide the retained feature set based on an adaptively estimated corrupted set. The procedure of robust online substitution is defined as follows:

• •

Update coefficients of retained features $\Psi$ based on the estimated uncorrupted set $S$ as follows: \bm{\beta}_{\Psi}:=$$\bm{\beta}_{\Psi}-\eta X_{\Psi,S}^{T}$$(X_{\Psi,S}^{T}\bm{\beta}_{\Psi}-\bm{y}_{S}), where $\eta$ is the step length.

• •

Retain the top $\mu$ largest (in magnitude) elements in $\bm{\beta}$ and set the rest to zero. Then all the non-zero features will be kept in the retained feature set $\Psi$.

• •

Compute the residual vector $\bm{r}\in\mathbbm{R}^{n\times 1}$ with the updated coefficients $\bm{\beta}$, then estimate the uncorrupted feature set $S$ via a thresholding operator $\mathcal{H}_{\tau}(\bm{r})$, where $\tau$ is the estimated size of uncorrupted set.

The procedure will be repeatedly executed until the residual vector $\bm{r}$ converges. The thresholding operator $\mathcal{H}_{\tau}(\cdot)$ is formally defined as follows:

Definition 1** (Thresholding Operator).**

Defining $\varphi^{-1}_{\bm{v}}(i)$ as the position of the $i$th element in input vector $\bm{v}$’s ascending order of magnitude and $\tau$ as the threshold parameter, the thresholding operator of $\bm{v}$ is defined as

[TABLE]

To estimate the uncorrupted set $S$, the thresholding operator $\mathcal{H}_{\tau}(\cdot)$ generally requires two inputs: residual vector $\bm{r}$ and the size of uncorrupted set $\tau$. The residual vector $\bm{r}$ can be computed with coefficients $\bm{\beta}$ as follows:

[TABLE]

For the size of uncorrupted set, two general cases are discussed. The first case is that the size can be estimated by users based on their prior knowledge on the data. For instance, if we know the data corruption happens rarely, then we can estimate the uncorrupted size as 95% of the entire data. However, it is hard to obtain prior knowledge on the data in the real-world. Thus, in the second case, we propose a method to adaptively estimate the uncorrupted size based on the residual vector $\bm{r}$. The method follows an intuition that when the coefficient $\bm{\beta}$ is close to $\bm{\beta}^{*}$, the residuals of uncorrupted samples are smaller than those of corrupted samples in strong possiblity. The intuition can be explained by the generative model in Equation (1), where the corrupted samples have the residual $\bm{r}\approx\bm{u}+\bm{\varepsilon}$, but the residual of uncorrupted samples only contains the white noise $\bm{\varepsilon}$.

The estimation of uncorrupted size can be formalized to solve the following problem:

[TABLE]

where $r_{\varphi(k)}$ represents the $k\textsuperscript{th}$ elements of residual vector $\bm{r}$ in ascending order of magnitude. The variable $\tau_{o}$ in the constraint is defined as an intermediate variable whose $r_{\varphi(\tau_{o})}^{2}$ has the closest value to $\frac{\lVert\bm{r}_{\mathcal{H}_{\tau^{\prime}}(\bm{r})}\rVert_{2}^{2}}{\tau^{\prime}}$, where $\tau^{\prime}=\tau-\left\lceil n/2\right\rceil$ and $\mathcal{H}_{\tau^{\prime}}(\bm{r})$ represent the position set containing the smallest $\tau^{\prime}$ elements in residual $\bm{r}$. The problem in Equation (5) can be solved by searching from $n$ to $\left\lceil n/2\right\rceil+1$ and return the first value $\hat{\tau}$ which satisfies the constraint. It is important to note that the estimation method in Equation (5) requires the coefficients $\bm{\beta}$ to be close to $\bm{\beta}^{*}$. Thus, we optimize the uncorrupted set $S$ along with coefficient $\bm{\beta}$ until both of them converge.

The details of RoOFS algorithm are presented in Algorithm 1. In Line 3, the algorithm receives data matrix $X_{\Psi^{k}}$ with the incoming feature set $\Psi^{k}$ at time $k$. The new feature set $\Psi^{k}$ is combined into the retained feature set $\Psi$ in Line 4. For each incoming feature set, the algorithm iteratively optimizes the regression coefficients $\bm{\beta}$ and the uncorrupted set $S$ until the value of residual vector $\bm{r}_{S_{t}}^{t}$ is converged in Line 14. Specifically, in Line 6, regression coefficients $\bm{\beta}$ are updated to a better fit for the current estimated feature set $\Psi$ and uncorrupted set $S_{t}$. In Line 8, feature set $\Omega$ that contains features with $\lvert\Psi\rvert-\mu$ smallest weights in $\bm{\beta}$ is selected. Then features in $\Omega$ are removed from the retained feature set $\Psi$ and the weights in $\bm{\beta}_{\Omega}$ are reset to zero in Lines 9 and 10. The residual vector $\bm{r}$ is updated in Line 11, while the uncorrupted set $S_{t+1}$ is estimated in Line 12 by the thresholding operator. Finally, both coefficients $\bm{\beta}$ and uncorrupted set $S$ are returned in Line 17.

V Theoretical Analysis

In this section, we show that the local optimal solution of our algorithm obtains a restricted error bound compared to ground truth solution. To prove the theoretical properties of our algorithm, we require that the least squares function satisfies the Subset Restricted Strong Convexity (SRSC), which is defined as follows:

Definition 2** (SRSC Property).**

The least squares function $f_{S}(\bm{\beta})=\lVert\bm{y}_{S}-X_{S}^{T}\bm{\beta}\rVert_{2}^{2}$ satisfies Subset Restricted Strong Convexity (SRSC) Property if the following holds for $\forall\bm{\beta}_{1},\bm{\beta}_{2}\in\Omega_{\mu}$ and $\forall S\in\mathcal{S}_{\gamma}$:

[TABLE]

To provide the local optimality property of our solution, the following two lemmas are first proved.

Lemma 1**.**

For a given least squares function $f(\bm{\beta})=\lVert\bm{y}-X^{T}\bm{\beta}\rVert_{2}^{2}$, let residual vector $\bm{r}=\bm{y}-X^{T}\bm{\beta}$ and $\delta(k)$ be the k-th position of the ascending order in vector $\bm{r}$, i.e. $r_{\delta(1)}\leq r_{\delta(2)}\leq...\leq r_{\delta(n)}$. For any $1\leq\tau_{1}<\tau_{2}\leq n$ and $\forall\bm{\beta}^{t}\in\Omega_{m}$, let $S_{1}=\{\delta(i)|1\leq i\leq\tau_{1}\}$ and $S_{2}=\{\delta(i)|1\leq i\leq\tau_{2}\}$. We then have $f_{S_{1}}(\bm{\beta}^{t})\leq f_{S_{2}}(\bm{\beta}^{t})$.

Proof.

Let $S_{3}=\{\delta(i):\tau_{1}+1\leq i\leq\tau_{2}\}$. Clearly, we have $f_{S_{2}}(\bm{\beta}^{t})=f_{S_{1}}(\bm{\beta}^{t})+f_{S_{3}}(\bm{\beta}^{t})$. Moreover, since each element in $S_{3}$ is larger than any of the element in $S_{1}$, we have $f_{S_{1}}(\bm{\beta}^{t})\leq f_{S_{2}}(\bm{\beta}^{t})+\frac{|S_{3}|}{|S_{1}|}f_{S_{1}}(\bm{\beta}^{t})$ $\leq\frac{|S_{1}|}{|S_{1}|+|S_{3}|}f_{S_{2}}(\bm{\beta}^{t})=\frac{\tau_{1}}{\tau_{2}}f_{S_{2}}(\bm{\beta}^{t})\leq f_{S_{2}}(\bm{\beta}^{t})$. ∎

Lemma 2**.**

Let $\tau_{*}=\gamma n$ be the true number of uncorrupted samples and $\tau_{t}$ be the estimated uncorrupted threshold at the $t$-th iteration. If $\tau_{t}\leq\tau_{*}$, then $f_{S_{t}}(\bm{\beta}^{t})\leq f_{S_{*}}(\bm{\beta}^{t})$. If $\tau_{t}>\tau_{*}$, then $f_{S_{t}}(\bm{\beta}^{t})\leq\lambda f_{S_{*}}(\bm{\beta}^{t})$, where \lambda=\Big{[}1+\frac{128(1-\gamma)}{2\gamma-1}\Big{]}.

Proof.

To simplify the notation, the subscripts $t$ that signify the $t$-th iteration will be omitted and the residual vector $\bm{r}$ is assumed to be sorted in ascending order of magnitude.

We will discuss the $\tau_{t}$ value in two different conditions compared to the value of $\tau_{*}$. In the first condition that $\tau_{t}\leq\tau_{*}$, let $S_{t}=\{\delta(i)|1\leq i\leq\tau_{t}\}$ and $S_{*}=\{\delta(i)|1\leq i\leq\tau_{*}\}$, we have $f_{S_{t}}(\bm{\beta}^{t})\leq f_{S_{*}}(\bm{\beta}^{t})$ according to Lemma 1. When $\tau_{t}>\tau_{*}$, we have the following properties according to the constraint specified in equation (5).

[TABLE]

The inequality (a) follows the definition of $\tau_{o}$ and the fact that $\lvert S_{*}\cap S_{t}\rvert\geq\tau^{\prime}$. The inequality (b) follows $\lvert S_{t}\setminus S_{*}\rvert\leq(1-\gamma)\cdot n$ and $\lVert\bm{r}_{S_{t}\setminus S_{*}}\rVert_{2}^{2}\leq\lvert S_{t}\setminus S_{*}\rvert r_{\tau}^{2}$. Then we have

[TABLE]

The inequality (c) follows $f_{S_{*}}(\bm{\beta})=f_{S_{*}\setminus S_{t}}(\bm{\beta})+f_{S_{*}\cap S_{t}}(\bm{\beta})$ and the inequality (d) follows $\tau^{\prime}=\tau_{t}-\frac{n}{2}$. ∎

Theorem 3**.**

Assume that least squares function $f_{S}(\bm{\beta})=\lVert\bm{y}_{S}-X_{S}^{T}\bm{\beta}\rVert_{2}^{2}$ satisfies Subset Restricted Strong Convexity (SRSC) Property for $\forall\bm{\beta}_{1},\bm{\beta}_{2}\in\Omega_{\mu}$ and $\forall S\in\mathcal{S}_{\gamma}$, then we have

[TABLE]

where \alpha=\big{(}\frac{1}{\eta\cdot\varphi_{\mu}}\big{)}^{2} and \lambda=\Big{[}1+\frac{128(1-\gamma)}{2\gamma-1}\Big{]}. Specifically, when the uncorrupted set size is less than ground truth, $\lambda=1$.

Proof.

As function $f_{S}(\bm{\beta})$ satisfies SRSC property, we have

[TABLE]

for $\forall\bm{\beta}_{1},\bm{\beta}_{2}\in\Omega_{\mu}$ and $\forall S\in\mathcal{S}_{\gamma}$. Let $\operatorname*{supp}(\bm{\beta}_{1})=\Omega_{1}$, then we have

[TABLE]

The equation (a) solves the minimum value of $\min_{\operatorname*{supp}(\bm{\beta})\subseteq\Omega_{1}}\{\cdot\}$ by setting its gradient to zero, and the inequality (b) follows $\lvert\Omega_{1}\rvert\leq\mu$. Let $\bm{\beta}_{1}$, $\bm{\beta}_{2}$ be ground truth coefficient $\bm{\beta}^{*}$ and estimated solution $\hat{\bm{\beta}}$ respectively, and set $S$ be ground truth uncorrupted set $S^{*}$, then we have

$\displaystyle f_{S^{*}}(\hat{\bm{\beta}})-f_{S^{*}}(\bm{\beta}^{*})\leq$ \displaystyle\frac{\mu}{2\varphi_{\mu}}{\Big{\|}\nabla f_{S^{*}}(\hat{\bm{\beta}})\Big{\|}}_{\infty}^{2}\leq \displaystyle\frac{\mu}{2\varphi_{\mu}}{\Big{(}\frac{1}{\eta}\hat{\bm{\beta}}_{\min}\Big{)}}^{2}

where $\hat{\bm{\beta}}_{\min}=\min_{i}\lvert\hat{\bm{\beta}}_{i}\rvert$. According to SRSC property, we have

[TABLE]

Because of $\langle\nabla f_{S^{*}}(\hat{\bm{\beta}}),-\hat{\bm{\beta}}\rangle\geq 0$, we have

[TABLE]

The inequality (c) follows the Equation (9). Let $\alpha$ be \big{(}\frac{1}{\eta\varphi_{\mu}}\big{)}^{2}, then we have

[TABLE]

According to Lemma 2, we have

[TABLE]

∎

Since the value of $f_{S^{*}}(\bm{0})$ and $f_{S^{*}}(\bm{\beta}^{*})$ are both constants and $f_{S^{*}}(\bm{\beta}^{*})$ is close to 0, the error bound of our solution is depended on the value of $\frac{\alpha\lambda}{1+\alpha}$. When the ratio of data corruption $\gamma$ is close to one, $\lambda$ is close to one according to its definition. In addition, the value of $\alpha$ is smaller when the parameter $\varphi_{\mu}$ of the SRSC property is smaller. Therefore, the error of our solution can be close to zero when both $\varphi_{\mu}$ and $\gamma$ are large enough.

VI Experimental Results

In this section, we report the extensive experimental evaluation performed to verify the robustness, effectiveness of feature selection, and efficiency of the proposed method. All the experiments were conducted on a 64-bit machine with Intel(R) core(TM) quad-core processor (i7CPU@3.6GHz) and 32.0GB memory. Details of both the source code and sample data used in the experiment can be downloaded here222https://goo.gl/C4HQjo.

VI-A Datasets and Metrics

To demonstrate the performance of our proposed method, comprehensive experiments are performed in synthetic datasets whose simulation samples were randomly generated according to the model in Equation (1). Specifically, we sample the regression coefficients $\bm{\beta}^{*}\in\mathbbm{R}^{p}$ as a random unit norm vector with feature ratio constraint $\lVert\bm{\beta}\rVert_{0}=\mu$. The data matrix $X$ was drawn independently and identically distributed from $\bm{x}_{i}\sim\mathcal{N}(\bm{0},I_{p})$ and the uncorrupted response variables were generated as $y_{i}^{*}=\bm{x}_{i}^{T}\bm{\beta}^{*}$. The set of uncorrupted samples $S$ was selected as a uniformly random $\tau_{*}$-sized subset of $[n]$. The response vector $\bm{y}$ containing corrupted samples was generated as $\bm{y}=\bm{y}^{*}+\bm{u}+\bm{\varepsilon}$, where the corruption vector $\bm{u}$ was sampled from the uniform distribution \big{[}-5\|\bm{y}^{*}\|_{\infty},5\|\bm{y}^{*}\|_{\infty}\big{]} and the additive dense noise was $\varepsilon_{i}\sim\mathcal{N}(0,\sigma^{2})$. For the real-world data set, we applied our methods on the IMDb reviews data set for the review score prediction. The data set contains 50,000 popular movie reviews with the review score from 1 to 10 provided by the IMDb website. The adversarial data corruption vector $\bm{u}$ was appended to its original review score, where $\bm{u}$ was also sampled from the range \big{[}-5\|\bm{y}^{*}\|_{\infty},5\|\bm{y}^{*}\|_{\infty}\big{]} randomly.

Following the setting in [5][10], we measured the performance of the regression coefficients recovery using the standard $L_{2}$ error $e=\lVert\hat{\bm{\beta}}-\bm{\beta}^{*}\rVert_{2}$, where $\hat{\bm{\beta}}$ represents the recovered coefficients for each method and $\bm{\beta}^{*}$ is the true regression coefficients. To validate the performance for corrupted set discovery, the F1 score is measured by comparing the discovered corrupted sets with the actual ones. Similarly, the F1 score is also used to measure the effectiveness of feature selection by comparing the selected feature set with actual ones. To compare the scalability of each method, the CPU running time for each of the competing methods was also measured.

VI-B Comparison Methods

The following methods are included in the performance comparison presented here: Grafting [19]. The Grafting method is an online version of $L_{1}$ regularization approach to selects features. Online Substitution (OS) [22] is a parameter-free online feature selection algorithm with limited-memory. Both Grafting and OS cannot handle the adversarial data corruption and train models without considering data corruption. We also compared our method to the robust regression methods [6] [7]. Homotopy and DALM are two $L_{1}$ based solvers that outperform other $L_{1}$ methods both in terms of recovery properties and running time [23]. A hard thresholding method, TORRENT (abbreviated ”TORR”) [5], developed for robust regression was also compared to our method. As the method requires a parameter for the corruption ratio, which is difficult to estimate in practice, we chose two versions of parameter settings: TORR* and TORR25. TORR* uses the true corruption ratio as its parameter, and TORR25 applies parameter that is uniformly distributed across the range of $\pm 25\%$ off the true value. Another recently proposed heuristic hard thresholding method, RLHH [10], is also compared in our experiment. The method is a parameter-free approach, where the data corruption is estimated by a heuristic hard thresholding method. As all these robust methods are not designed for online feature selection, we run them individually in different feature batches and select features with largest $\mu$ weights in regression coefficients when $\lVert\bm{\beta}\rVert_{0}=\mu$.

VI-C Recovery of regression coefficients

We selected 6 competing methods with which to evaluate the recovery performance of regression coefficients $\bm{\beta}$: Grafting, OS, Homotopy, DALM, TORR, RLHH. Figures 1(a) and 1(b) show the recovery performance for different feature numbers when the data size is fixed. The results show that 1) the proposed method, RoOFS, outperforms all the competing methods in all the setting of corruption ratios, and 2) The performance of RoOFS is very resistant to the corruption data because the error of RoOFS method increases much more slowly than others when corruption ratio increases from 5% to 40%. Figures 1(b) and 1(c) show that when data size increases, we have similar conclusion on the performance except the overall error is decreased since more data is applied. Figures 1(d) and 1(e) show that the result of coefficient recovery remains the same when the number of selected features increase. Figure 1(f) shows that almost all the methods without the dense noise setting perform more than 50% better than that in dense noise settings. Specifically, the error of RoOFS is close to zero which means it can almost exactly recover the ground true regression coefficients without the dense noise setting.

Figure 2 shows the recovery performance of regression coefficients in different ratios of feature sparsity. In general, the performance of RoOFS method outperforms all the other competing methods in all the data settings. Figure 2(a) and 2(b) show that 1) when the feature ratio increases, the recovery error of RoOFS method grows linearly with a small slope, which means our approach can be fitted into different settings of feature sparsity, and 2) the RoOFS method performs constantly well when the feature number increases from 2K to 4K. In addition, Figure 2(a) and 2(c) show that the RoOFS method is robust to corrupted data, because the error is not significantly impacted when the corruption ratio increases from 20% to 40%.

VI-D Recovery of Uncorrupted Set

As the online feature selection methods Grafting and OS do not explicitly estimate uncorrupted sets, we compared our proposed method with the robust methods: Homotopy, DALM, TORR, and RLHH. For the TORR algorithm, we use two parameter settings of TORR* and TORR25 for 0% and 25% deviation of true corrupted ratio, respectively. Table II shows the following: 1) RoOFS outperforms all the other methods up to 14.9% in different settings of data sizes, feature numbers and ratios of feature sparsity. 2) When increasing the corruption ratio, the F1 scores decrease for all the methods. Also, the F1 scores slightly increase 0.5% in average when data size become two times larger, which indicates the number of features has few influence on the estimation of uncorrupted set. 3) The result of TORR methods is highly dependent on the corruption ratio parameter: the results of TORR* is up to 26.1% better than TORR25. It is important to note that the true corruption ratio parameter used in TORR* cannot be estimated exactly in practice. 4) Without the dense noise settings, the F1 scores increase less than 1% compared to the F1 score based on the same setting with dense noise, which shows that dense noise has small impact on the performance of uncorrupted set recovery.

VI-E Performance of Feature Selection

We selected all the six competing methods to evaluate the performance of feature selection in different settings including data sizes, feature numbers, and dense noises. For each data setting, we chose different ratios of feature sparsity (also known as $\mu$/$p$) ranging from 10% to 60%. Table III shows the following: 1) the F1 scores of RoOFS method is up to 69.2% better than other methods, especially when the feature ratio is less than 40%. 2) Although the F1 scores of most methods such as Grafting and TORR are above 0.6 when the ratio is larger than 50%, the performance degraded significantly when the ratio decreased to 10%. However, the F1 score of RoOFS method is constantly higher than 0.85 in all the ratios of features. 3) OS method is very competitive in the task of feature selection; however, it still has lower F1 scores in all the settings when the ratio is less than 50%. 4) The setting of dense noise does not have significant impact on the performance of feature selection, since the F1 score without dense noise is only less than 1% larger than that in dense noise setting.

VI-F Performance in real-world data

To evaluate the robustness of our proposed methods in a real-world dataset, we compared the performance of sentiment prediction in different corruption settings, ranging from 5% to 40%. The dataset was first proposed by Maas et al. [24] as a benchmark for sentiment analysis. It consists of movie reviews taken from IMDB. One key aspect of this dataset is that each movie review has several sentences. The 100,000 movie reviews are divided into three datasets: 25,000 labeled training instances, 25,000 labeled test instances and 50,000 unlabeled training instances. The unlabeled data were designed as the additional corruption to the dataset: the score of sentiment were random number between one to ten. Table IV shows the mean absolute error of sentiment prediction in the IMDB datasets. From the result, we can conclude:

1. RoOFS method outperform all the other methods in different corruption settings.
2. Although the absolute error of the other methods such as Homotopy and DALM are above 4, the performance varied significantly when the ratio changed because these methods highly dependent on the parameters and it’s hard to estimate the feature sparsity ratio and true corruption ratio in the real-world data. However, the performance of RoOFS method is constantly above 3.10 in all the ratios of corruption.
3. It is true that OS has a very competitive performance in all the corruption settings because the deviation of corruption is small, which is less than 50% from the labeled data. But the running time of OS is too high to train the data which has 10k features.
4. When increasing the corruption ratio, the absolute error of RoOFS method decreased.

VI-G Efficiency

To evaluate the efficiency of our proposed method, we compared the performances of all the competing methods for two difference settings: data sizes and feature numbers. For Grafting and OS methods, the online features are handled individually due to their design. For the other methods, a hundred features are handled together as a batch. As Figure 3 shows, we found the following: 1) RoOFS algorithm has a very competitive efficiency compared to the thresholding based methods, TORR and RLHH, and significantly outperforms other four methods. 2) The running time of RoOFS algorithm increases linearly when both data size and feature number increase, which indicates that our algorithm can be scaled to massive datasets. 3) The running time of DALM method increases exponentially when data size increases, however, its efficiency has rarely impacted by increasing the number of features. 4) The efficiency of Grafting method fluctuates largely on the different data sizes, which indicates that its running time depends on the data size and content of data.

VII Conclusion

In this paper, a novel robust regression algorithm via online feature selection, RoOFS, is proposed to recover the regression coefficients and the uncorrupted set under the assumption that features cannot be accessed entirely at one time. To achieve this, we designed a robust online substitution method to alternately estimate the optimal uncorrupted set and substitute the retained feature set with newly updated features. We demonstrate that our algorithm can recover regression coefficients with a restricted error bound compared to ground truth. Extensive experiments on massive simulation data demonstrated that the proposed algorithm outperforms other competing methods in both effectiveness and efficiency.