Online Updating Huber Robust Regression for Big Data Streams

In recent years, big data has become a significant area of research and garnered extensive attention across various industries, including finance, manufacturing, and astronomy. Compared with traditional datasets, big data mainly has “5V” characteristics, i.e., high Volume, high Velocity, high Variety, low Veracity, and high Value(Jin et al., 2015). For big data streams, namely big data generated sequentially in a stream style, the characteristic of high Volume and high Velocity are particularly prominent, posing several challenges in its analysis. One such challenge concerns computer memory limitations. As big data streams are generated with such big volume, it is unrealistic to be completely stored in computer memory. Consequently, standard statistical analysis methods that assume adequate computer capacity are no longer applicable in the current settings. The second challenge pertains to the presence of outliers in big data streams, which can arise due to recording errors in the fast data-generating process. While linear regression-based methods are often used for big data analysis, they are highly sensitive to outliers, thereby limiting their effectiveness in this context. As such, it is imperative to develop computational efficient, memory-conserving and robust methods for big data stream analysis that can address these challenges.

To address the storage limitation problem in big data analysis, many scholars have designed different methodologies , which can be loosely divided into three categories, i.e. subsampling-based approaches, divide-and-conquer approaches and online updating approaches(Wang et al., 2016). According to Wang et al.(2016), classic subsampling methods include bags of little bootstrap(Kleiner et al., 2014), leveraging(Ping and Sun, 2015), mean log-likelihood(Faming et al., 2013) and subsample-based MCMC(Liang et al., 2016). Take Bags of Little Bootstrap as an example, Kleiner et al.(2014) combined bootstrap with subsampling, which addressed the high computation problem of bootstrap in big data occasions. As for divide-and-conquer, its most simple type, also called naïve-divide-and-conquer or one-shot, is gaining much popularity. Its main idea is dividing complete datasets into several blocks, obtaining the local estimators on each block and calculating the global one by simple averaging. Naïve-dc has many applications in regression. For example, to solve large-scale empirical risk minimization problem, Zhang et al.(2012) analyzed a communicational efficient average mixture algorithm based on naïve-dc which is robust to the amount of parallelization. Chen and Xie(2014) designed a split-and-conquer penalized regression, including $L_{1}$ penalty, SCAD and MCP. Their estimator is asymptotically equivalent to penalized estimator using entire datasets. Furthermore, some research also applied naïve-dc on high-dimensional big datasets. Lee et al.(2015) devised an averaging debiased LASSO estimator that can be obtained distributedly and communication efficient in high dimensional settings. Battey et al.(2018) cast their eyes on hypothesis testing. They utilized divide-and-conquer in the context of sparse high dimensional Rao and Wald test, solving computation complexity increase due to high dimensional regularization problems. We can also find some non-linear applications like ridge regression(Zhang et al., 2013), SVM(Lian and Fan, 2018) etc. More naïve-dc related literature can be seen in Nan and Xi(2011), Zhao et al.(2016),Zhu et al.(2021). However, complete datasets are necessary for naïve-dc before analyzing to find the appropriate block number. While this issue might not arise for conventional big datasets, it is a crucial challenge for big data streams since the latter are incessantly generated with significant volume and velocity, making it infeasible to obtain complete data.

Different from divide-and-conquer, online updating, aiming at big data streams analysis, has no requirement for storing historical data. Specifically, it uses key features extracted from new arriving data batches to continuously update historical data thus is able to achieve computational efficiency and minimal intensive storage. According to Schifano et al.(2016), they are the first group to investigate statistical inference on online updating scenario. Their online updating framework is able to solve the predictive residual tests in linear regression and also the estimating equation setting. However, this algorithm needs to start from scratch if new predictive variables are produced with data streams generating. To overcome this disadvantage, Wang et al.(2018) adjusted the online updating framework by bias correction. Through the use of nested models(Allison, 1995; Clogg et al., 1995) in regression coefficients comparison,their method improved estimation efficiency under a variety of practical situations. Lee et al.(2020) extended aforementioned work from the perspective of biased predictive variables in big data streams being corrected afterward under the linear model framework. However, their method is still biased if not knowing the exact data adjusted point. Recently, online updating has also been incorporated in other models for big data analysis, such as online updating GIC for variable selection in generalized linear model(Xue and Hu, 2021), Cox proportional hazards model in online updating framework(Wu et al., 2021), online nonparametric method for functional data(Fang Yao, 2021) etc. To conclude, online updating is accessible to not only ordinary big datasets, but also big data streams, thereby providing a wider range of applications than naïve-dc. Nevertheless, most existing online updating frameworks are combined with linear regression, which are known to be sensitive to outliers and heavy-tailed distributions. To address this issue, novel online updating frameworks that are designed to achieve robustness to outliers are urgently needed.

