Online Censoring for Large-Scale Regressions with Application to Streaming Big Data

I Introduction

Nowadays omni-present monitoring sensors, search engines, rating sites, and Internet-friendly portable devices generate massive volumes of typically dynamic data . The task of extracting the most informative, yet low-dimensional structure from high-dimensional datasets is thus of utmost importance. Fast-streaming and large in volume data, motivate well updating analytics rather than re-calculating new ones from scratch, each time a new observation becomes available. Redundancy is an attribute of massive datasets encountered in various applications , and exploiting it judiciously offers an effective means of reducing data processing costs.

In this regard, the notion of optimal design of experiments has been advocated for reducing the number of data required for inference tasks . In recent works, the importance of sequential optimization along with random sampling of Big Data has been highlighted . Specifically for linear regressions, random projection (RP)-based methods have been advocated for reducing the size of large-scale least-squares (LS) problems . As far as online alternatives, the randomized Kaczmarz’s (a.k.a. normalized least-mean-squares (LMS)) algorithm generates a sequence of linear regression estimates from projections onto convex subsets of the data . Sequential optimization includes stochastic approximation, along with recent advances on online learning . Frugal solvers of (possibly sparse) linear regressions are available by estimating regression coefficients based on (severely) quantized data ; see also for decentralized sparse LS solvers.

In this context, the idea here draws on interval censoring to discard “less informative” observations. Censoring emerges naturally in several areas, and batch estimators relying on censored data have been used in econometrics, biometrics, and engineering tasks , including survival analysis , saturated metering , and spectrum sensing . It has recently been employed to select data for distributed estimation of parameters and dynamical processes using resource-constrained wireless sensor networks, thus trading off performance for tractability . These works confirm that estimation accuracy achieved with censored measurements can be comparable to that based on uncensored data. Hence, censoring offers the potential to lower data processing costs, a feature certainly desirable in Big Data applications.

To this end, the present work employs interval censoring for large-scale online regressions. Its key novelty is to sequentially test and update regression estimates using censored data. Two censoring strategies are put forth, each tailored for mitigating different costs. In the first one, stochastic approximation algorithms are developed for sequentially updating the regression coefficients with low-complexity first- or second-order iterations to maximize the likelihood of censored and uncensored observations. This strategy is ideal when the number of observations are to be reduced, in order to lower the cost of storage or transmission to a remote estimation site. Relative to , the contribution here is a novel online scheme that greatly reduces storage requirements without requiring feedback from the estimator to sensors. Error bounds are derived, while simulations demonstrate performance close to estimation error limits.

The second censoring strategy focuses on reducing the complexity of large-scale linear regressions. The proposed methods are also online by design, but may also be readily used to reduce the complexity of solving a batch linear regression problem. The difference with dimensionality-reducing alternatives, such as optimal design of experiments, randomized Kaczmarz’s and RP-based methods, is that the introduced technique reduces complexity in a data-driven manner.

The rest of the paper is as follows. A formal problem description is in Section II, while the two censoring rules are introduced in Section II-A. First- and second-order stochastic approximation maximum-likelihood-based algorithms for censored observations are developed in Section III, along with threshold selection rules for controlled data reduction in Section III-B. Adaptive censoring algorithms for reduced-complexity linear regressions are in Section IV, with corresponding threshold selection rules given in Section IV-C, and robust versions of the algorithms outlined in Section IV-D. The proposed online-censoring and reduced-complexity methods are tested on synthetic as well as real data, and compared with competing alternatives in Section V. Finally, concluding remarks are made in Section VI.

II Problem Statement and Preliminaries

Consider a $p\times 1$ vector of unknown parameters $\boldsymbol{\theta}_{o}$ generating scalar streaming observations

where $\mathbf{x}_{n}$ is the $n$-th row of the $D\times{p}$ regression matrix $\mathbf{X}$, and the noise samples $\upsilon_{n}$ are assumed independently drawn from $\mathcal{N}(0,\sigma^{2})$. The high-level goal is to estimate $\boldsymbol{\theta}_{o}$ in an online fashion, while meeting minimal resource requirements. The term resources here refers to the total number of utilized observations and/or regression rows, as well as the overall computational complexity of the estimation task. Furthermore, the sought data-and complexity-reduction schemes are desired to be data-adaptive, and thus scalable to the size of any given dataset $\{y_{n},\mathbf{x}_{n}\}_{n=1}^{D}$. To meet such requirements, the proposed first- and second-order online estimation algorithms are based on the following two distinct censoring methods.

A generic censoring rule for the data in (1) is given by

