For a matrix $\bm{X}\in\mathbb{R}^{n\times p}$, we let $\bm{X}_{i\cdot},\bm{X}_{\cdot j}$ and $\bm{X}_{ij}$ denote the $i$-th row, $j$-th column and $(i,j)$-element of matrix $\bm{X}.$ $\bm{X}_{i,-j}$ is a sub-vector of $\bm{X}_{i\cdot}$ with the $j$-th element deleted and $\bm{X}_{-i,-j}$ is a sub-matrix of $\bm{X}$ with the $i$-th rows and the $j$-th column deleted while other elements remain unchanged. For a vector $\bm{x}\in\mathbb{R}^{m}$, we define its $\ell_{q}$-norm as $\lVert\bm{x}\lVert_{q}=(\sum_{i=1}^{p}x_{i}^{q})^{\frac{1}{q}}$. For a sequence of random variables $\{\xi_{n}\}_{n\in\mathbb{N}}$ and a corresponding sequence of constants $\{a_{n}\}_{n\in\mathbb{N}}$. We say that $\xi_{n}=\mathcal{O}_{p}(a_{n})$ if for any $\epsilon>0$, there exist two finite numbers $M,N>0$ such that $P(|\xi_{n}/a_{n}|>M)<\epsilon$ for any $n>N.$ Generally speaking, $\xi_{n}=\mathcal{O}_{p}(a_{n})$ denotes $\xi_{n}/a_{n}$ is stochastically bounded. $\xi_{n}=o_{p}(a_{n})$ means that $\xi_{n}/a_{n}$ converges to zero in probability. With the consideration of the streaming data, we use $\bm{X}^{(j)}$ and $\bm{Y}^{(j)}$ to stand for $\bm{X}$ and $\bm{y}$, namely explanatory variable and explained variables, arriving in $j$-th batch respectively. In addition, $\bm{X}^{(j)}_{\star}$ and $\bm{Y}^{(j)}_{\star}$ (with star index) are the cumulative variables of $\bm{X}^{(j)}$ and $\bm{Y}^{(j)}$. For example, $\bm{X}^{(j)}_{\star}=((\bm{X}^{(1)})^{\top},\ldots,(\bm{X}^{(j)})^{\top})^{\top}.$

We first consider the online lasso estimator through the gradient descent method. Define the score function $\bm{u}(y_{i};\bm{x}_{i},\bm{\beta}):=\partial d(y_{i};g(\bm{x}_{i}^{\top}\bm{\beta}))/\partial\bm{\beta}$. The negative gradient of the first term in (1) $\mathcal{L}(\bm{\beta}):=\sum_{i\in\mathcal{D}_{b}^{\star}}d(y_{i};g(\bm{x}_{i}^{\top}\bm{\beta}))$ is

Let $U^{(j)}(\bm{\beta}):=\sum_{i\in\mathcal{D}_{j}}\bm{u}(y_{i};\bm{x}_{i},\bm{\beta})$ be the score function pertaining to the data batch $\mathcal{D}_{j}$. Then $\bm{U}^{(b)}(\bm{\beta})$ in (2) can be rewritten as a linear aggregation form:

To pursue an online estimator upon the arrival of $\mathcal{D}_{b}$, a key step is to update ${\bm{U}}^{(b-1)}(\bm{\beta})$ to $\bm{U}^{(b)}(\bm{\beta})$ without re-accessing the whole historical raw data $\mathcal{D}_{b-1}^{\star}$.

To illustrate the idea, we consider a simple case with two data batches, i.e., $\mathcal{D}_{2}$ arrives subsequent to $\mathcal{D}_{1}$. The lasso estimator based on the first batch data is denoted by $\widehat{\bm{\beta}}^{(1)}$, which is the offline estimator $\bar{\bm{\beta}}^{(1)}$ according to (1). To avoid using individual-level raw data in $\mathcal{D}_{1}$, we approximate ${\bm{U}}^{(1)}(\bm{\beta})$ through a first-order Taylor expansion at $\widehat{\bm{\beta}}^{(1)}$. Note that $\widehat{\bm{\beta}}^{(1)}$ satisfies $\lVert\bm{U}^{(1)}(\widehat{\bm{\beta}}^{(1)})\lVert_{\infty}=\mathcal{O}_{p}(\lambda_{1}N_{1})$, which implies