Talking about the outlier problem, one popular solution is quantile regression(Bassett and Jr., 1978). It studies the conditional distribution of each quantile of response variables, thus overcoming the outlier problem and being able to explain potential heteroscedasticity of the datasets. Inspired by its good performance, loads of research incorporate quantile regression into big data analysis framework, especially in divide-and-conquer. Take Chen et al.(2019) as an example, they used a kernel function to smooth quantile loss function and operated divide-and-conquer iteratively to obtain a linear quantile regression estimator, handling with the constraint of machine number in naïve-dc. Moreover, Hu et al.(2021) addressed the nondifferentiable problem of quantile regression by developing CSL(Jordan et al., 2019) method to speed up big data algorithm. We can also find other research combining quantile regression with subsampling and online updating like Chen and Zhou(2020) and Wang et al.(2022). By using online updating, Wang’s(2022) method can be applied to big data streams. Particularly, Wang et al.(2022) analyzed quantile regression from the view of maximum likelihood and incorporated it in online updating algorithm. Without operations in check loss function, it is more convenient than the work of Chen et al.(2019), and Chen and Zhou(2020). Some variants of quantile regression like composite quantile regression and adaptive quantile regression are also investigated in big data analysis, which can be seen in Jiang et al.(2018),Jiang et al.(2021) for details.

While quantile regression is investigated thoroughly, other members in M-estimation family, such as Huber regression(Huber and Peter, 1964), have received litte attention. Huber regression transforms the loss of big residuals from quadratic form into linear form to cut down weights, thus achieving good robustness to outliers. Therefore, Huber regression presents a feasible alternative to quantile regression in solving outlier-related problems. Regrettably, there has been little research on the use of Huber regression in big data analysis. To date, the only literature used Huber loss in big data regression is Luo’s(2022) work. They embedded Huber loss function in Jordan’s CSL(Jordan et al., 2019) framework, adding two data-driven parameters in Huber regression to achieve a balance of statistical optimization and communication efficiency. Nonetheless, there is no literature focusing on combining Huber regression with online updating to solve analysis in big data streams.

In order to solve the storage limitation and outlier problems aforementioned in big data streams analysis, we propose a novel Online Updating Huber Regression algorithm, which is inspired by the maximum likelihood idea in Wang’s(2022) online updating algorithm. The resulting online updating estimator in our algorithm is renewed by current data and summary statistics of historical data, which is computationally efficient. Specifically, by assuming local estimators generated from a multivariate normal distribution and analyzing it from the view of maximum likelihood, a weighted least square type objective function can be obtained. In the objective function, two key features of each batch can be extracted, which can be updated and easily operated in big data streams. Meanwhile, we prove that the proposed estimator is asymptotically equivalent to Oracle estimator using the entire dataset and has lower computation complexity compared with Oracle one. Numerical experiments on both simulation and real data are also included to verify the theoretical results and illustrate the good performance of our new method.

Our study makes two primary contributions. Firstly, we employ online updating instead of the widely used divide-and-conquer method to address the limited computer memory constraint so that it can be applied on the analysis of big data streams. Furthermore, by integrating Huber regression into online updating framework, we broaden the usage of online updating in robust regression other than pure linear regression. The second contribution is about Huber regression. Based on the state-of-art in big data analysis using loss functions in M-estimation, we adopt Huber regression as an alternative to quantile regression to increase robustness. It not only fills certain research gaps in Huber regression, but also studies feasibility of generalizing quantile regression to other loss functions in M-estimation.

The rest of this paper is organized as follows. Section 2 explains the methodology of our Online Updating Huber Regression by obtaining the main online updating estimator. Section 3 proves some good properties of the Online Updating Huber Regression estimator, including asymptotic equivalency and low computation complexity. Section 4 introduces simulation and real data analysis results in order to verify the good performance of our algorithm. Section 5 is about conclusions and future work of our research.