where $\ast$ denotes an unknown value when the $n$-th datum has been censored (thus it is unavailable) - a case when we only know that $y_{n}\in\mathcal{C}_{n}$ for some set $\mathcal{C}_{n}$; otherwise, the actual measurement $y_{n}$ is observed. Given $\{z_{n},\mathbf{x}_{n}\}_{n=1}^{D}$, the goal is to estimate $\boldsymbol{\theta}_{o}$. Aiming to reduce the cost of storage and possible transmission, it is prudent to rely on innovation-based interval censoring of $y_{n}$. To this end, define per time $n$ the binary censoring variable $c_{n}=1$ if $y_{n}\in\mathcal{C}_{n}$; and zero otherwise. Each datum is decided to be censored or not using a predictor $\hat{y}_{n}$ formed using a preliminary (e.g., LS) estimate of $\boldsymbol{\theta}_{o}$ as

where $\{\tau_{n}\}_{n=1}^{D}$ are censoring thresholds, and as in (2), $*$ signifies that the exact value of $y_{n}$ is unavailable. The rule (4) censors measurements whose absolute normalized innovation is smaller than $\tau_{n}$; and it is non-adaptive in the sense that censoring depends on $\hat{\boldsymbol{\theta}}_{K}$ that has been derived from a fixed subset of $K$ measurements. Clearly, the selection of $\{\tau_{n}\}_{n=1}^{D}$ affects the proportion of censored data. Given streaming data $\{z_{n},c_{n},\mathbf{x}_{n}\}$, the next section will consider constructing a sequential estimator of $\boldsymbol{\theta}_{o}$ from censored measurements.

The efficiency of NAC in (4) in terms of selecting informative data depends on the initial estimate $\hat{\boldsymbol{\theta}}_{K}$. A data-adaptive alternative is to take into account all censored data $\{\mathbf{x}_{i},z_{i}\}_{i=1}^{n-1}$ available up to time $n$. Predicting data through the most recent estimate $\hat{\boldsymbol{\theta}}_{n-1}$ defines our data-adaptive censoring (AC) rule:

In Section IV, (5) will be combined with first- and second-order iterations to perform joint estimation and censoring online. Implementing the AC rule requires feeding back $\boldsymbol{\theta}_{n-1}$ from the estimator to the censor, a feature that may be undesirable in distributed estimation setups. Nonetheless, in centralized linear regression, AC is well motivated for reducing the problem dimension and computational complexity.

III Online Estimation with NAC

Since noise samples $\{\upsilon_{n}\}_{n=1}^{D}$ in (1) are independent and (4) applies independently over data, $\{z_{n},c_{n}\}_{n=1}^{D}$ are independent too. With $\mathbf{z}_{D}:=[z_{1},\ldots,z_{D}]^{T}$ and $\mathbf{c}_{D}:=[c_{1},\ldots,c_{D}]^{T}$, the joint pdf is $p(\mathbf{z}_{D},\mathbf{c}_{D};\boldsymbol{\theta})=\prod_{n=1}^{D}p(z_{n},c_{n};\boldsymbol{\theta})$ with

since $c_{n}=0$ means no censoring, and thus $z_{n}=y_{n}$ is Gaussian distributed; whereas $c_{n}=1$ implies $|y_{n}-\hat{y}_{n}|\leq{\tau_{n}\sigma}$, that is $\Pr\{c_{n}=1\}=\Pr\{\hat{y}_{n}-\tau_{n}\sigma-\mathbf{x}_{n}^{T}\boldsymbol{\theta}_{0}\leq{v_{n}}\leq{\hat{y}_{n}+\tau_{n}\sigma-\mathbf{x}_{n}^{T}\boldsymbol{\theta}_{0}}\}$, and after recalling that $v_{n}$ is Gaussian

where $z_{n}^{l}(\boldsymbol{\theta}):=-\tau_{n}-\frac{\mathbf{x}_{n}^{T}\boldsymbol{\theta}-\hat{y}_{n}}{\sigma}$ and $z_{n}^{u}(\boldsymbol{\theta}):=\tau_{n}-\frac{\mathbf{x}_{n}^{T}\boldsymbol{\theta}-\hat{y}_{n}}{\sigma}$. Then, the maximum-likelihood estimator (MLE) of $\boldsymbol{\theta}_{o}$ is

If the entire dataset $\{z_{n},c_{n},\mathbf{x}_{n}\}_{n=1}^{D}$ were available, the MLE could be obtained via gradient descent or Newton iterations.

Considering Big Data applications where storage resources are scarce, we resort to a stochastic approximation solution and process censored data sequentially. In particular, when datum $n$ becomes available, the unknown parameter is updated as

The overall scheme is tabulated as Algorithm 1.

Observe that when the $n$-th datum is not censored $(c_{n}=0)$, the second summand in the right-hand side (RHS) of (9) vanishes, and (8) reduces to an ordinary LMS update. When $c_{n}=1$, the first summand disappears, and the update in (8) exploits the fact that the unavailable $y_{n}$ lies in a known interval $(|y_{n}-\mathbf{x}_{n}^{T}\hat{\boldsymbol{\theta}}_{K}|\leq{\tau_{n}\sigma})$, information that would have been ignored by an ordinary LMS algorithm.