where $\bm{J}^{(1)}(\bm{\beta})=-\partial U^{(1)}(\bm{\beta})/\partial\bm{\beta}$ denotes the information matrix. Next, we propose to approximate $\bm{U}^{(2)}(\bm{\beta})$ by

Apparently, calculating $\widetilde{\bm{U}}^{(2)}(\bm{\beta})$ through (3) only requires access to the summary statistics $\{\widehat{\bm{\beta}}^{(1)},\bm{J}^{(1)}(\widehat{\bm{\beta}}^{(1)})\}$ rather than the full dataset $\mathcal{D}_{1}$.

Now we are ready to generalize the above approximation to an arbitrary data batch $\mathcal{D}_{b}$. Let $\widetilde{\bm{J}}^{(b-1)}=\sum_{j=1}^{b-1}\bm{J}^{(j)}(\widehat{\bm{\beta}}^{(j)})$ denote the aggregated information matrix. The approximation procedure in (3) becomes

The aggregated gradient $\widetilde{\bm{U}}^{(b)}(\bm{\beta})$ depends only on $\{\widehat{\bm{\beta}}^{(b-1)},\widetilde{\bm{J}}^{(b-1)}\}$. Hence, $\widehat{\bm{\beta}}^{(b)}$ can be computed through the following algorithm.

• •

  Step 1: update $\widehat{\bm{\beta}}^{(b)}$ through

  $$\widehat{\bm{\beta}}^{(b)}\leftarrow\widehat{\bm{\beta}}^{(b)}+\frac{\eta}{2N_{b}}\widetilde{\bm{U}}^{(b)}(\widehat{\bm{\beta}}^{(b)}),$$

  where $\eta$ is the learning rate in the gradient descent;

• •

  Step 2: for $r=1,\dots,p$, apply the soft-thresholding operator $\mathcal{S}(\widehat{\beta}_{r}^{(b)};\eta\lambda_{b})$ to the $r$-th component of $\widehat{\bm{\beta}}^{(b)}$ obtained in Step 1:

  $$\mathcal{S}(\widehat{\beta}^{(b)}_{r};\eta\lambda_{b})=\begin{cases}\widehat{\beta}^{(b)}_{r}+\eta\lambda_{b},\quad&\text{if}\ \widehat{\beta}^{(b)}_{r}<-\eta\lambda_{b}\\ 0,\quad&\text{if}\ |\widehat{\beta}^{(b)}_{r}|\leq\eta\lambda_{b}\\ \widehat{\beta}^{(b)}_{r}-\eta\lambda_{b},\quad&\text{if}\ \widehat{\beta}^{(b)}_{r}>\eta\lambda_{b}\end{cases}.$$

Note that $\lambda_{b}$ is the regularization parameter that is chosen adaptively for step $b$. More details on the adaptive tuning procedure is presented in Section 2.3.

The above two steps are carried out iteratively till convergence to obtain the online lasso estimator $\bm{\widehat{\beta}}^{(b)}$. In the implementation, we set the stopping criterion to be $\|{\eta}\times\widetilde{\bm{U}}^{(b)}(\widehat{\bm{\beta}}^{(b)})/{2N_{b}}\|_{2}\leq 10^{-6}$. In summary, our proposed online estimator $\widehat{\bm{\beta}}^{(b)}$ can be defined as

In contrast to the standard offline lasso estimator in (1), our proposed online estimator in (4) depends on the data only through the summary statistics $\{\widehat{\bm{\beta}}^{(b-1)},\widetilde{\bm{J}}^{(b-1)}\}$ instead of $\mathcal{D}_{b-1}^{\star}$.

We now proceed to study the online statistical inference and construct confidence intervals for the $r$-th component of the regression parameter vector, $r=1,\ldots,p$. We first define a low-dimensional projection

where $N_{b}$ and $\lambda_{b}$ are the same as in (4). Letting $\widehat{W}{{}^{(j)}}\in\mathbb{R}^{n_{j}\times n_{j}}$ be the diagonal matrix with diagonal elements $\sqrt{g^{\prime}(\bm{X}^{(j)}\bm{\widehat{\beta}}^{(j)})}$ and $\widehat{\bm{X}}^{{(j)}}$ be the weighted design matrix $\widehat{W}{{}^{(j)}}\bm{X}^{(j)}$, then (5) can be recast into