Suppose there are $N_{b}$ identically independent distributed samples $\{(y_{n},\textbf{x}_{n}),n=1,2,...,N_{b}\}$. $y\in\rm{R}$ is response variable and $\textbf{x}=(x^{1},x^{2},...,x^{p})^{T}\in\textbf{R}^{p\times 1}$ is p-dimensional covariates. In the big data streams scenario, suppose $N_{b}$ samples are generated in form of $b$ data subsets $\{D_{1},D_{2},...,D_{b}\}$ sequentially. $D_{t}=\{(y_{ti},\textbf{x}_{ti}),i=1,2,...,n_{t}\}$ stands for the $t$-th subset. Our aim is to fit a model $y_{ti}={\textbf{x}_{ti}}^{T}\bm{\theta}+{\epsilon_{ti}}$, where $\bm{\theta}=(\theta_{1},\theta_{2},...,\theta_{p})^{T}$ is a p-dimensional unknown parameter with true value $\bm{\theta}_{0}=(\theta_{01},\theta_{02},...,\theta_{0p})^{T}$, and $\epsilon_{ti}$ is random error which can be heavy-tailed distributed, heteroscedastic distributed or contaminated by outliers. $f(\cdot)$ is the probability density distribution of $\epsilon_{ti}$.

As we want to prove that our Online Updating estimator is asymptotic equivalent to Oracle one using entire dataset follow-up, we first point the expression of Oracle Huber Regression estimator here. According to Huber(1964), the expression of Oracle Huber Regression estimator of $\bm{\theta}_{0}$ is as follows:

where $\rho(u)=\frac{u^{2}}{2}I(|u|<k)+\{k|u|-\frac{k^{2}}{2}\}I(|u|\geq k)$.

Moreover, under following regularity conditions Huber and Peter (1964, 1981); Wang et al. (2022):

(N1) For $\lambda(\bm{\theta})=E(\textbf{x}\psi(y-\textbf{x}^{T}\bm{\theta}))$, where $\psi(\cdot)=\rho^{{}^{\prime}}(\cdot)$, $\bm{\theta}_{0}$, the true value of $\bm{\theta}$, satisfies $\lambda(\bm{\theta}_{0})=0$.

(N2)$E(\psi^{2})<\infty,0<E(\psi^{{}^{\prime}})<\infty$.

(N3) Covariates $\textbf{x}_{ti}$ satisfies $\frac{max_{t=1,...,b,i=1,...,n_{t}\|{\textbf{x}_{ti}}\|}}{\sqrt{N_{b}}}\rightarrow\infty,\|{\textbf{x}_{ti}}\|=({\textbf{x}_{ti}^{T}\textbf{x}_{ti}})^{\frac{1}{2}}$.

(N4)Exist a positive matrix $\bm{\Sigma}$ satisfying $\frac{1}{N_{b}}\sum_{t=1}^{b}\sum_{i=1}^{n_{t}}\textbf{x}_{ti}\textbf{x}_{ti}^{T}\xrightarrow{p}\bm{\Sigma},N_{b}\rightarrow\infty$, where $\xrightarrow{p}$ stands for convergence in probability.

We can obtain the asymptotic normal distribution of Oracle Huber Regression estimator $\widehat{\bm{\theta}}_{N_{b}}$ as follows:

where $\frac{E(\psi^{2})}{[E(\psi^{\prime})]^{2}}$ is a function of $f(\cdot)$, $\xrightarrow{d}$ stands for convergence in distribution.

In this part, we mainly deduce the expression of Online Updating Huber Robust Regression (UHR) estimator. Firstly, we can calculate the local Huber estimator in each subset. Suppose $\widehat{\bm{\theta}}_{t}$ is the local Huber estimator in the $t$-th subset. According to Huber(1964), its formula is like Eq.(3).

As the entire dataset is divided into several subsets, subset size $n_{t}$ can be much smaller than total sample size $N_{b}$. Therefore, the calculation of Eq.(3) can be more computationally efficient and less storage intensive than Eq.(1), so that we can solve it faster. Furthermore, under regularity conditions (N1)-(N4), local estimator $\widehat{\bm{\theta}}_{t}$ follows asymptotic normal distribution in Eq.(4), where $\bm{\Sigma}_{t}=E(\textbf{x}_{ti}\textbf{x}_{ti}^{T})$.

Taking the idea in Wang’s(2022) work as an inspiration, we consider the local estimators from the perspective of maximum likelihood. Due to identical independent distribution among total samples, local estimator $\widehat{\bm{\theta}}_{1},...,\widehat{\bm{\theta}}_{b}$ are mutually independent. Thus, we can take $\{\widehat{\bm{\theta}}_{1},...,\widehat{\bm{\theta}}_{b}\}$ as a sample generated from the following multivariate normal distribution.

Based on Eq.(5), the maximum likelihood function of $\bm{\theta}$ is as follows.

By turning Eq.(6) into logarithmic form, we can get Eq.(7)