Since the SA-MLE is in fact a Robbins-Monroe iteration on the sequence $\{\mathbf{g}(\boldsymbol{\theta})\}_{n=1}^{D}$, it inherits related convergence properties. Specifically, by selecting $\mu_{n}=1/(nM)$ (for an appropriate $M$), the SA-MLE algorithm is asymptotically efficient and Gaussian [21, pg. 197]. Performance guarantees also hold with finite samples. Indeed, with $D$ finite, the regret attained by iterates $\{\boldsymbol{\theta}_{n}\}$ against a vector $\boldsymbol{\theta}$ is defined as

Selecting $\mu$ properly, Algorithm 1 can afford bounded regret as asserted next; see Appendix for the proof.

Suppose $\|\mathbf{x}_{n}\|_{2}\leq\bar{x}$ and $|\beta_{n}(\boldsymbol{\theta})|\leq\bar{\beta}$ for $n=1,\ldots,D$, and let $\boldsymbol{\theta}^{\ast}$ be the minimizer of (7). By choosing $\mu=\|\boldsymbol{\theta}^{\ast}-\hat{\boldsymbol{\theta}}_{K}\|_{2}/(\sqrt{2D}\bar{\beta}\bar{x})$, the regret of the SA-MLE satisfies

Proposition 1 assumes bounded $\mathbf{x}_{n}$’s and noise. Although the latter is not satisfied by e.g., the Gaussian distribution, appropriate bounds ensure that (1) holds with high probability.

If extra complexity can be afforded, one may consider incorporating second-order information in the SA-MLE update to improve its performance. In practice, this is possible by replacing scalar with matrix step-sizes $\mathbf{M}_{n}$. Thus, the first-order stochastic gradient descent (SGD) update in (8) is modified as follows

Due to the rank-one update $\mathbf{M}_{n}=((n-1)/n)\mathbf{M}_{n-1}+(1/n)\gamma_{n-1}(\boldsymbol{\theta}_{n-1})$ $\mathbf{x}_{n-1}\mathbf{x}_{n-1}^{T}$, the matrix step size $\mathbf{C}_{n}:=\mathbf{M}_{n}^{-1}$ can be obtained efficiently using the matrix inversion lemma as

Similar to its first-order counterpart, the algorithm is initialized by the preliminary estimate $\boldsymbol{\theta}_{0}=\hat{\boldsymbol{\theta}}_{K}$, and $\mathbf{C}_{0}=\sigma^{2}(\mathbf{X}_{K}^{T}\mathbf{X}_{K})^{-1}$. The second-order SA-MLE method is summarized as Algorithm 2, while the numerical tests of Section V-A confirm its faster convergence at the cost of $\mathcal{O}(p^{2})$ complexity per update.

III-B Controlling Data Reduction via NAC

To apply the NAC rule of (4) for data reduction at a controllable rate, a relation between thresholds $\{\tau_{n}\}$ and the censoring rate must be derived. Furthermore, prior knowledge of the problem at hand (e.g., observations likely to contain outliers) may dictate a specific pattern of censoring probabilities $\{\pi_{n}^{\ast}\}_{n=1}^{D}$. If $d$ is the number of uncensored data after NAC is applied on a dataset of size $D\geq{d}$, then $(D-d)/D$ is the censoring ratio. Since $\{y_{n}\}$ are generated randomly according to (1), it is clear that $d$ is itself a random variable. The analysis is thus focused on the average censoring ratio

where $\pi_{n}:=\Pr(c_{n}=1)$ is the probability of censoring datum $n$, that as a function of $\tau_{n}$ is given by [cf. (4)]

By the properties of the LSE, $\hat{\boldsymbol{\theta}}_{K}\sim{\mathcal{N}(\boldsymbol{\theta}_{o},\sigma^{2}(\mathbf{X}_{K}^{T}\mathbf{X}_{K})^{-1})}$, it follows that

Thus, the censoring probabilities in (III-B) simplify to

Solving (16) for $\tau_{n}$, one arrives for a given ${\pi}^{\star}_{n}=\pi_{n}(\tau_{n}^{\star})$ at

Hence, for a prescribed $\bar{c}$, one can select a desired censoring probability pattern $\{\pi^{\star}_{n}\}_{n=1}^{D}$ to satisfy (14), and corresponding $\{\tau_{n}^{\star}\}_{n=1}^{D}$ in accordance with (17).

As expected, due to the normalization by $\sigma$ in (4), ${\pi}$ does not depend on $\sigma$. Interestingly, it does not depend on $\mathbf{R}_{x}$ either. Having expressed $\pi$ as a function of $\tau$, the latter can be tuned to achieve the desirable data reduction. Following the law of large numbers and given parameters $p$ and $K$, to achieve an average censoring ratio of $\bar{c}={\pi}^{\star}=(D-d)/D$, the threshold can be set to