With this view, $\bm{\widehat{\gamma}}_{r}^{(b)}$ can be computed in a similar fashion as the online lasso estimator. Specifically, (6) has the same form as (1) if we choose the deviance function $d(\cdot;\cdot)$ as the square error. Then, it is straightforward to propose the online algorithm to find $\bm{\widehat{\gamma}}_{r}^{(b)}$ by following the procedure in Section 2.1. At this time, the summary statistics is $(\widetilde{\bm{J}}^{(b)}_{r,-r},\widetilde{\bm{J}}_{-r,-r}^{(b)}),$ which has been already stored in $\widetilde{\bm{J}}^{(b)}.$ As a result, we find a solution of (5). After obtaining $\widehat{\bm{\gamma}}_{r}^{(b)}$ in (5), we compute

Denote $\widetilde{\bm{\gamma}}^{(j)}_{r}=(\widehat{\bm{\gamma}}^{(j)}_{r,1},\ldots,-1,\ldots,\widehat{\bm{\gamma}}^{(j)}_{r,p})^{\top}\in\mathbb{R}^{p}$, whose $r$-th element is $-1$. Then, upon the arrival of the batch data $\mathcal{D}_{b}$, we define our proposed online debiased lasso estimator as

It is worth noting that the debiased lasso estimator involves the initial lasso estimator defined in (4), as well as two additional terms: a debiasing term and an online error correction term. van de Geer et al. , (2014) studied the offline version of debiased lasso for GLM. Our debiased term could be viewed as an online generalization of the offline counterpart in van de Geer et al. , (2014). However, they are rather different in the sense that, the debiased term alone will not suffice to ensure the asymptotic normality. As $\widetilde{\bm{J}}^{(b)}$ is used to approximate the information matrix $\bm{J}(\bm{\beta}^{0})$, the approximation error accumulates even though each ${\widehat{\bm{\beta}}}^{(j)},j=1,\ldots,b,$ is consistent for $\bm{\beta}^{0}$. The correction term in (2.2) is used to eliminate the approximation error arising from the online update.

Meanwhile, the proposed debiased lasso estimator with the online error correction term aligns with the online learning framework, as (2.2) only requires the following summary statistics rather than the entire dataset $\mathcal{D}_{b}^{\star}$:

which keep the same size when new data arrive, and can be easily updated.

The asymptotic normality of the online debiased lasso and the consistency of two lasso-typed estimators in (4) and (5) are established in Section 3. For variance estimation, let

The estimated standard error $\widehat{\sigma}^{(b)}_{r}:={\sqrt{v^{(b)}}}/{\widehat{\tau}_{r}^{(b)}}$ can also be updated online accordingly.

In an offline setting, the regularization parameter $\lambda$ is typically determined by cross-validation where the entire dataset is split into training and test sets multiple times. However, since the full dataset is not accessible in an online setting, such procedure is not feasible. To align with the nature of streaming datasets, we use the “rolling-original-recalibration” procedure with mean squared prediction error (MSE) being the cross-validation criterion (Han et al. ,, 2021). Specifically, at time point $b$, cumulative dataset up to time point $b-1$ serves as the training set while the new data batch $\mathcal{D}_{b}$ is the test set. It is worth noting that instead of re-accessing raw data $\{\mathcal{D}_{1},\dots,\mathcal{D}_{b-1}\}$, we plug in $\widehat{\bm{\beta}}^{(b-1)}$ to evaluate the prediction error for a sequence of $\lambda$ in a candidate set $T_{\lambda}$:

and choose $\lambda$ such that $\lambda_{b}:=\underset{\lambda\in T_{\lambda}}{\arg\min}\ PE_{b}(\lambda)$. The initial $\lambda_{1}$ is selected by the classical offline cross-validation.

We summarize the results in Sections 2.1-2.3 in Figure 1 and Algorithm 1. It consists of two main blocks: one is online lasso estimation and the other is online low-dimensional projection. Outputs from both blocks are used to compute the online debiased lasso estimator as well as the construction of confidence intervals in real-time. In particular, when a new data batch $\mathcal{D}_{b}$ arrives, it is first sent to the online lasso estimation block, where the summary statistics $\left\{\widehat{\bm{\beta}}^{(b-1)},\widetilde{\bm{J}}^{(b-1)}\right\}$ are used to compute $\widetilde{\bm{U}}^{(b)}$. Then we use gradient descent to update the lasso estimator $\widehat{\bm{\beta}}^{(b-1)}$ to $\widehat{\bm{\beta}}^{(b)}$ at a sequence of tuning parameters values without retrieving the whole dataset. At the same time, regarding the cumulative dataset that produces the old lasso estimate $\widehat{\bm{\beta}}^{(b-1)}$ as training set and the newly arrived $\mathcal{D}_{b}$ as testing set, we can choose the tuning parameter $\lambda_{b}$ that gives the smallest prediction error. Now, the selected $\lambda_{b}$ and sub-matrices of $\widetilde{\bm{J}}^{(b)}$ are passed to the low-dimensional projection block for the calculation of $\widehat{\bm{\gamma}}_{r}^{(b)}(\lambda_{b})$. The resulting projection $\widehat{\bm{z}}_{r}^{(b)}$ and $\widehat{\tau}_{r}^{(b)}$ from the low-dimensional projection block together with the lasso estimator $\bm{\widehat{\beta}}^{(b)}$ will be used to compute the debiased lasso estimator $\widehat{\beta}^{(b)}_{\text{on},r}$ and its estimated standard error $\widehat{\sigma}_{r}^{(b)}$.