where $C$ is a constant, $\frac{[E(\psi^{\prime})]^{2}}{E(\psi^{2})}$ is a function of $f(\cdot)$, which is irrelevant with $\bm{\theta}$ as well. By maximizing Eq.(7), we can obtain the UHR estimator. Meanwhile, base on non-negativity of $\frac{[E(\psi^{\prime})]^{2}}{E(\psi^{2})}$, maximizing Eq.(7) is equivalent to minimize the following loss function.

By using $\frac{1}{n_{t}}\sum_{i=1}^{n_{t}}\textbf{x}_{ti}\textbf{x}_{ti}^{T}$ to replace $\bm{\Sigma}_{t},t=1,...,b$ according to $\bm{\Sigma}_{t}=E(\textbf{x}_{ti}\textbf{x}_{ti}^{T})$, we can obtain the final objective function as follows, where $\textbf{X}_{t}=(\textbf{x}_{t1},...,\textbf{x}_{tn_{t}})^{T}$.

Thus Online Updating Huber Robust Regression estimator of $\bm{\theta}_{0}$, $\widehat{\bm{\theta}}_{N_{b}}^{uhr}$, can be obtained by minimizing Eq.(9).

According to Eq.(10), we can only extract two key features from each subset, i.e. $\textbf{X}_{t}^{T}\textbf{X}_{t}$ and $\widehat{\bm{\theta}}_{t}$, with no need to store entire dataset. This enables us to save computer memory and speed up computation process, thus more available to big data streams.

The UHR algorithm implementation process can be seen in Figure 1. To start the algorithm, we import the first subset $D_{1}=\{(y_{1i},\textbf{x}_{1i}),i=1,...,n_{1}\}$. By using function rlm of package MASS in R, the local estimator of the first subset $\widehat{\bm{\theta}}_{1}$ can be calculated. Then two features in the first subset $\textbf{X}_{1}^{T}\textbf{X}_{1}$ and $\textbf{X}_{1}^{T}\textbf{X}_{1}\widehat{\bm{\theta}}_{1}$ can be obtained. For the second subset $D_{2}=\{(y_{2i},\textbf{x}_{2i}),i=1,...,n_{2}\}$, use the same method to calculate $\widehat{\bm{\theta}}_{2}$ and $\textbf{X}_{2}^{T}\textbf{X}_{2}$, then add them into $\textbf{X}_{1}^{T}\textbf{X}_{1}$ and $\textbf{X}_{1}^{T}\textbf{X}_{1}\widehat{\bm{\theta}}_{1}$. The sum obtained, $\textbf{X}_{1}^{T}\textbf{X}_{1}+\textbf{X}_{2}^{T}\textbf{X}_{2}$ and $\textbf{X}_{1}^{T}\textbf{X}_{1}\widehat{\bm{\theta}}_{1}+\textbf{X}_{2}^{T}\textbf{X}_{2}\widehat{\bm{\theta}}_{2}$ are new results of subset $D_{1}$ updated by subset $D_{2}$. The rest subsets can be processed in the same manner. After the last subset $D_{b}=\{(y_{bi},\textbf{x}_{bi}),i=1,...,n_{1}\}$ is imported, we can obtain key elements, i.e. $\sum_{t=1}^{b}\textbf{X}_{t}^{T}\textbf{X}_{t}$ and $\sum_{t=1}^{b}\textbf{X}_{t}^{T}\textbf{X}_{t}\widehat{\bm{\theta}_{t}}$, and finally obtain the Online Updating estimator $\widehat{\bm{\theta}}_{N_{b}}^{uhr}$ updated by the entire dataset. The pseudo-code can be seen in Algorithm 1.

In order to prove the asymptotic normality of our UHR estimator, we first add some assumptions into regularity conditions (N1)-(N4)(Wang et al., 2022).

(N5) Take $n=\frac{N_{b}}{b}$ as the average size of each subset. Suppose $\frac{n}{\sqrt{N_{b}}}\rightarrow\infty$, and all $n_{t}$ diverge in the same order $O(n)$, i.e., $c_{1}\leq{\rm{min}}_{t}n_{t}/n\leq{\rm{max}}_{t}n_{t}/n\leq c_{2}$ for some positive constants $c_{1}$ and $c_{2}$.

(N6) Exist positive definite matrix $\bm{\Sigma}_{1},...,\bm{\Sigma}_{b}$ which satisfy $\frac{1}{n_{t}}\sum_{i=1}^{n_{t}}\textbf{x}_{ti}\textbf{x}_{ti}^{T}\xrightarrow{p}\bm{\Sigma}_{t},n_{t}\rightarrow\infty$ and converge to $\bm{\Sigma}_{t}$ in the same rate, i.e. convergence rate is irrelevant with $t$.