Figure 1 depicts $\pi$ as a function of $\tau$ for $p=100$ and $K=200$. Function (III-B) is compared with the simulation-based estimate of $\pi_{n}$ using 100 Monte Carlo runs, confirming that (III-B) offers a reliable approximation of ${\pi}$, which improves as $p$ grows. However, for the approximation $(\mathbf{X}_{K}^{T}\mathbf{X}_{K})^{-1}\approx\mathbf{R}_{x}^{-1}/{K}$ to be accurate, $K$ should be large too. Figure 1 shows the probability of censoring for varying $K$ with fixed $p=100$ and $\tau=1$. Approximation (III-B) yields a reliable value for $\pi$ for as few as $K\approx 200$ preliminary data.

IV Big Data Streaming Regression with AC

The NAC-based algorithms of Section III emerge in a wide range of applications for which censoring occurs naturally as part of the data acquisition process; see e.g., the Tobit model in economics , and survival data analytics in . Apart from these applications where data are inherently censored, our idea is to employ censoring deliberately for data reduction. Leveraging NAC for data reduction decouples censoring from estimation, and thus eliminates the need for obtaining further information. However, one intuitively expects improved performance with a joint censoring-estimation design.

In this context, first- and second-order sequential algorithms will be developed in this section for the AC in (5). Instead of $\hat{\boldsymbol{\theta}}_{K}$, AC is performed using the latest estimate of $\boldsymbol{\theta}$. Apart from being effective in handling streaming data, AC can markedly lower the complexity of a batch LS problem. Section IV-A introduces an AC-based LMS algorithm for large-scale streaming regressions, while Section IV-B puts forth an AC-based recursive least-squares (RLS) algorithm as a viable alternative to random projections and sampling.

A first-order AC-based algorithm is presented here, inspired by the celebrated LMS algorithm. Originally developed for adaptive filtering, LMS is well motivated for low-complexity online estimation of (possibly slow-varying) parameters. Given $(y_{n},\mathbf{x}_{n})$, LMS entails the simple update

for a given $\tau_{n}>0$. For the sake of analysis, a common threshold will be adopted; that is, $\tau_{n}=\tau$ $\forall n$. The truncated cost can be also expressed as $f_{n}^{(\tau)}(\boldsymbol{\theta})=\max\{0,(e_{n}^{2}(\boldsymbol{\theta})-\tau^{2}\sigma^{2})/2\}$. Being the pointwise maximum of two convex functions, $f_{n}^{(\tau)}(\boldsymbol{\theta})$ is convex, yet not everywhere differentiable. From standard rules of subdifferential calculus, its subgradient is

An SGD iteration for the instantaneous cost in (23) with $\tau_{n}=\tau$, performs the following AC-LMS update per datum $n$

where $\mu>0$ can be either constant for tracking a time-varying parameter, or, diminishing over time for estimating a time-invariant $\boldsymbol{\theta}_{o}$. Different from SA-MLE, the AC-LMS does not update $\boldsymbol{\theta}$ if datum $n$ is censored. The intuition is that if $y_{n}$ can be closely predicted by $\hat{y}_{n}:=\mathbf{x}_{n}^{T}\boldsymbol{\theta}_{n-1}$, then $(y_{n},\mathbf{x}_{n})$ can be censored (small innovation is indeed ‘not much informative’). Extracting interval information through a likelihood function as in Algorithm 1 appears to be challenging here. This is because unlike NAC, the AC data $\{z_{n}\}_{n=1}^{D}$ are dependent across time.

Interestingly, upon invoking the “independent-data assumption” of SA , following the same steps as in Section III, and substituting $\hat{\boldsymbol{\theta}}_{K}=\boldsymbol{\theta}_{n-1}$ into (9), the interval information term is eliminated. This is a strong indication that interval information from censored observations may be completely ignored without the risk of introducing bias. Indeed, one of the implications of the ensuing Proposition 2 is that the AC-LMS is asymptotically unbiased. Essentially, in AC-LMS as well as in the AC-RLS to be introduced later, both $\mathbf{x}_{n}$ and $y_{n}$ are censored – an important feature effecting further data reduction and lowering computational complexity of the proposed AC algorithms. The mean-square error (MSE) performance of AC-LMS is established in the next proposition proved in the Appendix.

Proposition 2 asserts that AC-LMS achieves a bounded MSE. It also links MSE with the AC threshold $\tau$ that can be used to adjust the censoring probability. Closer inspection reveals that the MSE bound decreases with $\tau$. In par with intuition, lowering $\tau$ allows the estimator to access more data, thus enhancing estimation performance at the price of increasing the data volume processed.

IV-B AC-RLS