In this section, we conduct simulation studies to examine the finite-sample performance of the proposed ODL procedure in high-dimensional GLMs. We randomly generate a total of $N_{b}$ samples arriving in a sequence of $b$ data batches, denoted by $\{\mathcal{D}_{1},\dots,\mathcal{D}_{b}\}$, from the logistics regression model. Specifically,

where $\bm{x}_{i}^{(j)}\sim\mathcal{N}(\bm{0},\bm{\Sigma})$ and $\bm{\beta}^{0}\in\mathbb{R}^{p}$ is a $p$-dimensional sparse parameter vector. Recall that $s_{0}$ is the number of non-zero components of $\bm{\beta}^{0}$. We set half of the nonzero coefficients to be $1$ (relatively strong signals), and another half to be $0.01$ (weak signals). We consider the following settings: (i) $N_{b}=120$, $b=12$, $n_{j}=10$ for $j=1,\ldots,12$, $p=100$ and $s_{0}=6$; (ii) $N_{b}=624$, $b=12$, $n_{j}=52$ for $j=1,\ldots,12$, $p=600$ and $s_{0}=10$. Under each setting, two types of $\bm{\Sigma}$ are considered: (a) $\bm{\Sigma}=\bm{I}_{p}$; (b) $\bm{\Sigma}=\{0.5^{|i-j|}\}_{i,j=1,\dots,p}$. We set the step size in gradient descent $\eta=0.005$.

The objective is to conduct both estimation and inference along the arrival of a sequence of data batches. The evaluation criteria include: averaged absolute bias in estimating $\bm{\beta}^{0}$ (A.bias); averaged estimated standard error (ASE); empirical standard error (ESE); coverage probability (CP) of the 95% confidence intervals; averaged length of the 95% confidence interval (ACL). These metrics will be evaluated separately for three groups: (i) $\beta^{0}_{r}=0$, (ii) $\beta^{0}_{r}=0.01$ and (iii) $\beta^{0}_{r}=1$. Comparison is made among (i) the maximum likelihood estimator obtained by fitting the conventional GLM at the terminal point $b$ where $N_{b}>p$, (ii) the offline debiased $\ell_{1}$-penalized estimator at the terminal point $b$ which is also the benchmark method, and (iii) our proposed ODL estimator at several intermediate points from $j=1,\dots,b$. We include two offline methods in comparison using R packages hdi (Dezeure et al. ,, 2015) and glm, respectively. The results are reported in Tables 1-4.

It can be seen from Tables 1-4 that the estimation bias of the ODL estimator decreases rapidly as the number of data batches $b$ increasing from $2$ to $12$. Both the estimated standard errors and averaged length of 95% confidence intervals exhibit similar decreasing trend over time, and almost coincide with those by the offline benchmark method at the terminal points. Moreover, ODL enjoys great computation efficiency and it is almost 16 times faster than its offline counterpart.

It is worth noting that even though at the terminal point where the cumulative sample size $N_{b}$ is slightly larger than $p$, we can still fit the conventional GLM to obtain the MLE. It fails to provide reliable coverage probabilities due to severely large biases and estimated standard errors. In particular, the estimation bias of MLE becomes 50 times that of offline or online debiased lasso when $p=100$ as shown in Tables 1 and 2, and it further increases to 80 times the debiased estimators when $p=600$. Furthermore, as clearly indicated by the large empirical standard errors, MLE under this setting suffers from severe instability. Such an invalid estimation and inference results by MLE further demonstrates the advantage of our proposed online debiased method under the high-dimensional sparse logistic regression setting with streaming datasets.