Obviously, $\frac{1}{n_{t}}\sum_{i=1}^{n_{t}}\textbf{x}_{ti}\textbf{x}_{ti}^{T}-\bm{\Sigma}_{t}=\textbf{O}(\frac{1}{\sqrt{n_{t}}})$ due to $\bm{\Sigma}_{t}=E(\textbf{x}_{ti}\textbf{x}_{ti}^{T})$. Moreover, according to (N5) and (N6), all subsets have the same convergence rate, which is $\textbf{O}(\frac{1}{\sqrt{n_{t}}})=\textbf{O}(\frac{1}{\sqrt{n}})$.

Under regularity conditions (N1)-(N6), UHR estimator $\widehat{\bm{\theta}}_{N_{b}}^{uhr}$ follows asymptotic distribution as follows:

where $\Phi(f)=\frac{E(\psi^{2})}{[E(\psi^{\prime})]^{2}}\Big{(}\sum_{t=1}^{b}\frac{n_{t}}{N_{b}}\bm{\Sigma}_{t}\Big{)}^{-1}$. The proof is as follows.

For local Huber estimator $\bm{\widehat{\theta}}_{t}$, it follows asymptotic condition:

where $\psi(\epsilon_{ti})=\rho^{{}^{\prime}}(\epsilon_{ti})=\left\{\begin{array}[]{ll}-k,&\epsilon_{ti}<-k\\ \epsilon_{ti},&\lvert\epsilon_{ti}\lvert\leq k\\ k,&\epsilon_{ti}>k\end{array}.\right.$

According to Eq.(10), the relationship with $\widehat{\bm{\theta}}_{N_{b}}^{uhr}$ and $\widehat{\bm{\theta}}_{t}$ is

where $w_{t}=\frac{n_{t}}{N_{b}}$. By direct calculation, we can obtain that

By $\sqrt{N_{b}}\Bigg{(}\sum_{t=1}^{b}w_{t}\Big{(}\frac{\textbf{X}_{t}^{T}\textbf{X}_{t}}{n_{t}}-\bm{\Sigma}_{t}\Big{)}(\widehat{\bm{\theta}_{t}}-\bm{\theta}_{0})\Bigg{)}=O_{p}\Big{(}\frac{b}{\sqrt{N_{b}}}\Big{)}$ and regularity conditions (N1)-(N6), we can obtain

As $\sum_{t=1}^{b}w_{t}=1$, we have

Then we can obtain that

If the covariates of each batch are homogeneous, i.e., $\Sigma_{1}=\cdots=\Sigma_{b}=\Sigma$, the $\Phi(f)=\frac{E(\psi^{2})}{[E(\psi^{\prime})]^{2}}\bm{\Sigma}^{-1}$, which coincides the asymptotic variance of oracle estimator in Eq.(2). The updating estimator $\widehat{\bm{\theta}}_{N_{b}}^{uhr}$ is asymptotically equivalent with the oracle estimator $\widehat{\bm{\theta}}_{N_{b}}$.

In this part, we do a rough calculation on the computation complexity of Oracle Huber Regression and Online Updating Huber Regression, attesting our Online Updating algorithm has a lower computation complexity.

In our work, we use function rlm in MASS package of R to calculate the local estimator and Oracle estimator obtained by entire dataset. Function rlm is implemented by Iteratively Reweighted Least Squares. For the sake of calculation, suppose rlm’s maximum iterative time is $K$, total sample size is $N$ and is divided into $b$ subsets evenly with $n$ samples each. Therefore, the computation complexity of solving the local estimator and the Oracle Huber estimator are about $O(K(pn^{2}+pn+np^{3}))$ and $O(KpN^{2}+(Kp+Kp^{2})N+Kp^{3})$ respectively. As Online Updating algorithm contains $b$ loops and calculating the key features of each subset cost about $O(p^{2}n)$, total computation complexity of Online Updating algorithm is about $O(Kb(pn^{2}+pn+n+p^{3})+bp^{2}n)$. By $N=nb$, we can simplify it into $O(KpNn+(Kp^{2}+Kp+1)N+Kbp^{3})$. Although $K$ is set with a large number(mostly 1000), it can achieve convergence within 10 times of iteration in practice, while total sample $N$ can reach a million or more, thus $K\ll N$. As for $p$, we only consider low-dimension cases where $p\ll N$. To sum up, the dominants part of computation complexity of Oracle and Online Updating ones are $O(KpN^{2})$ and $O(KpNn)$. By thorough comparison, we can find the only difference between these two are $N$ and $n$, where $n\ll N$. In conclusion, theoretically, our Online Updating Huber Regression has lower computation complexity than Oracle Huber Regression.