A second-order AC algorithm is introduced here for the purpose of sequential estimation and dimensionality reduction. It is closely related to the RLS algorithm, which per time $n$ implements the updates; see e.g.,

where $\mathbf{C}_{n}$ is the sample estimate for $\mathbf{R}_{x}^{-1}$ and is typically initialized to $\mathbf{C}_{0}=\epsilon\mathbf{I}$, for some small positive $\epsilon$, e.g., . The RLS estimate at time $n$ can be also obtained as

To obtain a second-order counterpart of AC-LMS, we replace the quadratic instantaneous cost of RLS with the truncated quadratic in (23). The matrix step-size is further surrogated by

Applying the matrix inversion lemma to find $\mathbf{M}_{n}^{-1}$ yields the next AC-RLS updates

where $c_{n}$ is decided by (5). For $c_{n}=1$, the parameter vector is not updated, while costly updates of $\mathbf{C}_{n}$ are also avoided. In addition, different from the iterative expectation-maximization algorithm in , AC-RLS skips completely covariance updates. Its performance is characterized by the following proposition shown in the Appendix.

As corroborated by Proposition 3, the AC-RLS estimates are guaranteed to converge to $\boldsymbol{\theta}_{o}$ for any choice of $\tau$. Overall, the novel AC-RLS algorithm offers a computationally-efficient and accurate means of solving large-scale LS problems encountered with Big Data applications.

At this point, it is useful to contrast and compare AC-RLS with RP and random sampling methods that have been advocated as fast LS solvers . In practice, RP-based schemes first premultiply data $(\mathbf{y},\mathbf{X})$ with a random matrix $\mathbf{R}=\mathbf{HD}$, where $\mathbf{H}$ is a $D\times D$ Hadamard matrix and $\mathbf{D}$ is a diagonal matrix whose diagonal entries take values $\{-1/\sqrt{D},+1/\sqrt{D}\}$ equiprobably. Intuitively, $\mathbf{R}$ renders all rows of “comparable importance” (quantified by the leverage scores ), so that the ensuing random matrix $\mathbf{S}_{d}$ exhibits no preference in selecting uniformly a subset of $d$ rows. Then, the reduced-size LS problem can be solved as $\check{\boldsymbol{\theta}}_{d}=\arg\min_{\boldsymbol{\theta}}{\|{\mathbf{S}_{d}\mathbf{HD}(\mathbf{y}-\mathbf{X}\boldsymbol{\theta})}\|_{2}^{2}}$. For a general preconditioning matrix $\mathbf{HD}$, computing the products $\mathbf{HDy}$ and $\mathbf{HDX}$ requires a prohibitive number of $\mathcal{O}(D^{2}p)$ computations. This is mitigated by the fact that $\mathbf{H}$ has binary $\{+1,-1\}$ entries and thus multiplications can be implemented as simple sign flips. Overall, the RP method reduces the computational complexity of the LS problem from $\mathcal{O}(Dp^{2})$ to $\scriptstyle{\mathcal{O}}$$(Dp^{2})$ operations.

By setting $\tau=Q^{-1}(d/(2D))$, our AC-RLS Algorithm 3 achieves an average reduction ratio $d/D$ by scanning the observations, and selecting only the most informative ones. The same data ratio can be achieved more accurately by choosing a sequence of data-adaptive thresholds $\{\tau_{n}\}_{n=1}^{D}$, as described in the next subsection. As will be seen in Section V-C, AC-RLS achieves significantly lower estimation error compared to RP-based solvers. Intuitively, this is because unlike RPs that are based solely on $\mathbf{X}$ and are thus observation-agnostic, AC extracts the most informative in terms of innovation subset of rows for a given problem instance $(\mathbf{y},\mathbf{X})$.

Regarding the complexity of AC-RLS, if the pair $(y_{n},\mathbf{x}_{n})$ is not censored, the cost of updating $\boldsymbol{\theta_{n}}$ and $\mathbf{C}_{n}$ is $\mathcal{O}(p^{2})$ multiplications. For a censored datum, there is no such cost. Thus, for $d$ uncensored data the overall computational complexity is $\mathcal{O}(dp^{2})$. Furthermore, evaluation of the absolute normalized innovation requires $\mathcal{O}(p)$ multiplications per iteration. Since this operation takes place at each of the $D$ iterations, there are $\mathcal{O}(Dp)$ computations to be accounted for. Overall, AC-RLS reduces the complexity of LS from $\mathcal{O}(Dp^{2})$ to $\mathcal{O}(dp^{2})+\mathcal{O}(Dp)$. Evidently, the complexity reduction is more prominent for larger model dimension $p$. For $p\gg{1}$, the second term may be neglected, yielding an $\mathcal{O}(dp^{2})$ complexity for AC-RLS.

The novel AC-LMS and AC-RLS algorithms bear structural similarities to sequential set-membership (SM)-based estimation . However, the model assumptions and objectives of the two are different. SM assumes that the noise distribution in (1) has bounded support, which implies that $\boldsymbol{\theta}_{o}$ belongs to a closed set. This set is sequentially identified by algorithms interpreted geometrically, while certain observations may be deemed redundant and thus discarded by the SM estimator. In our Big Data setup, an SA approach is developed to deliberately skip updates of low importance for reducing complexity regardless of the noise pdf.

IV-C Controlling Data Reduction via AC

A clear distinction between NAC and AC is that the latter depends on the estimation algorithm used. As a result, threshold design rules are estimation-driven rather than universal. In this section, threshold selection strategies are proposed for AC-RLS. Recall the average reduction ratio $\bar{c}$ in (14), and let $\boldsymbol{\zeta}_{n}:=(\boldsymbol{\theta}_{o}-\boldsymbol{\theta}_{n})/\sigma\sim{\mathcal{N}(\mathbf{0},\mathbf{K}_{n})}$ denote the normalized error at the $n-$th iteration. Similar to (14)–(III-B), it holds that

For $n\gg{p}$, estimates $\boldsymbol{\theta}_{n}$ are sufficiently close to $\boldsymbol{\theta}_{o}$ and thus $\mathbf{K}_{n}\approx\mathbf{0}$. Then, the data-agnostic $\tau_{n}\approx Q^{-1}(\frac{1-\pi_{n}}{2})$ attains an average censoring probability $\bar{\pi}$, while its asymptotic properties have been studied in . For finite data, this simple rule leads to under-censoring by ignoring appreciable values of $\mathbf{K}_{n}$, which can increase computational complexity considerably. This consideration motivates well the data-adaptive threshold selection rules designed next.

AC-RLS updates can be seen as ordinary RLS updates on the subsequence of uncensored data. After ignoring the transient error due to initialization, it holds that $\mathbf{K}_{n}\approx\left[\sum_{i=1}^{n}(1-c_{i})\mathbf{x}_{i}\mathbf{x}_{i}^{T}\right]^{-1}$. The term $\mathbf{x}_{n}^{T}\mathbf{K}_{n-1}\mathbf{x}_{n}$ is encountered as $\mathbf{x}_{n}^{T}\mathbf{C}_{n-1}\mathbf{x}_{n}/n$ in the updates of Alg. 3, but it is not computed for censored measurements. Nonetheless, $\mathbf{x}_{n}^{T}\mathbf{C}_{n-1}\mathbf{x}_{n}/n$ can be obtained at the cost of $p(p+1)$ multiplications per censored datum. Then, the exact censoring probability at AC-RLS iteration $n$ can be tuned to a prescribed $\pi^{\star}_{n}$ by selecting

Given $\{\pi^{\star}_{n}\}_{n=1}^{D}$ satisfying (14), an average censoring ratio of $(D-d)/D$ is thus achieved in a controlled fashion.

To attain $\pi^{\star}_{n}$, the threshold per datum $n$ is selected as

It is well known that for large $n$, the RLS error covariance matrix $\mathbf{K}_{n}$ converges to $\frac{\sigma^{2}}{n}\mathbf{R}_{x}^{-1}$. Specifying $\{\pi^{\star}_{n}\}_{n=1}^{D}$ is equivalent to selecting an average number of $\sum_{i=1}^{n}(1-\pi^{\star}_{i})$ RLS iterations until time $n$. Thus, the AC-RLS with controlled selection probabilities yields an error covariance matrix $\mathbf{K}_{n}\approx\left(\sum_{i=1}^{n}(1-\pi^{\star}_{i})\right)^{-1}\sigma^{2}\mathbf{R}_{\mathbf{x}}^{-1}$. Combined with (30), the latter leads to

Plugging $\sigma_{e_{n}}$ into (30) yields the simple threshold selection

Unlike (29), where thresholds are decided online at an additional computational cost, (31) offers an off-line threshold design strategy for AC-RLS. Based on (31), to achieve $\bar{c}={\pi}^{\star}=(D-d)/D$, thresholds are chosen as

which attains a constant $\pi^{\ast}$ across iterations.

IV-D Robust AC-LMS and AC-RLS

AC-LMS and AC-RLS were designed to adaptively select data with relatively large innovation. This is reasonable provided that (1) contains no outliers whose extreme values may give rise to large innovations too, and thus be mistaken for informative data. Our idea to gain robustness against outliers is to adopt the modified AC rule

Similar to (5), a nominal censoring variable $c_{n}$ is activated here too for observations with absolute normalized innovation less than $\tau$. To reveal possible outliers, a second censoring variable $c_{n}^{o}$ is triggered when the absolute normalized innovation exceeds threshold $\tau_{o}>\tau.$

Having separated data-censoring from outlier identification in (33), it becomes possible to robustify AC-LMS and AC-RLS against outliers. Towards this end, one approach is to completely ignore $y_{n}$ when $c_{n}^{o}=1$. Alternatively, the instantaneous cost function in (23) can be modified to a truncated Huber loss (cf. )

Applying the first-order SGD iteration on the cost $f^{o}(e_{n})$, yields the robust (r) AC-LMS iteration

Similarly, the second-order SGD yields the rAC-RLS

Observe that when $c_{n}^{o}=1$, only $\boldsymbol{\theta}_{n}$ is updated, and the computationally costly update of (35b) is avoided.

V Numerical Tests

To further evaluate the efficacy of the proposed methods, additional simulations were run for different levels of censoring by adjusting $\tau$. Plotted in Figs. 3 and 3 are the MSE curves of the first- and second-order SA-MLE respectively, for different values of $\tau$. Notice that censoring up to $50\%$ of the data (green solid curve) incurs negligible estimation error compared to the full-data case (blue solid curve). In fact, even when operating on data reduced by $95\%$ (red dashed curve) the proposed algorithms yield reliable online estimates.

V-B AC-LMS comparison with Randomized Kaczmarz

V-C AC-RLS

The AC-RLS algorithm developed in Section IV-B was tested on synthetic data. Specifically, the AC-RLS is treated here as an iterative method that sweeps once through the entire dataset, even though more sweeps can be performed at the cost of additional runtime. Its performance in terms of relative MSE was compared with the Hadamard (HD) preconditioned randomized LS solver, while plotted as a function of the compression ratio $d/D$. Parallel to the two methods, a uniform sampling randomized LSE was run as a simple benchmark. Measurements were generated according to (1) with $p=300$, $D=10,000$, and $v_{n}\sim{\mathcal{N}(0,9)}$. Regarding the data distribution, three different scenario’s were examined. In Figure 5, $\mathbf{x}_{n}$’s were generated according to a heavy tailed multivariate $t-$distribution with one degree of freedom, and covariance matrix with $(i,j)$-th entry $\boldsymbol{\Sigma}_{i,j}=2\times{0.5}^{|i-j|}$. Such a data distribution yields matrices $\mathbf{X}$ with highly non-uniform leverage scores, thus imitating the effect of a subset of highly “important” observations randomly scattered in the dataset. In such cases, uniform sampling without preconditioning performs poorly since many of those informative measurements are missed. As seen in the plot, preconditioning significantly improves performance, by incorporating “important” information through random projections. Further improvement is effected by our data-driven AC-RLS through adaptively selecting the most informative measurements and ignoring the rest, without overhead in complexity.

The experiment was repeated (Fig. 5) for $\mathbf{x}_{n}$ generated from a multivariate $t-$distribution with 3 degrees of freedom, and $\boldsymbol{\Sigma}$ as before. Leverage scores for this dataset are moderately non-uniform, thus inducing more redundancy and resulting in lower performance for all algorithms, while closing the “gap” between preconditioned and non-preconditioned random sampling. Again, the proposed AC-RLS performs significantly better in estimating the unknown parameters for the entire range of data size reduction.

Finally, Fig. 5 depicts related performance for Gaussian $\mathbf{x}_{n}\sim{\mathcal{N}(\mathbf{0},\boldsymbol{\Sigma})}$. Compared to the previous cases, normally distributed rows yield a highly redundant set of measurements with $\mathbf{X}$ having almost uniform leverage scores. As seen in the plots, preconditioning offers no improvement in random sampling for this type data, whereas the AC-RLS succeeds in extracting more information on the unknown $\boldsymbol{\theta}$.

To further assess efficacy of the AC-RLS algorithm, real data tests were performed. The Protein Tertiary Structure dataset from the UCI Machine Learning Repository was tested. In this linear regression dataset, $p=9$ attributes of proteins are used to predict a value related to protein structure. A total of $D=45,730$ observations are included. Since the true $\boldsymbol{\theta}_{o}$ is unknown, it is estimated by solving LS on the entire dataset. Subsequently, the noise variance is also estimated via sample averaging as $\sigma^{2}=(1/D)\sum_{n=1}^{D}{(y_{n}-\mathbf{x}_{n}^{T}\boldsymbol{\theta}_{o})^{2}}$. Figure 6 depicts relative squared-error (RSE) with respect to the data reduction ratio $d/D$. The RSE curve for the HD-preconditioned LS corresponds to the average RSE across 50 runs, while the size of the vertical bars is proportional to its standard deviation. Different from RP-based methods, the RSE for AC-RLS does not entail standard deviation bars, because for a given initialization and data order, the output of the algorithm is deterministic. It can be observed that for $d/D\geq{0.25}$ the AC-RLS outperforms RPs in terms of estimating $\boldsymbol{\theta}$, while for very small $d/D$, RPs yield a lower average RSE, at the cost however of very high error uncertainty (variance).

V-D Robust AC-RLS

To test rAC-LMS and rAC-RLS of Section IV-D, datasets were generated with $D=10,000$, $p=30$ and $\mathbf{x}_{n}\sim{\mathcal{N}(\mathbf{0},\boldsymbol{\Sigma})}$, where $\boldsymbol{\Sigma}_{i,j}=2\times{0.5}^{|i-j|}$; noise was i.i.d. Gaussian $v_{n}\sim{\mathcal{N}(0,9)}$; meanwhile measurements $y_{n}$ were generated according to (1) with random and sporadic outlier spikes $\{o_{n}\}_{n=1}^{D}$. Specifically, we generated $o_{n}=\alpha_{n}\beta_{n}$, where $\alpha_{n}\sim{\textrm{Bernoulli}}(0.05)$, and $\beta_{n}\sim{\mathcal{N}(0,25\times{9})}$, thus resulting in approximately $5\%$ of the data effectively being outliers. Similar to previous experiments, our novel algorithms were run once through the set selecting $d$ out of $D$ data to update $\boldsymbol{\theta}_{n}$. Plotted in Fig. 7 is the RSE averaged across 100 runs as a function of $d/D$ for the HD-preconditioned LS, the plain AC-RLS, and the rAC-RLS with a Huber-like instantaneous cost. As expected, the performance of AC-RLS is severely undermined especially when tuned for very small $d/D$, exhibiting higher error than the RP-based LS. However, our rAC-RLS algorithm offers superior performance across the entire range of $d/D$ values.

VI Concluding Remarks

We developed online algorithms for large-scale LS linear regressions that rely on censoring for data-driven dimensionality reduction of streaming Big Data. First, a non-adaptive censoring setting was considered for applications where observations are censored – possibly naturally – separately and prior to estimation. Computationally efficient first- and second-order online algorithms were derived to estimate the unknown parameters, relying on stochastic approximation of the log-likelihood of the censored data. Performance was bounded analytically, while simulations demonstrated that the second-order method performs close to the CRLB.

Furthermore, online data reduction occurring parallel to estimation was also explored. For this scenario, censoring is performed deliberately and adaptively based on estimates provided by first- and second-order algorithms. Robust versions were also developed for estimation in the presence of outliers. Studied under the scope of stochastic approximation, the proposed algorithms were shown to enjoy guaranteed MSE performance. Moreover, the resulting recursive methods were advocated as low-complexity recursive solvers of large LS problems. Experiments run on synthetic and real datasets corroborated that the novel AC-LMS and AC-RLS algorithms outperformed competing randomized algorithms.

Our future research agenda includes approaches to nonlinear (e.g., kernel-based) parametric and nonparametric large-scale regressions, along with estimation of dynamical (e.g., state-space) processes using adaptively censored measurements.

where $\{\boldsymbol{\theta}_{n}\}_{n=1}^{D}$ is any sequence of estimates produced by the SA-MLE. By choosing $\mu=\|\boldsymbol{\theta}^{\ast}-\hat{\boldsymbol{\theta}}_{K}\|_{2}/(\sqrt{2D}\bar{\beta}\bar{x})$, the aforementioned bound leads to Proposition 1. ∎

Under a3), there exists a constant $\alpha>0$ such that $\nabla^{2}F(\boldsymbol{\theta})\succeq{\alpha\mathbf{I}}$ $\forall\boldsymbol{\theta}$. Interchanging differentiation with expectation yields

It can be verified that the function $g(z):=Q(\tau+z)+Q(\tau-z)$ is minimized for $z=0$ when $\tau>0$. To see this, observe that its derivative $g^{\prime}(z)=-\phi(\tau+z)+\phi(\tau-z)$ vanishes when $|\tau+z|=|\tau-z|$. Therefore, $g(z)\geq g(0)=2Q(\tau)$ for all $z$; and hence,

for all $\mathbf{x}$ and $\boldsymbol{\theta}$. The latter implies

showing that $F(\boldsymbol{\theta})$ is $\alpha-$strongly convex with $\alpha=2Q(\tau)\lambda_{\min}(\mathbf{R}_{x})$. As expected, $\alpha$ reduces for increasing $\tau$.

It can be verified that since the cross-terms in (36) can be bounded from below and above as

The last expression reveals that the average distance between gradients can be decomposed into two terms. The first term can be bounded using the fourth-order moment. The second term appears due to data censoring and clearly depends on $\tau$, while it is assumed bounded as

Finally, the expected norm of the gradient at $\boldsymbol{\theta}=\boldsymbol{\theta}_{o}$ is bounded and equal to

Since $\mathbf{C}_{n}$ converges monotonically to $\mathbf{C}_{\infty}$, there exists $k>0$ such that for all $n>k$

References