# Online Debiased Lasso for Streaming Data

Ruijian Han 1, Lan Luo∗2, Yuanyuan Lin1 and Jian Huang2
1Department of Statistics, The Chinese University of Hong Kong, Hong Kong SAR,
China
2Department of Statistics and Actuarial Science, University of Iowa, Iowa
City, Iowa, USA  Ruijian Han and Lan Luo contributed equally to this work.

Abstract

We propose an online debiased lasso (ODL) method for statistical inference in
high-dimensional linear models with streaming data. The proposed ODL consists
of an efficient computational algorithm for streaming data and approximately
normal estimators for the regression coefficients. Its implementation only
requires the availability of the current data batch in the data stream and
sufficient statistics of the historical data at each stage of the analysis. A
dynamic procedure is developed to select and update the tuning parameters upon
the arrival of each new data batch so that we can adjust the amount of
regularization adaptively along the data stream. The asymptotic normality of
the ODL estimator is established under the conditions similar to those in an
offline setting and mild conditions on the size of data batches in the stream,
which provides theoretical justification for the proposed online statistical
inference procedure. We conduct extensive numerical experiments to evaluate
the performance of ODL. These experiments demonstrate the effectiveness of our
algorithm and support the theoretical results. An air quality dataset and an
index fund dataset from Hong Kong Stock Exchange are analyzed to illustrate
the application of the proposed method.

Key words: Adaptive tuning; Confidence interval; Debiased lasso; Gradient
descent; Online algorithm.

##  1 Introduction

The advent of distributed online learning systems such as Apache Flink
(Carbone et al. ,, 2015) has motivated new developments in data analytics for
streaming processing. Such systems enable efficient analyses of massive
streaming data assembled through, for example, mobile or web applications
(Jiang et al. ,, 2018), e-commerce purchases (Akter & Wamba,, 2016),
infectious disease surveillance programs (Choi et al. ,, 2016; Samaras et al.
,, 2020), mobile health consortia (Shameer et al. ,, 2017; Kraft et al. ,,
2020), and financial trading floors (Das et al. ,, 2018). Streaming data
refers to a data collection scheme where observations arrive sequentially and
perpetually over time, making it challenging to fit into computer memory for
static analyses. Researchers would query such continuous and unbounded data
streams in real-time to answer questions of interest including assessing
disease progression, monitoring product safety, and validating drug efficacy
and side effects. In these scenarios, it is essential for practitioners to
process data streams sequentially and incrementally as part of online
monitoring and decision-making procedures. Additionally, data streams from
various fields such as bioinformatics, medical imaging, and computer vision
are usually high-dimensional in nature.

In this paper, we consider the problem of online statistical inference in
high-dimensional linear regression with streaming data and propose an online
debiased lasso (ODL) estimator. While substantial advancements have been made
in online learning and associated optimization problems, the existing works
focus on online computational algorithms for point estimation and their
numerical convergence properties (Langford et al. ,, 2009; Duchi et al. ,,
2011; Tarrès & Yao,, 2014; Sun et al. ,, 2020). However, these works did not
consider the statistical distribution properties of the online point
estimators, which are needed for making statistical inference. Online
statistical inference methods have been mostly developed under low-dimensional
settings where $p\ll n$ (Schifano et al. ,, 2016; Luo & Song,, 2020). The goal
of this paper is to develop an online algorithm and statistical inference
procedure for analyzing high-dimensional streaming data.

###  1.1 Related work

The last decade has witnessed enormous progress on statistical inference in
high-dimensional models, see, for example, Zhang & Zhang, (2014); van de Geer
et al. , (2014); Javanmard & Montanari, (2014), as well as the review paper
Dezeure et al. , (2015) and the references therein. Most of the advancements
such as the novel debiased lasso have been developed for an offline setting. A
major difficulty in an online setting with streaming data is that one does not
have full access to the entire dataset as new data arrives on a continual
basis. To tackle the computational and inference problems due to the evolving
nature of the high dimensional stream, it is desirable to develop an algorithm
and statistical inference procedure in an online mode by updating the
regression parameters sequentially with newly arrived data batch and summary
statistics of historical raw data.

In recent years, there has been an ever-increasing interest in developing
online variable selection methods for high-dimensional streaming data. Most of
the work is along the line of lasso (Tibshirani,, 1996). For example, Langford
et al. , (2009) proposed an online $\ell_{1}$-regularized method via a variant
of the truncated SGD. Fan et al. , (2018) adopted the diffusion approximation
techniques to characterize the dynamics of the sparse online regression
process. Comprehensive development of online counterparts of popular offline
variable selection algorithms such as lasso, Elastic Net (Zou & Hastie,,
2005), Minimax Convex Penalty (MCP) (Zhang,, 2010), and Feature Selection with
Annealing (FSA) (Duchi et al. ,, 2011) has been studied by Sun et al. ,
(2020). Nevertheless, it is known that variable selection methods focus on
point estimation, but do not provide any uncertainty assessment. There is no
systematic study on statistical inference, including interval estimation and
hypothesis testing, with high-dimensional streaming data. Another complication
in dealing with high-dimensional streaming data is that, the regularization
parameter $\lambda$ that controls the sparsity level can no longer be
determined by the traditional cross-validation. Instead, small coefficients
are rounded to zero with a certain threshold or a pre-specified sparsity level
(Sun et al. ,, 2020). Recently, Deshpande et al. , (2019) considered a class
of online estimators in a high-dimensional auto-regressive model and studied
the asymptotic properties via martingale theories. Shi et al. , (2020)
proposed an inference procedure for high-dimensional linear models via
recursive online-score estimation. In both works, it is assumed that the
entire dataset is available at the initial stage for computing an initial
estimator (e.g. the lasso estimator) and the information in the streaming data
is used to reduce the bias of the initial estimator. However, the assumption
that the full dataset is available at the initial stage is not realistic in an
online learning setting.

###  1.2 Our contributions

The goal of this work is to develop an online debiased lasso estimator for
statistical inference with high-dimensional streaming data. Our proposed ODL
differs from the aforementioned works on online inference in two crucial
aspects. First, we do not assume the availability of the full dataset at the
initial stage. Second, at each stage of the analysis, we only require the
availability of the current data batch and sufficient statistics of historical
data. Therefore, ODL achieves statistical efficiency without accessing the
entire dataset. Furthermore, we propose a new approach for tuning parameter
selection that is naturally suited to the streaming data structure. In
addition, we provide a detailed theoretical analysis of the proposed ODL
estimator. The main contributions of the paper are as follows.

  * •

We introduce a new approach for online statistical inference in high-
dimensional linear models. Our proposed ODL consists of two main ingredients:
online lasso estimation and online debiasing lasso. Instead of re-accessing
the entire dataset, we utilize only sufficient statistics and the current data
batch.

  * •

We propose a new adaptive procedure to determine and update the tuning
parameter $\lambda$ dynamically upon the arrival of a new data batch, which
enables us to adjust the amount of regularization adaptively along with the
data accumulation. Our proposed online tuning parameter selector aligns with
the online estimation and debiasing procedures and involves summary statistics
only.

  * •

We establish the asymptotic normality of the proposed ODL estimator under the
conditions similar to those in an offline setting and mild conditions on the
batch sizes. We show that the asymptotic normality result holds as the
cumulative sample size goes to infinity, regardless of finite data batch
sizes. These results provide a theoretical basis for constructing confidence
intervals and conducting hypothesis tests with approximately correct
confidence levels and test sizes, respectively.

  * •

Extensive numerical experiments with simulated data demonstrate that ODL
algorithm is computationally efficient and strongly support the theoretical
properties of the ODL estimator. An air pollution dataset is also used to
illustrate the application of ODL.

The rest of the paper is organized as follows. Section 2 presents the model
formulation and our proposed ODL procedure. Section 3 includes the theoretical
properties of the proposed ODL estimator. Simulation experiments are given in
Section 4 to evaluate the performance of our proposed ODL in comparison to the
offline ordinary least square estimator. In Section 5.1 we demonstrate the
application of the proposed method on an air pollution dataset. Concluding
remarks are given in Section 6. Detailed proof of the theoretical properties
is included in the appendix.

##  2 Online debiased lasso

Consider a time point $b\geq 2$ with a total of $N_{b}$ samples arriving in a
sequence of $b$ data batches, denoted by
$\\{\mathcal{D}_{1},\dots,\mathcal{D}_{b}\\}$. Samples in each data batch
$\mathcal{D}_{j}=\\{\bm{y}^{(j)},\bm{X}^{(j)}\\}$ satisfy:

$\bm{y}^{(j)}=\bm{X}^{(j)}\bm{\beta}_{0}+\bm{\epsilon}^{(j)},\ j=1,\dots,b,$
(1)

where $\bm{y}^{(j)}=(y^{(j)}_{1},\dots,y^{(j)}_{n_{j}})^{\top}$ is the
response vector and
$\bm{X}^{(j)}=(\bm{x}^{(j)}_{1},\dots,\bm{x}^{(j)}_{n_{j}})^{\top}$ is an
$n_{j}\times p$ design matrix with $n_{j}$ being the data batch size. Here the
regression coefficient
$\bm{\beta}_{0}=(\beta_{0,1},\ldots,\beta_{0,p})^{\top}\in\mathbb{R}^{p}$ is
an unknown but sparse vector, and the error terms ${\epsilon}^{(j)}_{i}$,
$i=1,\dots,n_{j}$, are independent and identically distributed (i.i.d) with
mean zero and finite but unknown variance $\sigma^{2}_{\epsilon}$. Throughout
the paper, we consider a high-dimensional linear model, in particular, we
allow $p\geq N_{b}\equiv\sum_{j=1}^{b}n_{j}$.

In a streaming data setting where data volume accumulates fast over time,
individual-level raw data may not be stored in memory for a long time, making
it impossible to implement the offline debiased algorithms (Zhang & Zhang,,
2014; van de Geer et al. ,, 2014; Javanmard & Montanari,, 2014) that require
access to the entire dataset. To address this issue, we develop an online
debiasing procedure for each component of $\bm{\beta}_{0}$ in model (1).
Without loss of generality, our discussion in the following focuses on the
estimation and inference of $\beta_{0,r}$, the $r$-th component of
$\bm{\beta}_{0}$, $r=1,\ldots,p$.

When the first data batch $\mathcal{D}_{1}=\\{\bm{y}^{(1)},\bm{X}^{(1)}\\}$
arrives, we start off by applying the offline debiased lasso to obtain the
initial estimator. Specifically, let $\bm{x}_{r}^{(1)}$ be the $r$-th column
of $\bm{X}^{(1)}$ and $\bm{X}_{-r}^{(1)}$ be the sub-matrix of $\bm{X}^{(1)}$
excluding the $r$-th column. An initial lasso estimator is given by

$\bm{\widehat{\beta}}^{(1)}:=\underset{\bm{\beta}\in\mathbb{R}^{p}}{\arg\min}\
\left\\{\frac{1}{2n_{1}}\|\bm{y}^{(1)}-\bm{X}^{(1)}\bm{\beta}\|_{2}^{2}+\lambda_{1}\|\bm{\beta}\|_{1}\right\\},$
(2)

where $\lambda_{1}\geq 0$ is a regularization parameter. Following Zhang &
Zhang, (2014), to construct a confidence interval for $\beta_{0,r}$, a low-
dimensional projection $\bm{\widehat{z}}^{(1)}_{r}$ that acts as the
projection of $\bm{x}_{r}^{(1)}$ to the orthogonal complement of the column
space of $\bm{X}^{(1)}_{-r}$, is defined as
$\bm{\widehat{z}}^{(1)}_{r}:=\bm{x}_{r}^{(1)}-\bm{X}^{(1)}_{-r}\bm{\widehat{\gamma}}_{r}^{(1)}$,
where

$\bm{\widehat{\gamma}}^{(1)}_{r}:=\underset{\bm{\gamma}\in\mathbb{R}^{(p-1)}}{\arg\min}\left\\{\frac{1}{2n_{1}}\|\bm{x}^{(1)}_{r}-\bm{X}^{(1)}_{-r}\bm{\gamma}\|_{2}^{2}+\lambda_{1}\|\bm{\gamma}\|_{1}\right\\},$
(3)

and $\lambda_{1}$ is taken to be the same as in (2) for simplicity. Then, the
offline debiased lasso estimator of $\bm{\beta}_{0}$, $r=1,\ldots,p$, is
defined as

${\widehat{\beta}}^{(1)}_{\text{off},r}:={\widehat{\beta}}^{(1)}_{r}-\frac{(\bm{\widehat{z}}^{(1)}_{r})^{\top}(\bm{y}^{(1)}-\bm{X}^{(1)}\bm{\widehat{\beta}}^{(1)})}{(\bm{\widehat{z}}^{(1)}_{r})^{\top}\bm{x}_{r}^{(1)}}.$
(4)

Later on, when the second batch
$\mathcal{D}_{2}=\\{\bm{y}^{(2)},\bm{X}^{(2)}\\}$ arrives, the offline
debiased lasso algorithm would replace $\bm{y}^{(1)}$ and $\bm{X}^{(1)}$ by
the augmented full dataset $\left\\{\bm{y}^{(1)},\bm{y}^{(2)}\right\\}$ and
$\left\\{\bm{X}^{(1)},\bm{X}^{(2)}\right\\}$ respectively in (2)-(4). However,
$\\{\bm{y}^{(1)},\bm{X}^{(1)}\\}$ may no longer be available in an online
setting. To address this issue, we propose an online estimation and inference
procedure that utilize the information in historical raw data via summary
statistics. We present the three main steps of ODL in Subsections 2.1-2.3
below.

###  2.1 Online lasso

Upon the arrival of data batch
$\mathcal{D}_{b}=\\{\bm{y}^{(b)},\bm{X}^{(b)}\\}$ with $b\geq 2$, if all
previous data batches $\\{\mathcal{D}_{1},\ldots,\mathcal{D}_{b-1}\\}$ are
available, we can use the offline lasso method that solves the following
optimization problem:

$\bm{\widehat{\beta}}^{(b)}(\lambda_{b}):=\underset{\bm{\beta}\in\mathbb{R}^{p}}{\arg\min}\
\left\\{\frac{1}{2N_{b}}\sum_{j=1}^{b}\|\bm{y}^{(j)}-\bm{X}^{(j)}\bm{\beta}\|_{2}^{2}+\lambda_{b}\|\bm{\beta}\|_{1}\right\\},$
(5)

where $N_{b}=\sum_{j=1}^{b}n_{j}$ is the cumulative sample size and
$\lambda_{b}$ is the regularization parameter adaptively chosen for step $b$.
For simplicity, we use $\bm{\widehat{\beta}}^{(b)}$ to denote
$\bm{\widehat{\beta}}^{(b)}(\lambda_{b})$ except in the discussion on the
choice of $\lambda_{b}$ in Section 2.4. However, since we only assume the
availability of summary statistics of historical data, we cannot use the
algorithms such as coordinate descent (Friedman et al. ,, 2007) that requires
the availability of the whole dataset.

Note that the objective function in (5) depends on the data only through the
following summary statistics:

$\bm{{S}}^{(b)}\equiv\sum_{j=1}^{b}(\bm{X}^{(j)})^{\top}\bm{X}^{(j)},\
\bm{U}^{(b)}\equiv\sum_{j=1}^{b}(\bm{X}^{(j)})^{\top}\bm{y}^{(j)}.$ (6)

Therefore, it is reasonable to refer to these summary statistic as sufficient
statistics relative to this objective function, or simply sufficient
statistics without confusing with the concept of sufficient statistic in a
parametric model. Based on these two summary statistics, one can obtain the
solution to (5) by the gradient descent algorithm.

Let
$\mathcal{L}(\bm{\beta})=\sum_{j=1}^{b}\|\bm{y}^{(j)}-\bm{X}^{(j)}\bm{\beta}\|_{2}^{2}/(2N_{b})$,
and the gradient of $\mathcal{L}(\bm{\beta})$ is given by

$\frac{\partial\mathcal{L}(\bm{\beta})}{\partial\bm{\beta}}=\frac{1}{N_{b}}(\bm{{S}}^{(b)}\bm{\beta}-\bm{U}^{(b)}).$
(7)

Notably, the gradient depends on the historical raw data only through the
summary statistics $\bm{{S}}^{(b)}$ and $\bm{U}^{(b)}$. We update the solution
iteratively by combining a gradient descent step and soft thresholding
(Daubechies et al. ,, 2004; Donoho & Johnstone,, 1994). Specifically, each
iteration consists of the following two steps:

  * •

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

$\widehat{\bm{\beta}}^{(b)}\leftarrow\widehat{\bm{\beta}}^{(b)}-\eta\frac{\partial\mathcal{L}(\bm{\beta})}{\partial\bm{\beta}}=\widehat{\bm{\beta}}^{(b)}-\frac{\eta}{N_{b}}(\bm{S}^{(b)}\widehat{\bm{\beta}}^{(b)}-\bm{U}^{(b)}),$
(8)

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

  * •

Step 2 (Soft thresholding): apply the soft-thresholding operator
${\mathcal{S}}(\widehat{\beta}^{(b)}_{r};\eta\lambda_{b})$ to the $r$-th
component in $\widehat{\bm{\beta}}^{(b)}$ in step 1, for $r=1,\dots,p$, where
$\mathcal{S}(x,\lambda)=\text{sgn}(x)(|x|-\lambda)_{+}$.

These two steps are carried out iteratively till convergence. In the
implementation, the stopping criterion is set as
$\|{\partial\mathcal{L}(\bm{\beta})}/{\partial\bm{\beta}}\|_{2}\leq 10^{-6}$.

The size of the sufficient statistics will not increase as more data batches
arrive. For example, when a new batch $\mathcal{D}_{b}$ arrives, we update the
summary statistics by

$\bm{S}^{(b)}=\bm{S}^{(b-1)}+(\bm{X}^{(b)})^{\top}\bm{X}^{(b)},\
\bm{U}^{(b)}=\bm{U}^{(b-1)}+(\bm{X}^{(b)})^{\top}\bm{y}^{(b)}.$

incrementally, which are matrices of fixed dimensions even if $b\to\infty$. In
addition, a consistent estimator of $\sigma_{\epsilon}^{2}$ by the method of
moments is given by

$(\widehat{\sigma}^{2}_{\epsilon})^{(b)}\equiv\frac{N_{b-1}}{N_{b}}(\widehat{\sigma}_{\epsilon}^{2})^{(b-1)}+\frac{n_{b}}{N_{b}}(\bm{y}^{(b)}-\bm{X}^{(b)}\widehat{\bm{\beta}}^{(b)})^{\top}(\bm{y}^{(b)}-\bm{X}^{(b)}\widehat{\bm{\beta}}^{(b)}),$
(9)

which will be used in constructing the confidence intervals.

###  2.2 Online low-dimensional projection

Next, we obtain an online estimator for the low-dimensional projection. Let

$\bm{\widehat{\gamma}}_{r}^{(b)}:=\underset{\bm{\gamma}\in\mathbb{R}^{(p-1)}}{\arg\min}\left\\{\frac{1}{2N_{b}}\sum_{j=1}^{b}\|\bm{x}^{(j)}_{r}-\bm{X}^{(j)}_{-r}\bm{\gamma}\|_{2}^{2}+\lambda_{b}\|\bm{\gamma}\|_{1}\right\\},$
(10)

where $N_{b}$ and $\lambda_{b}$ are the same as in (5). We can summarize the
data information in the following two statistics:
$\bm{R}^{(b)}=\sum_{j=1}^{b}(\bm{X}^{(j)}_{-r})^{\top}\bm{X}^{(j)}_{-r},\
\bm{T}^{(b)}=\sum_{j=1}^{b}(\bm{X}^{(j)}_{-r})^{\top}\bm{x}^{(j)}_{r}.$
Repeating similar procedure in the online lasso, we can obtain
$\bm{\widehat{\gamma}}_{r}^{(b)}$ and further define a low-dimensional
projection
$\bm{\widehat{z}}^{(b)}_{r}:=\bm{x}^{(b)}_{r}-\bm{X}^{(b)}_{-r}\bm{\widehat{\gamma}}^{(b)}_{r}.$
It is worth mentioning that $\bm{R}^{(b)}$ and $\bm{T}^{(b)}$ are obtained
from $\bm{S}^{(b)}$ directly with $\bm{R}^{(b)}=\bm{S}^{(b)}_{-r,-r}$ and
$\bm{T}^{(b)}=\bm{S}^{(b)}_{-r,r}$, where $\bm{S}^{(b)}_{-r,-r}$ is a sub-
matrix of $\bm{S}^{(b)}$ excluding the $r$-th row and the $r$-th column, and
$\bm{S}^{(b)}_{-r,r}$ is a sub-matrix of $\bm{S}^{(b)}$ with the $r$-th row
being deleted but the $r$-th column being kept. The low-dimensional projection
$\bm{\widehat{z}}^{(b)}_{r}$ will be used in constructing the debiased
estimator in Subsection 2.3.

###  2.3 Online debiased lasso estimator

When data batch $\mathcal{D}_{b}$ arrives, the ODL estimator for
$\beta_{0,r}$, $r=1,\ldots,p$, is defined as

$\begin{split}{\widehat{\beta}}^{(b)}_{\text{on},r}:=\widehat{{\beta}}^{(b)}_{r}+\left\\{\sum_{j=1}^{b}(\bm{\widehat{z}}^{(j)}_{r})^{\top}\bm{x}^{(j)}_{r}\right\\}^{-1}\left\\{\sum_{j=1}^{b}(\bm{\widehat{z}}^{(j)}_{r})^{\top}\bm{y}^{(j)}-\sum_{j=1}^{b}(\bm{\widehat{z}}^{(j)}_{r})^{\top}\bm{X}^{(j)}\bm{\widehat{\beta}}^{(b)}\right\\}.\end{split}$
(11)

Although the historical data are still involved in (11), we only need to store
the following statistics rather than the entire dataset to compute
${\widehat{\beta}}^{(b)}_{\text{on},r}$. Specifically, let
$a_{1}^{(b)}:=\sum_{j=1}^{b}(\bm{\widehat{z}}^{(j)}_{r})^{\top}\bm{x}^{(j)}_{r},\
a_{2}^{(b)}:=\sum_{j=1}^{b}(\bm{\widehat{z}}^{(j)}_{r})^{\top}\bm{y}^{(j)},\
\bm{A}_{1}^{(b)}:=\sum_{j=1}^{b}(\bm{\widehat{z}}^{(j)}_{r})^{\top}\bm{X}^{(j)},$
which have the same dimensions when the new data arrives and can be easily
updated. For example, we update $a_{1}^{(b)}$ by
$a_{1}^{(b)}=a_{1}^{(b-1)}+(\bm{\widehat{z}}^{(b)}_{r})^{\top}\bm{x}^{(b)}_{r}.$
Consequently, by substituting the online lasso estimator and low-dimensional
projection into (11), we obtain the ODL estimator
${\widehat{\beta}}^{(b)}_{\text{on},r}.$ Meanwhile, as discussed in Section 3,
the estimated standard error is given by
$\widehat{\sigma}^{(b)}_{\epsilon}\widehat{\tau}^{(b)}_{r}$, where

$\widehat{\tau}^{(b)}_{r}=\frac{\sqrt{\sum_{j=1}^{b}(\bm{\widehat{z}}^{(j)}_{r})^{\top}\bm{\widehat{z}}^{(j)}_{r}}}{\sum_{j=1}^{b}(\bm{\widehat{z}}^{(j)}_{r})^{\top}\bm{x}^{(j)}_{r}},$
(12)

and $(\widehat{\sigma}^{2}_{\epsilon})^{(b)}$ is given in (9) in Subsection
2.1.

###  2.4 Tuning parameter selection

In an offline setting, the tuning parameter $\lambda$ can be chosen from a
candidate set via cross-validation where the entire dataset is split into
training and testing sets multiple times. However, such a sample-splitting
scheme is not applicable in an online setting since we do not have the full
dataset at hand. A natural sample-splitting idea that aligns with the
streaming data structure originates from the forecasting accuracy evaluation
in time series; see Figure 1. At time point $b$, those sequentially arrived
data batches up to time point $b-1$, denoted by
$\\{\mathcal{D}_{1},\dots,\mathcal{D}_{b-1}\\}$, serve as the training set,
and the current data batch $\mathcal{D}_{b}$ is the testing set. This
procedure is also known as “rolling-original-recalibration” (Tashman,, 2000).

Specifically, with only the first data batch $\mathcal{D}_{1}$,
$\bm{\widehat{\beta}}^{(1)}$ is a standard lasso estimator with $\lambda_{1}$
selected by the classical offline cross-validation. When the $b$-th data batch
$\mathcal{D}_{b}$ arrives, we calculate

$PE_{b}(\lambda)=\frac{1}{n_{b}}\|(\bm{y}^{(b)}-\bm{X}^{(b)}\bm{\widehat{\beta}}^{(b-1)}(\lambda)\|_{2}^{2},~{}\lambda\in
T_{\lambda},$

and define

$\lambda_{b}:=\arg\underset{\lambda\in T_{\lambda}}{\min}\ PE_{b}(\lambda).$

In such a way, we are able to determine $\lambda_{b}$ upon the arrival of a
new data batch $\mathcal{D}_{b}$ adaptively, and extract the corresponding
lasso estimator $\widehat{\bm{\beta}}^{(b)}(\lambda_{b})$ as the starting
point for ODL.

Figure 1: A diagram illustrates the series of training and test sets, where
the black dots form the training sets, and the gray dots form the test sets.
At a time point $b$, the ODL estimator $\widehat{\bm{\beta}}^{(b-1)}(\lambda)$
is obtained based on the training set
$\\{\mathcal{D}_{1},\dots,\mathcal{D}_{b-1}\\}$ and the current data batch
$\mathcal{D}_{b}=\\{\bm{y}^{(b)},\bm{X}^{(b)}\\}$ is the testing set.

###  2.5 Summary

We now summarize our proposed ODL procedure for the statistical inference of
$\beta_{0,r},r=1,\ldots,p,$ using a flowchart in Figure 2. 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\\{\bm{S}^{(b-1)},\bm{U}^{(b-1)}\right\\}$ are
updated to $\left\\{\bm{S}^{(b)},\bm{U}^{(b)}\right\\}$. These summary
statistics facilitate the updating of the lasso estimator
$\widehat{\bm{\beta}}^{(b-1)}$ to $\widehat{\bm{\beta}}^{(b)}$ at some grid
values of the tuning parameters 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}$ is passed to the low-dimensional projection block for the
calculation of $\widehat{\bm{\gamma}}_{r}^{(b)}(\lambda_{b})$. The idea of
online updating is the same as in lasso estimation, except the relevant
summary statistics are the sub-matrices of $\bm{S}^{(b)}$. The resulting
projection $\widehat{\bm{z}}_{r}^{(b)}$ output from the low-dimensional
projection block together with the lasso estimator
$\bm{\widehat{\beta}}^{(b)}(\lambda_{b})$ will be used to compute the debiased
lasso estimator $\widehat{\beta}^{(b)}_{\text{on},r}$ and its estimated
standard error.

Figure 2: Flowchart of the online debiasing algorithm. When a new data batch
$\mathcal{D}_{b}$ arrives, it is sent to the lasso estimation block for
updating $\widehat{\bm{\beta}}^{(b-1)}$ to $\widehat{\bm{\beta}}^{(b)}$. At
the same time, it is also viewed as test set for adaptively choosing tuning
parameter $\lambda_{b}$. In the low-dim projection block, we extract sub-
matrices from the updated summary statistic $\bm{S}^{(b)}$ to compute
$\widehat{\bm{\gamma}}_{r}^{(b)}(\lambda_{b})$ and the corresponding low-
dimensional projection $\widehat{\bm{z}}_{r}^{(b)}$. Outputs
$\widehat{\beta}_{r}^{(b)}(\lambda_{b})$ and $\widehat{\bm{z}}_{r}^{(b)}$ from
the two blocks are further used to compute the debiased lasso estimator
$\widehat{\beta}_{\text{on},r}^{(b)}$ and its estimated standard error
$\widehat{\text{SE}}(\widehat{\beta}_{\text{on},r}^{(b)})$.

######  Remark 1.

When $p$ is large, the online algorithm presented in Figure 2 requires a
sufficient large storage capacity, since sample covariance matrix
$\bm{S}^{(b)}$ requires $\mathcal{O}(p^{2})$ space complexity. To reduce
memory usage, we can apply the eigenvalue decomposition (EVD) of
$\bm{S}^{(b)}=Q_{b}\Lambda_{b}Q_{b}^{\top}$, where $Q_{b}$ is the $p\times
N_{b}$ orthogonal matrix combined by the eigenvectors, $\Lambda_{b}$ is the
$N_{b}\times N_{b}$ diagonal matrix whose diagonal elements are the
eigenvalues of $\bm{S}^{(b)}$. We only need to store $Q_{b}$ and
$\Lambda_{b}.$ Since $r_{b}=$ rank$(\Lambda_{b})\leq\min\\{N_{b},p\\},$ we can
use an incremental EVD approach (Cardot & Degras,, 2018) to update $Q_{b}$ and
$\Lambda_{b}$. Then the space complexity reduces to $\mathcal{O}(r_{b}p)$. The
space complexity can be further reduced by setting a threshold. For example,
select the principal components which explain most of the variations in the
predictors. However, incremental EVD could increase the computational cost
since it requires additional $\mathcal{O}(r_{b}^{2}p)$ computational
complexity. Indeed, there is a trade-off between the space complexity and
computational complexity. How to balance this trade-off is an important
computational issue and deserves careful analysis, but is beyond the scope of
this study.

##  3 Theoretical properties

To establish the asymptotic properties of the ODL estimator proposed in
Section 2, we first introduce some notation. Consider a random design matrix
$\bm{X}$ with i.i.d rows. Let $\bm{\Sigma}$ be the covariance matrix of each
row of $\bm{X}$. Denote the inverse of $\bm{\Sigma}$ by
$\bm{\Theta}=\bm{\Sigma}^{-1}.$ For $r=1,\ldots,p$, define

$\bm{\gamma}_{r}:=\arg\min_{\bm{\gamma}\in\mathbb{R}^{p-1}}\mathbb{E}[\lVert\bm{x}_{r}-\bm{X}_{-r}\bm{\gamma}\lVert_{2}^{2}],$

and the corresponding residual vector is
$\bm{z}_{r}:=\bm{x}_{r}-\bm{X}_{-r}\bm{\gamma}_{r}.$ Let
$s_{0}=|\\{j:{\bm{\beta}}_{j}\neq 0\\}|$ and $s_{r}=|\\{k\neq
r:\bm{\Theta}_{k,r}\neq 0\\}|$ be two sparsity levels.

The following regularity conditions on the design matrix $\bm{X}$ and the
error terms are imposed to establish the asymptotic results. Specifically, we
assume that $\bm{X}$ has either i.i.d sub-Gaussian or bounded rows. We first
consider the sub-Gaussian case.

######  Assumption 1.

Suppose that

(A1) The design matrix $\bm{X}$ has i.i.d sub-Gaussian rows.

(A2) The smallest eigenvalue $\Lambda^{2}_{\min}$ of $\bm{\Sigma}$ is strictly
positive and $1/\Lambda^{2}_{\min}=O(1).$ In addition, the largest diagonal
element of $\bm{\Sigma}$, $\max_{j}\Sigma_{j,j}={O}(1).$

(A3) The error terms ${\epsilon}^{(j)}_{i}$, $i\in\mathcal{D}_{j}$,
$j=1,\dots,b$ are sub-exponential.

######  Theorem 1.

Assume Assumption 1 holds. For the $j$-th data batch, suppose that the tuning
parameter $\lambda_{j}$ satisfies $\lambda_{j}=C\sqrt{{\log p}/{N_{j}}}$,
$j=1,\ldots,b$. If the first batch size ${n_{1}}\geq cs_{r}\log p$, the
subsequent batch size $n_{j}\geq c\log p,j=2,\ldots,b$, for some constant $c$,
and

$\displaystyle s_{0}s_{r}\frac{\log(p)}{\sqrt{N_{b}}}=o(1),\
s_{r}^{2}{\frac{\log(p)}{N_{b}}\log{\frac{N_{b}}{n_{1}}}}=o(1),$ (13)

then, for any $r=1,\ldots,p$ and sufficiently large $N_{b}$,

$\displaystyle({\widehat{\beta}}^{(b)}_{\text{on},r}-{\beta}_{0,r})/\widehat{\tau}^{(b)}_{r}=W_{r}+\Delta_{r},$
$\displaystyle
W_{r}=\frac{\bm{\widehat{z}}_{r}^{\top}\bm{\epsilon}}{\lVert\bm{\widehat{z}}_{r}\lVert_{2}},\
|\Delta_{r}|=o_{\mathbb{P}}(1),$

where $\widehat{\tau}^{(b)}_{r}$ is defined in (12).

######  Remark 2.

Similar to the offline debiased lasso estimator (Zhang & Zhang,, 2014; van de
Geer et al. ,, 2014), Theorem 1 implies that the dimensionality $p$ could be
at the exponential rate of the data size. However, the problem here is more
difficult than that in the offline setting and the proofs for the properties
of the offline debiased estimator do no apply here. Specifically, let
$\bm{\widetilde{z}}_{r}=((\bm{\widetilde{z}}^{(1)}_{r})^{\top},\ldots,(\bm{\widetilde{z}}^{(b)}_{r})^{\top})^{\top}$
be the low-dimensional projection in the offline case, where
$\bm{\widetilde{z}}^{(j)}_{r}=\bm{x}_{r}^{(j)}-\bm{X}_{-r}^{(j)}\bm{\widehat{\gamma}}^{(b)}_{r}$
is computed based on the $j$-th batch data, $j=1,\ldots,b$. Here,
$\bm{\widehat{\gamma}}^{(b)}_{r}$ depends on the historical data
$\\{\mathcal{D}_{1},\dots,\mathcal{D}_{b}\\}$. In contrast, the proposed
online low-dimensional projection
$\bm{\widehat{z}}^{(j)}_{r}=\bm{x}^{(j)}_{r}-\bm{X}^{(j)}_{-r}\bm{\widehat{\gamma}}^{(j)}_{r}$,
where $\bm{\widehat{\gamma}}^{(j)}_{r}$ is obtained solely based on the
$\\{\mathcal{D}_{1},\dots,\mathcal{D}_{j}\\}$. Thus the KKT condition for the
lasso minimization problem does not hold for the online estimator. The
arguments for the asymptotic properties of the debiased estimator in Zhang &
Zhang, (2014) and van de Geer et al. , (2014) heavily use the KKT condition.
Therefore, different arguments are needed to establish Theorem 1.

######  Remark 3.

Theorem 1 is established for the proposed online debiased lasso estimators
based on the algorithm described in Section 2. Indeed, the proof of Theorem 1
uses the specific form of the algorithm. Therefore, this result does not apply
to other online estimators computed using a different algorithm. For example,
it is not clear whether the estimators based on the online algorithms in
Langford et al. , (2009) and Fan et al. , (2018) will have similar asymptotic
distributional properties.

######  Remark 4.

The error terms are assumed to have sub-exponential tails in (A3). For sub-
Gaussian design matrix $\bm{X}$ in (A1), the assumption (A3) is the same as
that in the offline setting for the asymptotic properties of the debiased
lasso estimator in van de Geer et al. , (2014).

The requirement on the minimum batch size in Theorem 1 indicates that, one may
apply the online lasso algorithm once the sample size of the first data batch
reaches the order of $s_{r}\log p$. After that, we update the lasso estimators
when the size of the newly arrived batch is at the order of $\log p$. The next
theorem justifies that the order of the subsequent batch size $O(\log p)$
could be relaxed to $O(1)$, at the price of a relatively stronger condition on
$N_{b}$.

######  Theorem 2.

Assume Assumption 1 holds. When the $j$-th batch data arrives, suppose that
the tuning parameter $\lambda_{j}$ satisfies $\lambda_{j}=C\sqrt{{\log
p}/{N_{j}}}$. If the first batch size ${n_{1}}\geq cs_{r}\log p$ for some
constant $c$ and

$\displaystyle s_{0}s_{r}\sqrt{\frac{\log^{3}(p)}{N_{b}}}=o(1),\
s_{r}^{2}{\frac{\sqrt{\log^{3}(p)}}{N_{b}}\log{\frac{N_{b}}{n_{1}}}}=o(1),$

then, for any $r=1,\ldots,p$ and sufficiently large $N_{b}$,

$\displaystyle({\widehat{\beta}}^{(b)}_{\text{on},r}-{\beta}_{0,r})/\widehat{\tau}^{(b)}_{r}=W_{r}+\Delta_{r},$
$\displaystyle
W_{r}=\frac{\bm{\widehat{z}}_{r}^{\top}\bm{\epsilon}}{\lVert\bm{\widehat{z}}_{r}\lVert_{2}},\
|\Delta_{r}|=o_{\mathbb{P}}(1),$

where $\widehat{\tau}^{(b)}_{r}$ is defined in (12).

######  Remark 5.

The requirement of the first batch size $n_{1}\geq c_{1}s_{r}\log p$ in
Theorems 1 and 2 is needed to establish the consistency of the lasso-typed
estimator (Bühlmann & van de Geer,, 2011); otherwise, the error bound of
$\bm{\widehat{\gamma}}^{(1)}_{r}$ defined in (10) in the first step cannot be
controlled, resulting in large error (diverges as $N_{b}\to\infty$) in the
projection $\bm{\widehat{z}}^{(b)}_{r}$. When there is not enough data at the
initial stage, e.g., $n_{1}=\log p$, the error bound of
$\bm{\widehat{\gamma}}^{(1)}_{r}$ can also be controlled by considering some
bounded parameter space such as $\\{\bm{\gamma}:\|\bm{\gamma}\|_{1}\leq C\\}$
for some large constant $C$ rather than
$\\{\bm{\gamma}:\bm{\gamma}\in\mathbb{R}^{p-1}\\}$.

We now consider the case when the covariates are bounded. For a matrix
$\bm{A}=(a_{ij})$, let $\|\bm{A}\|_{\infty}$ be the largest absolute value of
its elements, that is, $\|\bm{A}\|_{\infty}=\max_{i,j}|a_{ij}|$.

######  Assumption 2.

Suppose that

(B1) The covariates are bounded by a finite constant $K>0$, namely,
$\lVert\bm{X}\lVert_{\infty}\leq K,$ where $\bm{X}$ is the design matrix.

(B2) The smallest eigenvalue $\Lambda^{2}_{\min}$ of $\bm{\Sigma}$ is strictly
positive and $1/\Lambda^{2}_{\min}=O(1).$ Moreover,
$\max_{j}\Sigma_{j,j}=O(1).$

(B3) $\lVert\bm{X}_{-r}\bm{\gamma}_{r}\lVert_{\infty}=O(K)$ and
$\max_{r}\mathbb{E}({z}^{4}_{r,1})=O(K^{4}),$ where ${z}_{r,1}$ is the first
element of $\bm{z}_{r}:=(\bm{x}_{r}-\bm{X}_{-r}\bm{\gamma}_{r})$.

######  Theorem 3.

Assume Assumption 2 holds. When the $j$-th batch data arrives, suppose that
the tuning parameter $\lambda_{j}$ satisfies $\lambda_{j}=C\sqrt{{\log
p}/{N_{j}}}$. If the first batch size ${n_{1}}\geq cs_{r}^{2}\log p$ for some
constant $c$ and

$\displaystyle s_{0}s_{r}\frac{\log(p)}{\sqrt{N_{b}}}=o(1),\
s_{r}^{2}{\frac{\log(p)}{N_{b}}\log{{N_{b}}}}=o(1),$

then, for any $r=1,\ldots,p$ and sufficiently large $N_{b}$,

$\displaystyle({\widehat{\beta}}^{(b)}_{\text{on},r}-{\beta}_{0,r})/\widehat{\tau}^{(b)}_{r}=W_{r}+\Delta_{r},$
$\displaystyle
W_{r}=\frac{\bm{\widehat{z}}_{r}^{\top}\bm{\epsilon}}{\lVert\bm{\widehat{z}}_{r}\lVert_{2}},\
|\Delta_{r}|=o_{\mathbb{P}}(1),$

where $\widehat{\tau}^{(b)}_{r}$ is defined in (12).

Theorems 2 and 3 are established without specific assumptions on data batch
sizes except for the first batch. Comparing to Theorem 2, Theorem 3 requires a
relatively stronger condition on $n_{1}$, but a more relaxed condition on the
cumulative sample size $N_{b}$. Furthermore, rewriting
$(\widehat{\sigma}^{2}_{\epsilon})^{(b)}$ as
$(\widehat{\sigma}^{2}_{\epsilon})^{(b)}=({1}/{N_{b}})\sum_{j=1}^{b}(\bm{y}^{(j)}-\bm{X}^{(j)}\widehat{\bm{\beta}}^{(j)})^{\top}(\bm{y}^{(j)}-\bm{X}^{(j)}\widehat{\bm{\beta}}^{(j)}),$
we can see that $\widehat{\sigma}^{(b)}_{\epsilon}$ is consistent for
$\sigma_{\epsilon}$ in view of the consistency of $\widehat{\bm{\beta}}^{(b)}$
in Lemma 2.

The proofs of Theorems 1–3 are included in the appendix. According to Theorems
1 \- 3, $W_{r}$ is asymptotically normal through verifying the conditions of
the Lindeberg central limit theorem. As a result, for any $0<\alpha<1$, a
$(1-\alpha)\%$ confidence interval for $\beta_{0,r}$ is

${\widehat{\beta}}^{(b)}_{\text{on},r}\pm\Phi^{-1}(1-\frac{\alpha}{2})(\widehat{\sigma}^{(b)}_{\epsilon}\widehat{\tau}^{(b)}_{r}),$

where $\widehat{\sigma}^{(b)}_{\epsilon}$ is defined in (9), $\Phi(\cdot)$ is
the cumulative distribution function of the standard normal distribution and
$\Phi^{-1}$ is its inverse function.

##  4 Simulation experiments

###  4.1 Setup

In this section, we conduct simulation studies to examine the finite-sample
performance of the proposed online debiasing procedure in high-dimensional
linear models. 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

$y^{(j)}_{i}=\bm{\beta}_{0}^{\top}\bm{x}_{i}^{(j)}+\epsilon_{i}^{(j)},\
i=1,\dots,n_{j};\ j=1,\dots,b,$

where $\epsilon_{i}\overset{iid}{\sim}\mathcal{N}({0},\sigma^{2}_{\epsilon})$,
$\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}=420$, $b=12$, $n_{j}=35$ for
$j=1,\ldots,12$, $p=400$ and $s_{0}=6$; (ii) $N_{b}=1200$, $b=12$, $n_{j}=100$
for $j=1,\ldots,12$, $p=1000$ and $s_{0}=20$. 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$ in case (i) and $\eta=0.05$ in case (ii).

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 quantities 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 between our proposed online debiased lasso at several intermediate
points from $j=1,\dots,b$ and the ordinary least squares (OLS) estimator at
the terminal point $b$ where $N_{b}>p$. We include the OLS method using R
package lm as a benchmark for comparison. The results are reported in Tables
1-4.

###  4.2 Bias and coverage probability

It can be seen from Tables 1-4 that the estimation bias of the online debiased
lasso 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. Even though
the coverage probabilities of the confidence intervals by the OLS at the end
point are around the nominal 95% level, both the estimation bias and standard
errors of OLS estimator are much larger than those of online debiased lasso.
In particular, the estimation bias of OLS could even be 10 times that of
online debiased lasso when $p=400$ as shown in Tables 1 and 2. Furthermore, it
is worth noting that even though the coverage probability of both estimators
reaches the nominal level at the terminal point, the ACL of OLS is about 2 to
4 times the one of online debiased lasso. Such a loss of statistical
efficiency by OLS further demonstrates the advantage of our proposed online
debiased method under the high-dimensional sparse regression setting with
streaming datasets.

To visualize the asymptotic normality of our proposed online debiasing
estimator, we plot the proposed online debiased lasso estimates at several
intermediate points $b=2,6,10$ against the theoretical quantiles of a standard
normal distribution in Figure 3. Similar plots for more settings with
$N_{b}=1200$ and $p=1000$ are provided in the Appendix. In these Q-Q plots,
the scattered points summarized from 200 replications stay closely along the
$45^{\circ}$ diagonal blue line, indicating the validity of asymptotic normal
distribution. Furthermore, such trend becomes clearer as $b$ increases.

By comparing across different signal groups, i.e. $\beta_{0,r}=0,0.01,1$, we
observe that both ASE and ACL are quite close to each other and even coincide
when $N_{b}=1200$, as shown in Tables 3-4. We believe this is reasonable, as
each column in the design matrix, denoted by
$\bm{x}_{r}\in\mathbb{R}^{N_{b}\times 1}$, is of the same marginal
distribution, and thus the estimated standard errors computed according to
(12) in the simulations are identical up to a certain decimal for every
component in $\bm{\beta}$, regardless of signal strength.

Table 1:  $N_{b}=420$, $b=12$, $p=400$, $s_{0}=6$, $\bm{\Sigma}=\bm{I}_{p}$. Performance on statistical inference. Tuning parameter $\lambda$ is chosen from $T_{\lambda}=\\{0.15,0.20,0.25,0.30\\}$ using the adaptive method in Section 2.4. Simulation results are summarized over 200 replications. In the table, we report the $\lambda$ selected with highest frequency among 200 replications. | $\beta_{0,r}$ | OLS | online debiased lasso
---|---|---|---
data batch index |  |  | 2 | 4 | 6 | 8 | 10 | 12
$\lambda$ |  |  | 0.30 | 0.30 | 0.20 | 0.15 | 0.15 | 0.15
A.bias | 0 | 0.013 | 0.008 | 0.005 | 0.004 | 0.004 | 0.003 | 0.003
0.01 | 0.026 | 0.007 | 0.009 | 0.008 | 0.003 | 0.002 | 0.002
1 | 0.016 | 0.152 | 0.022 | 0.009 | 0.005 | 0.002 | 0.001
ASE | 0 | 0.223 | 0.119 | 0.091 | 0.074 | 0.063 | 0.056 | 0.051
0.01 | 0.226 | 0.119 | 0.091 | 0.074 | 0.063 | 0.056 | 0.051
1 | 0.222 | 0.118 | 0.091 | 0.074 | 0.063 | 0.056 | 0.051
ESE | 0 | 0.229 | 0.140 | 0.096 | 0.073 | 0.062 | 0.055 | 0.051
0.01 | 0.235 | 0.014 | 0.094 | 0.070 | 0.060 | 0.057 | 0.051
1 | 0.226 | 0.171 | 0.102 | 0.079 | 0.067 | 0.057 | 0.053
CP | 0 | 0.934 | 0.901 | 0.936 | 0.953 | 0.953 | 0.952 | 0.951
0.01 | 0.947 | 0.903 | 0.943 | 0.952 | 0.955 | 0.943 | 0.943
1 | 0.940 | 0.683 | 0.902 | 0.933 | 0.935 | 0.948 | 0.948
ACL | 0 | 0.874 | 0.465 | 0.356 | 0.289 | 0.247 | 0.220 | 0.199
0.01 | 0.886 | 0.467 | 0.356 | 0.289 | 0.248 | 0.221 | 0.200
1 | 0.871 | 0.464 | 0.356 | 0.289 | 0.247 | 0.220 | 0.199
Table 2:  $N_{b}=420$, $b=12$, $p=400$, $s_{0}=6$, $\bm{\Sigma}=\\{0.5^{|i-j|}\\}_{i,j=1,\dots,p}$. Performance on statistical inference. Tuning parameter $\lambda$ is chosen from $T_{\lambda}=\\{0.15,0.20,0.25,0.30\\}$ using the adaptive method in Section 2.4. Simulation results are summarized over 200 replications. In the table, we report the $\lambda$ selected with highest frequency among 200 replications. | $\beta_{0,r}$ | OLS | online debiased lasso
---|---|---|---
data batch index |  |  | 2 | 4 | 6 | 8 | 10 | 12
$\lambda$ |  |  | 0.30 | 0.30 | 0.20 | 0.15 | 0.15 | 0.15
A.bias | 0 | 0.018 | 0.011 | 0.009 | 0.006 | 0.005 | 0.005 | 0.004
0.01 | 0.015 | 0.004 | 0.004 | 0.002 | 0.003 | 0.004 | 0.004
1 | 0.019 | 0.179 | 0.048 | 0.022 | 0.007 | 0.004 | 0.004
ASE | 0 | 0.288 | 0.120 | 0.093 | 0.076 | 0.066 | 0.060 | 0.054
0.01 | 0.287 | 0.120 | 0.092 | 0.076 | 0.066 | 0.060 | 0.054
1 | 0.287 | 0.121 | 0.093 | 0.076 | 0.066 | 0.060 | 0.054
ESE | 0 | 0.295 | 0.139 | 0.096 | 0.075 | 0.065 | 0.059 | 0.054
0.01 | 0.290 | 0.134 | 0.098 | 0.075 | 0.066 | 0.060 | 0.056
1 | 0.306 | 0.161 | 0.104 | 0.079 | 0.067 | 0.059 | 0.055
CP | 0 | 0.934 | 0.904 | 0.937 | 0.952 | 0.953 | 0.951 | 0.950
0.01 | 0.945 | 0.925 | 0.945 | 0.958 | 0.956 | 0.941 | 0.946
1 | 0.931 | 0.655 | 0.899 | 0.918 | 0.945 | 0.958 | 0.955
ACL | 0 | 1.127 | 0.472 | 0.363 | 0.299 | 0.261 | 0.233 | 0.213
0.01 | 1.126 | 0.472 | 0.362 | 0.298 | 0.261 | 0.234 | 0.213
1 | 1.127 | 0.474 | 0.363 | 0.299 | 0.261 | 0.233 | 0.213
Figure 3: QQ plots of standardized $\widehat{\beta}_{\text{on},r}^{(b)}$ with total sample size $N_{b}=420$, $p=400$ and $\bm{\Sigma}=\bm{I}_{p}$. Each column represents the estimated parameter $\widehat{\beta}_{\text{on},r}^{(b)}$ at data batches $b=2,6,10$. Each row corresponds to a true value of parameter $\bm{\beta}_{0}$, i.e. $\beta_{0,r}=0,0.01,1$.  Table 3:  $N_{b}=1200$, $b=12$, $p=1000$, $s_{0}=20$, $\bm{\Sigma}=\bm{I}_{p}$. Performance on statistical inference. Tuning parameter $\lambda$ is chosen from $T_{\lambda}=\\{0.15,0.20,0.25,0.30\\}$ using the adaptive method in Section 2.4. Simulation results are summarized over 200 replications. In the table, we report the $\lambda$ selected with highest frequency among 200 replications. | $\beta_{0,r}$ | OLS | online debiased lasso
---|---|---|---
data batch index |  |  | 2 | 4 | 6 | 8 | 10 | 12
$\lambda$ |  |  | 0.25 | 0.15 | 0.15 | 0.15 | 0.15 | 0.15
A.bias | 0 | 0.004 | 0.005 | 0.003 | 0.002 | 0.002 | 0.002 | 0.001
0.01 | 0.004 | 0.005 | 0.004 | 0.002 | 0.002 | 0.001 | 0.001
1 | 0.003 | 0.019 | 0.006 | 0.004 | 0.003 | 0.003 | 0.002
ASE | 0 | 0.071 | 0.083 | 0.056 | 0.045 | 0.039 | 0.035 | 0.032
0.01 | 0.071 | 0.083 | 0.056 | 0.045 | 0.039 | 0.035 | 0.032
1 | 0.071 | 0.083 | 0.056 | 0.045 | 0.039 | 0.035 | 0.032
ESE | 0 | 0.071 | 0.088 | 0.055 | 0.045 | 0.039 | 0.035 | 0.032
0.01 | 0.072 | 0.087 | 0.054 | 0.046 | 0.039 | 0.036 | 0.033
1 | 0.072 | 0.095 | 0.056 | 0.046 | 0.039 | 0.035 | 0.032
CP | 0 | 0.947 | 0.933 | 0.956 | 0.952 | 0.951 | 0.950 | 0.950
0.01 | 0.946 | 0.939 | 0.961 | 0.947 | 0.948 | 0.946 | 0.946
1 | 0.947 | 0.906 | 0.941 | 0.940 | 0.950 | 0.950 | 0.953
ACL | 0 | 0.277 | 0.325 | 0.220 | 0.178 | 0.154 | 0.137 | 0.125
0.01 | 0.277 | 0.325 | 0.220 | 0.178 | 0.154 | 0.137 | 0.125
1 | 0.278 | 0.325 | 0.220 | 0.178 | 0.154 | 0.137 | 0.125
Table 4:  $N_{b}=1200$, $b=12$, $p=1000$, $s_{0}=20$, $\bm{\Sigma}=\\{0.5^{|i-j|}\\}_{i,j=1,\dots,p}$. Performance on statistical inference. Tuning parameter $\lambda$ is chosen from $T_{\lambda}=\\{0.15,0.20,0.25,0.30\\}$ using the adaptive method in Section 2.4. Simulation results are summarized over 200 replications. In the table, we report the $\lambda$ selected with highest frequency among 200 replications. | $\beta_{0,r}$ | OLS | online debiased lasso
---|---|---|---
data batch index |  |  | 2 | 4 | 6 | 8 | 10 | 12
$\lambda$ |  |  | 0.30 | 0.25 | 0.15 | 0.15 | 0.15 | 0.15
A.bias | 0 | 0.005 | 0.007 | 0.004 | 0.004 | 0.004 | 0.003 | 0.003
0.01 | 0.005 | 0.006 | 0.003 | 0.003 | 0.003 | 0.002 | 0.002
1 | 0.005 | 0.017 | 0.005 | 0.002 | 0.002 | 0.002 | 0.001
ASE | 0 | 0.091 | 0.087 | 0.060 | 0.049 | 0.043 | 0.038 | 0.035
0.01 | 0.091 | 0.086 | 0.060 | 0.049 | 0.043 | 0.038 | 0.035
1 | 0.091 | 0.087 | 0.060 | 0.049 | 0.043 | 0.038 | 0.035
ESE | 0 | 0.092 | 0.091 | 0.059 | 0.049 | 0.043 | 0.038 | 0.035
0.01 | 0.090 | 0.091 | 0.059 | 0.049 | 0.043 | 0.039 | 0.035
1 | 0.090 | 0.095 | 0.060 | 0.049 | 0.042 | 0.038 | 0.034
CP | 0 | 0.946 | 0.935 | 0.953 | 0.949 | 0.948 | 0.947 | 0.946
0.01 | 0.956 | 0.934 | 0.952 | 0.949 | 0.946 | 0.947 | 0.948
1 | 0.949 | 0.921 | 0.948 | 0.950 | 0.956 | 0.950 | 0.958
ACL | 0 | 0.358 | 0.339 | 0.236 | 0.193 | 0.168 | 0.150 | 0.137
0.01 | 0.358 | 0.339 | 0.236 | 0.193 | 0.168 | 0.150 | 0.137
1 | 0.357 | 0.339 | 0.236 | 0.193 | 0.168 | 0.150 | 0.137

##  5 Applications

###  5.1 Analysis of Beijing PM 2.5 data

We apply the proposed ODL to analyze the Beijing PM2.5 Data by Liang et al. ,
(2015), which is available in UC Irvine Machine Learning Repository. Fine
particulate matter less than 2.5 microns (PM2.5) is an air pollutant that
threatens human health. Therefore, understanding the changes of PM 2.5 level
is an important issue. The dataset contains hourly PM2.5 records from 1
January 2010 to 31 December 2014, in conjunction with 5 meteorological
features. We are interested in whether the meteorological variables, such as
the wind direction, have an influence on PM 2.5 level.

Before applying our online algorithm, we first preprocess the original raw
data. We transform the categorical predictors into the one-hot vector. We also
include interaction terms, which are coded as products of all pairs of the
original features. As a result, the dimension of the feature vector is
$p=296$. Since the curve of an exponential distribution fits the PM 2.5 data
well, we use the logarithm of PM 2.5 as the response variable. In addition, we
split the data into $b=120$ batches fairly by its chronological order. Each
batch contains half-month data with size $n_{j}=348,j=1,\ldots,b$.

First, we examine the influence of wind direction. There are 4 types of wind
directions: northwest (NW), northeast (NE), southeast (SE), and calm and
variable (cv). The results are shown in Figure 4. From the left panel, we can
observe that in the most cases, SE wind has a positive influence on PM 2.5
while NE and NW has a negative impact. This observation is consistent with the
statement in Liang et al. , (2015). The major heavily polluting industries are
located at the south and east of Beijing, but the north region lacks
industries of this kind. Besides, another interesting observation is that
comparing to other seasons, all wind directions are not significant on
decreasing the level of PM 2.5 in winter and the SE wind even has positive
effect on increasing the level of PM 2.5. One possible explanation is the
heating supply in northern China in winter. At that time, coals are burned to
provide the heat which significantly increases the PM 2.5 level in the whole
region. In the presence of coal burning, wind directions are insignificant
variables. In the middle panel, we present the estimated standard errors. As
expected, the standard errors decrease as the number of batches increases.
Combining the results in the left and middle panels, we present the
$t$-statistic on the right panel.

Figure 4:  The influence of wind direction in different seasons. Left panel:
the heat map of the estimated coefficients at the end of a year. Middle panel:
the corresponding estimated standard errors. Right panel: the $t$-statistic.

Next, we focus on another two variables: pressure and dew point. The results
are shown in Figure 5. It can be seen that the increase of dew point is
associated with the increase of PM 2.5 except in the summer. On the contrary,
apart from the summer, the pressure itself has a negative impact on the PM
2.5. This finding also agrees with the study in Liang et al. , (2015). The
main difference is on the influence of dew point in summer time. We believe
the difference arises from the interaction terms. Actually, the coefficient of
the square of the dew point is significantly positive, and its estimated
standard error is similar to the middle panel in Figure 4. Both of them have a
decreasing trend. The values of $t$-statistic are also presented on the right
panel. For ease of illustration, we also present the trace of the outcomes on
the wind direction, pressure and dew point to show the trend of the estimation
with the influx of new data. The result is presented in Figure 6. Moreover, we
identify other significant variables such as the wind speed and some
interaction terms in this analysis. These findings suggest some interesting
covariates that warrant further investigation and validation.

Figure 5:  The influence of pressure and dew point in different seasons. Left
panel: the heat map of the estimated coefficients at the end of a year. Middle
panel: the corresponding estimated standard errors. Right panel: the
$t$-statistic. Figure 6:  The trace plots on the influence of wind direction,
pressure and dew point in different seasons. Each vertical bar corresponds to
the 95% confidence interval. Left panel: wind direction in Summer and Winter.
Middle panel: wind direction in Spring and Autumn. Right panel: pressure and
dew point in four seasons.

###  5.2 Hang Seng Index fund data

We next illustrate the application of ODL with an index fund dataset. This
dataset consists of the returns of $1148$ stocks listed in Hong Kong Stock
Exchange and the Hang Seng Index (HSI, a freefloat-adjusted market-
capitalization-weighted stock-market index in Hong Kong) during the period
from January 2010 to December 2020. The response variable is the return of the
HSI for every three days, and the predictors are every-three-day returns of
the $1148$ stocks. We partition the data into batches according to
chronological order. Specifically, the first batch consists of a two-year
dataset from 2010-2011 to ensure sufficient sample size at the initial stage
and each subsequent batch contains one-year data. Hence, $b=10,n_{1}=164$ and
$n_{j}=82$ for $j=2,\ldots,10.$ Similar to Lan et al. , (2016), the goal of
this study is to identify the most relevant stocks that can be used to create
a portfolio for index tracking.

The proposed ODL method is applied and the coefficients of 19 stocks are
identified to be significant at a significance level $\alpha=0.05.$ The
estimates of the $19$ regression coefficients and their standard errors are
presented in Figure 7. Among these selected stocks, only three of them (with
stocks code 0004.HK, 1088.HK and 3988.HK) are not constituent stocks of the
current HSI (June 2021), but they are highly associated with the constituent
stocks of the HSI. For example, 0004.HK (Wharf Holdings) is the parent company
of 1997.HK (Wharf Real Estate Investment Company Ltd), a constituent stock of
the current HSI. For the other 16 selected stocks, they cover all sub-indexes
of HSI, including Finance Sub-index, Utilities Sub-index, Properties Sub-index
and Commerce & Industry Sub-index. Our analysis demonstrates the importance of
selecting diversified stocks to establish a portfolio for tracking HSI.

In Figure 7, the left panel displays the estimated coefficients of the 19
stocks, among which the significance of many stocks does not change much in
past years except for 0700.HK (Tencent Holding Ltd). Tencent was listed on the
Hong Kong Stock Exchange in 2004 and was added as a Hang Seng Index
Constituent Stock in 2008. The Chinese tech giant has become the most valuable
publicly traded company in China in 2018 and thus its weight in HSI was
increasing in past years, which is consistent with our analysis. In the right
panel, as expected, the standard errors decrease as more and more data are
collected. However, the standard errors of several stocks increase in
2019-2020. We believe this might be related to the impact of the unprecedented
COVID-19 pandemic. COVID-19 virus has ravaged economies all over the world and
changed consumer behavior and preferences. In the COVID-19 pandemic, there are
many losers in traditional industries, but the tech giants are thriving, as
demand for online services and digital utilities has exploded. The shift in
market may have caused extra uncertainty in statistical analysis. In summary,
our proposed ODL performs reasonably well in analyzing this financial dataset.

Figure 7:  Analysis results of the Hang Seng Index fund data. Left panel: the
heat map of the estimated coefficients. Right panel: the corresponding
estimated standard errors. Note that the results in the first column of each
graph is based on the data collected from 2010-2011.

##  6 Discussion

In this paper we developed an online debiased lasso estimator for statistical
inference in linear models with high-dimensional streaming data. The proposed
method does not assume the availability of the full dataset at the initial
stage and only requires the availability of the current batch of the data
stream and the sufficient statistics of the historical data. A natural dynamic
tuning parameter selection procedure that takes advantage of streaming data
structure is developed as an important ingredient of the proposed algorithm.
The proposed online inference procedure is justified theoretically under
regularity conditions similar to those in the offline setting and mild
conditions on the batch size.

There are several other interesting questions that deserve further study.
First, we focused on the problem of making statistical inference about
individual regression coefficients, the proposed method can be extended to the
case of making inference about a fixed and low-dimensional subvector of the
coefficient. Second, we did not address the problem of variable selection in
the online learning setting consider here. This is apparently different from
the variable selection problem in the offline setting. The main issue is how
to recover the variables that are dropped at the early stages of the stream
but may be important as more data come in. Third, it would be interesting to
generalize the proposed method to generalized linear and nonlinear models.
These questions warrant thorough investigation in the future.

## Acknowledgment

The work of J. Huang is partially supported by the U.S. NSF grant DMS-1916199.
The work of Y. Lin is supported by the Hong Kong Research Grants Council
(Grant No. 14306219 and 14306620), the National Natural Science Foundation of
China (Grant No. 11961028) and Direct Grants for Research, The Chinese
University of Hong Kong.

## References

  * Akter & Wamba, (2016) Akter, S., & Wamba, S. F. 2016.  Big data analytics in E-commerce: a systematic review and agenda for future research.  Electronic Markets, 26, 173–194.
  * Bühlmann & van de Geer, (2011) Bühlmann, P., & van de Geer, S. 2011.  Statistics for high-dimensional data: methods, theory and applications.  Berlin: Springer.
  * Carbone et al. , (2015) Carbone, P., Katsifodimos, A., Ewen, S., Markl, V., Haridi, S., & Tzoumas, K. 2015\.  Apache Flink: stream and batch processing in a single engine.  IEEE Data Engineering Bulletin, 38, 28–38.
  * Cardot & Degras, (2018) Cardot, H., & Degras, D. 2018.  Online Principal Component Analysis in High Dimension: Which Algorithm to Choose?  International Statistical Review, 86, 29 – 50.
  * Choi et al. , (2016) Choi, J., Cho, Y., Shim, E., & Woo, H. 2016.  Web-based infectious disease surveillance systems and public health perspectives: a systematic review.  BMC Public Health, 16(1238).
  * Das et al. , (2018) Das, S., Behera, R. K., & Rath, S. K. 2018.  Real-Time Sentiment Analysis of Twitter Streaming data for Stock Prediction.  Procedia Computer Science, 132, 956–964.  International Conference on Computational Intelligence and Data Science.
  * Daubechies et al. , (2004) Daubechies, I., Defrise, M., & De Mol, C. 2004.  An iterative thresholding algorithm for linear inverse problems with a sparsity constraint.  Communications on Pure and Applied Mathematics, 57(11), 1413–1457.
  * Deshpande et al. , (2019) Deshpande, Y., Javanmard, A., & Mehrabi, M. 2019.  Online Debiasing for Adaptively Collected High-dimensional Data.  arXiv preprint arXiv:1911.01040.
  * Dezeure et al. , (2015) Dezeure, R., Bühlmann, P., Meier, L., & Meinshausen, N. 2015.  High-dimensional inference: confidence intervals, $p$-values and R-software hdi.  Statistical Science, 30(4), 533–558.
  * Donoho & Johnstone, (1994) Donoho, D. L., & Johnstone, J. M. 1994.  Ideal spatial adaptation by wavelet shrinkage.  Biometrika, 81, 425–455.
  * Duchi et al. , (2011) Duchi, J., Hazan, E., & Singer, Y. 2011.  Adaptive subgradient methods for online learning and stochastic optimization.  The Journal of Machine Learning Research, 12, 2121–2159.
  * Fan et al. , (2018) Fan, J., Gong, W., Li, C. J., & Sun, Q. 2018.  Statistical sparse online regression: a diffusion approximation perspective.  Pages 1017–1026 of: Storkey, Amos, & Perez-Cruz, Fernando (eds), Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics.  Proceedings of Machine Learning Research, vol. 84.  PMLR.
  * Friedman et al. , (2007) Friedman, J., Hastie, T., Höfling, H., & Tibshirani, R. 2007.  Pathwise coordinate optimization.  The Annals of Applied Statistics, 1(2), 302–332.
  * Javanmard & Montanari, (2014) Javanmard, A., & Montanari, A. 2014.  Confidence intervals and hypothesis testing for high-dimensional regression.  The Journal of Machine Learning Research, 15(1), 2869–2909.
  * Jiang et al. , (2018) Jiang, T., Yang, J., Yu, C., & Sang, Y. 2018.  A clickstream data analysis of the differences between visiting behaviors of desktop and mobile users.  Data and Information Management, 2(3), 130–140.
  * Kraft et al. , (2020) Kraft, R., Birk, F., Reichert, M., Deshpande, A., Schlee, W., Langguth, B., Baumeister, H., Probst, T., Spiliopoulou, M., & Pryss, R. 2020.  Efficient processing of geospatial mHealth data using a scalable crowdsensing platform.  Sensors, 20(12), 3456.
  * Lan et al. , (2016) Lan, W., Zhong, P. S., Li, R., Wang, H., & Tsai, C. L. 2016.  Testing a single regression coefficient in high dimensional linear models.  Journal of Econometrics, 195(1), 154–168.
  * Langford et al. , (2009) Langford, J., Li, L., & Zhang, T. 2009.  Sparse online learning via truncated gradient.  Journal of Machine Learning Research, 10, 777–801.
  * Liang et al. , (2015) Liang, X., Zou, T., Guo, B., Li, S., Zhang, H., Zhang, S., Huang, H., & Chen, S. 2015.  Assessing Beijing’s PM2. 5 pollution: severity, weather impact, APEC and winter heating.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 471(2182), 2015–0257.
  * Luo & Song, (2020) Luo, L., & Song, P. X.-K. 2020.  Renewable Estimation and Incremental Inference in Generalized Linear Models with Streaming Datasets.  Journal of the Royal Statistical Society: Series B, 82, 69–97.
  * Raskutti et al. , (2010) Raskutti, G., Wainwright, M. J., & Yu, B. 2010.  Restricted eigenvalue properties for correlated Gaussian designs.  The Journal of Machine Learning Research, 11, 2241–2259.
  * Samaras et al. , (2020) Samaras, L., García-Barriocanal, E., & Sicilia, M. A. 2020.  Syndromic surveillance using web data: a systematic review.  Innovation in Health Informatics, 39–77.
  * Schifano et al. , (2016) Schifano, E. D., Wu, J., Wang, C., Yan, J., & Chen, M. H. 2016.  Online updating of statistical inference in the big data setting.  Technometrics, 58(3), 393–403.
  * Shameer et al. , (2017) Shameer, K., Badgeley, M.A., Miotto, R., Glicksberg, B.S., Morgan, J.W., & Dudley, J.T. 2017.  Translational bioinformatics in the era of real-time biomedical, health care and wellness data streams.  Briefings in Bioinformatics, 18(1), 105–124.
  * Shi et al. , (2020) Shi, C., Song, R., Lu, W., & Li, R. 2020.  Statistical inference for high-dimensional models via recursive online-score estimation.  Journal of the American Statistical Association, 1–12.
  * Sun et al. , (2020) Sun, L., Wang, M., Guo, Y., & Barbu, A. 2020.  A novel framework for online supervised learning with feature selection.  arXiv preprint arXiv:1803.11521.
  * Tarrès & Yao, (2014) Tarrès, P., & Yao, Y. 2014.  Online Learning as Stochastic Approximation of Regularization Paths: Optimality and Almost-Sure Convergence.  IEEE Transactions on Information Theory, 60(9), 5716–5735.
  * Tashman, (2000) Tashman, L. J. 2000.  Out-of-sample tests of forecasting accuracy: an analysis and review.  International Journal of Forecasting, 16, 437–450.
  * Tibshirani, (1996) Tibshirani, R. 1996.  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society: Series B, 58(1), 267–288.
  * van de Geer et al. , (2014) van de Geer, S., Bühlmann, P., Ritov, Y. A., & Dezeure, R. 2014.  On asymptotically optimal confidence regions and tests for high-dimensional models.  The Annals of Statistics, 42(3), 1166–1202.
  * Zhang, (2010) Zhang, C. H. 2010.  Nearly unbiased variable selection under minimax concave penalty.  The Annals of Statistics, 38(2), 894–942.
  * Zhang & Zhang, (2014) Zhang, C. H., & Zhang, S. S. 2014.  Confidence intervals for low dimensional parameters in high dimensional linear models.  Journal of the Royal Statistical Society: Series B, 76(1), 217 – 242.
  * Zou & Hastie, (2005) Zou, H., & Hastie, T. 2005.  Regularization and variable selection via the elastic net.  Journal of the Royal Statistical Society: Series B, 67(2), 301–320.

##  Appendix A Appendix

In the appendix, we prove Theorems 1-3 and include an additional figure of Q-Q
plots from the simulation studies in Section 4.2. We first define the notation
needed below. For any sequences $\\{X_{N}\\}_{N\in\mathbb{N}}$ and
$\\{Y_{N}\\}_{N\in\mathbb{N}}$, we say that $X_{N}=O_{\mathbb{P}}(Y_{N})$ if
for any $\epsilon>0$, there exists $M_{1},M_{2}>0$ such that
$\mathbb{P}(|X_{n}/Y_{n}|\geq M_{1})<\epsilon$ for any $n>M_{2}$. Roughly
speaking, $X_{N}=O_{\mathbb{P}}(Y_{N})$ means that $X_{N}/Y_{N}$ is
stochastically bounded. Besides, $X_{N}=o_{\mathbb{P}}(Y_{N})$ means that
$X_{N}/Y_{N}$ converges to zero in probability. Particularly,
$X_{N}=\Omega(Y_{N})$ if $X_{N}=O_{\mathbb{P}}(Y_{N})$ and
$Y_{N}=O_{\mathbb{P}}(X_{N})$. In addition, we use $c_{1},c_{2},c_{3}\ldots$
to stand for the constants which do not depend on $N_{b}$.

Lemmas 1 \- 5 are needed to prove Theorem 1. The first two lemmas show the
consistency of the lasso estimators.

######  Lemma 1.

Suppose that Assumption 1 holds and $n_{1}\geq c_{1}s_{r}\log p$ for some
constant $c_{1}$. Then, for any $j=1,\ldots,b$, with probability at least
$1-p^{-3}$, low-dimensional projection defined in (10) with
$\lambda_{j}=c_{2}\sqrt{{\log p}/{N_{j}}}$ satisfies,

$\lVert\bm{\widehat{\gamma}}^{(j)}_{r}-\bm{{\gamma}}_{r}\lVert_{1}\leq
c_{3}s_{r}\lambda_{j}.$ (14)

###### Proof of Lemma 1.

For notational convenience, we suppress the subcript $r$ of $\bm{\gamma}$ in
this proof. For any fixed $j=1,\ldots,b,$

$\bm{\widehat{\gamma}}^{(j)}:=\underset{\bm{\gamma}\in\mathbb{R}^{(p-1)}}{\arg\min}\left\\{\frac{1}{2N_{j}}\sum_{i=1}^{j}\|\bm{x}^{(i)}_{r}-\bm{X}^{(i)}_{-r}\bm{\gamma}\|_{2}^{2}+\lambda_{j}\|\bm{\gamma}\|_{1}\right\\},$
(15)

where $N_{j}$ is the cumulative size of data at step $j$. Note that
$S_{r}=\\{k:\bm{\Theta}_{k,r}\neq 0,k\neq r\\}$ and $s_{r}=|S_{r}|.$ To
establish the consistency of the lasso estimator, we first show that
$\widetilde{\bm{\Sigma}}^{(j)}_{-r}:{=}\sum_{i=1}^{j}(\bm{X}^{(i)}_{-r})^{\top}\bm{X}^{(i)}_{-r}/N_{j}$
satisfies the compatibility condition for the set $S_{r}$. Namely, there is a
constant $\phi_{j}$ such that for all $\bm{\gamma}$ satisfying
$\|\bm{\gamma}_{S_{r}^{c}}\|_{1}\leq 3\|\bm{\gamma}_{S_{r}}\|_{1},$ it holds
that

$\|\bm{\gamma}_{S_{r}}\|_{1}^{2}\leq(\bm{\gamma}^{\top}\widetilde{\bm{\Sigma}}^{(j)}_{-r}\bm{\gamma})s_{r}/\phi_{j}^{2},$
(16)

where the $i$-th element of $\bm{\gamma}_{S_{r}}$ is denoted by
$\bm{\gamma}_{S_{r},i}=\bm{\gamma}_{i}1_{\\{i\in S_{r}\\}}$ for
$i=1,\ldots,{p-1}.$ It follows from an extension of Corollary 1 in Raskutti et
al. , (2010) from Gaussian case to sub-Gaussian case that, with probability at
least $1-p^{4}$, the above inequality holds as long as $N_{j}\geq
c_{1}s_{r}\log p$ and ${\bm{\Sigma}}^{(j)}_{-r}$ meets the compatibility
condition, which hold under the assumption that $n_{1}\geq c_{1}s_{r}\log p$
and (A2) in Assumption 1 respectively.

The remaining proof follows standard arguments. More notations are introduced.
Recall that

$\bm{\gamma}:=\underset{\bm{\gamma}\in\mathbb{R}^{(p-1)}}{\arg\min}\
\mathbb{E}\left\\{\frac{1}{2N_{j}}\sum_{i=1}^{j}\|\bm{x}^{(i)}_{r}-\bm{X}^{(i)}_{-r}\bm{\gamma}\|_{2}^{2}\right\\}$

and let $\bm{v}^{(j)}=\bm{\widehat{\gamma}}^{(j)}-\bm{{\gamma}}.$ Since
$\bm{\widehat{\gamma}}^{(j)}$ is the lasso estimator defined in (15), it
follows that

$\displaystyle\frac{1}{2N_{j}}\sum_{i=1}^{j}\|\bm{x}^{(i)}_{r}-\bm{X}^{(i)}_{-r}\bm{\widehat{\gamma}}^{(j)}\|_{2}^{2}+\lambda_{j}\|\bm{\widehat{\gamma}}^{(j)}\|_{1}\leq\frac{1}{2N_{j}}\sum_{i=1}^{j}\|\bm{x}^{(i)}_{r}-\bm{X}^{(i)}_{-r}\bm{{\gamma}}\|_{2}^{2}+\lambda_{j}\|\bm{{\gamma}}\|_{1}.$

Then,

$\displaystyle(\bm{v}^{(j)})^{\top}\widetilde{\bm{\Sigma}}^{(j)}_{-r}\bm{v}^{(j)}$
$\displaystyle\leq$
$\displaystyle\frac{2}{N_{j}}\sum_{i=1}^{j}(\bm{x}^{(i)}_{r}-\bm{X}^{(i)}_{-r}\bm{{\gamma}})^{\top}\bm{X}^{(i)}_{-r}\bm{v}^{(j)}+2\lambda_{j}(\|\bm{{\gamma}}\|_{1}-\|\bm{\widehat{\gamma}}^{(j)}\|_{1})$
$\displaystyle\leq$ $\displaystyle\frac{2}{N_{j}}\max_{k\neq
r}\left\lvert\sum_{i=1}^{j}(\bm{x}^{(i)}_{r}-\bm{X}^{(i)}_{-r}\bm{{\gamma}})^{\top}\bm{x}^{(i)}_{k}\right\lvert\|\bm{v}^{(j)}\|_{1}+2\lambda_{j}(\|\bm{{\gamma}}\|_{1}-\|\bm{\widehat{\gamma}}^{(j)}\|_{1}).$

Note that

$\displaystyle\sum_{i=1}^{j}(\bm{x}^{(i)}_{r}-\bm{X}^{(i)}_{-r}\bm{{\gamma}})^{\top}\bm{x}^{(i)}_{k}=\sum_{i=1}^{j}\sum_{l=1}^{n_{i}}({x}^{(i)}_{r,l}-\bm{X}^{(i)}_{-r,l}\bm{{\gamma}})^{\top}{x}^{(i)}_{k,l},$
(17)

where ${x}^{(i)}_{r,l}$ and $\bm{X}^{(i)}_{-r,l}$ are the explanatory
variables from the $l$-th observation in $i$-th data batch. Then, (17) is
written as the sum of i.i.d random variables. Specifically,
$\mathbb{E}\\{({x}^{(i)}_{r,l}-\bm{X}^{(i)}_{-r,l}\bm{{\gamma}})^{\top}{x}^{(i)}_{k,l}\\}=0$
and $({x}^{(i)}_{r,l}-\bm{X}^{(i)}_{-r,l}\bm{{\gamma}})^{\top}{x}^{(i)}_{k,l}$
is sub-exponential distributed by the definition of $\bm{\gamma}$ and (A1) in
Assumption 1 respectively. By Bernstein inequality, we obtain that

$\displaystyle\mathbb{P}\left(\left|\sum_{i=1}^{j}\sum_{l=1}^{n_{i}}({x}^{(i)}_{r,l}-\bm{X}^{(i)}_{-r,l}\bm{{\gamma}})^{\top}{x}^{(i)}_{k,l}\right|\geq
2c_{2}\sqrt{{N_{j}}{\log p}}\right)\leq p^{-5}$

for some constant $c$ which does not depend on $p$ and $N_{i}$. Since the
above inequality holds for any $k\neq r$, by Bonferroni inequality, it holds
that

$\displaystyle\mathbb{P}\left(\max_{k\neq
r}\left|\sum_{i=1}^{j}\sum_{l=1}^{n_{i}}({x}^{(i)}_{r,l}-\bm{X}^{(i)}_{-r,l}\bm{{\gamma}})^{\top}{x}^{(i)}_{k,l}\right|<2c_{2}\sqrt{{N_{j}}{\log
p}}\right)>1-p^{-4}.$

By choosing $\lambda_{j}=c_{2}\sqrt{{\log p}/{N_{j}}}$, then, with probability
at least $1-p^{-4},$

$\displaystyle(\bm{v}^{(j)})^{\top}\widetilde{\bm{\Sigma}}^{(j)}_{-r}\bm{v}^{(j)}$
$\displaystyle\leq$
$\displaystyle\lambda_{j}\|\bm{v}^{(j)}\|_{1}+2\lambda_{j}(\|\bm{{\gamma}}\|_{1}-\|\bm{\widehat{\gamma}}^{(j)}\|_{1})$
$\displaystyle=$
$\displaystyle\lambda_{j}\|\bm{v}^{(j)}\|_{1}+2\lambda_{j}(\|\bm{{\gamma}}_{S_{r}^{c}}\|_{1}-\|\bm{\widehat{\gamma}}^{(j)}_{S_{r}}\|_{1}-\|\bm{\widehat{\gamma}}^{(j)}_{S_{r}^{c}}\|_{1})$
$\displaystyle\leq$
$\displaystyle\lambda_{j}\|\bm{v}^{(j)}\|_{1}+2\lambda_{j}(\|\bm{v}^{(j)}_{S_{r}}\|_{1}-\|\bm{v}^{(j)}_{S_{r}^{c}}\|_{1})$
$\displaystyle=$
$\displaystyle\lambda_{j}(3\|\bm{v}^{(j)}_{S_{r}}\|_{1}-\|\bm{v}^{(j)}_{S_{r}^{c}}\|_{1}),$

where (A) holds due to $\|\bm{{\gamma}}_{S_{r}}\|_{1}=0.$ It further implies
that $\|\bm{v}^{(j)}_{S_{r}^{c}}\|_{1}\leq 3\|\bm{v}^{(j)}_{S_{r}}\|_{1}.$
Together with the compatibility condition, we have

$\displaystyle\frac{\|\bm{v}_{S_{r}}^{(j)}\|_{1}^{2}\phi_{j}^{2}}{s_{r}}\leq(\bm{v}^{(j)})^{\top}\widetilde{\bm{\Sigma}}^{(j)}_{-r}\bm{v}^{(j)}\leq
3\lambda_{j}\|\bm{v}^{(j)}_{S_{r}}\|_{1}.$

Consequently,

$\displaystyle\|\bm{v}^{(j)}\|_{1}\leq 4\|\bm{v}_{S_{r}}^{(j)}\|_{1}\leq
12s_{r}\lambda_{j}/\phi_{j}^{2}.$ (19)

Since (19) holds for any $j=1,\ldots,b$, the proof of Lemma 1 is complete.
$\Box$

######  Lemma 2.

Suppose that Assumption 1 hold and $N_{b}\geq c_{1}s_{0}\log p$ for some
constant $c_{1}$. Then, with probability at least $1-p^{-4}$, the lasso
estimator in (5) with $\lambda_{b}=c_{2}\sqrt{{\log p}/{N_{b}}}$ satisfies,

$\lVert{\widehat{\bm{\beta}}}^{(b)}-{{\bm{\beta}}}_{0}\lVert_{1}\leq
c_{3}s_{0}\lambda_{b}.$

The proof of Lemma 2 is structurally similar to the proof of Lemma 1 by
letting $j=b$. $(A3)$ in Assumption 1 is used to obtain the concentration
inequality as Bernstein inequality in Lemma 1. We omit the details here. The
next lemma is used to estimate the cumulative terms in the online learning.

######  Lemma 3.

Recall that $n_{j}$ and $N_{j}$ are the batch size and the cumulative batch
size respectively when the $j$-th data arrives, $j=1,\ldots,b.$ Then,

$\displaystyle\sum_{j=1}^{b}\frac{n_{j}}{N_{j}}$ $\displaystyle\leq$
$\displaystyle 1+\log\frac{N_{b}}{n_{1}},$ (20)
$\displaystyle\sum_{j=1}^{b}\frac{n_{j}}{\sqrt{N_{j}}}$ $\displaystyle\leq$
$\displaystyle 2\sqrt{N_{b}}.$ (21)

###### Proof of Lemma 3.

We first prove (20). Let

$f(t)=\log(1+t)-\frac{t}{1+t},\ t>0.$

Since $f(0)=0$ and $f^{\prime}(t)>0$ for $t>0$, we have $f(t)\geq 0.$ Choosing
$t=n_{b}/N_{b-1}$ yields
$\log\left({N_{b}}/{N_{b-1}}\right)\geq{n_{b}}/{N_{b}},$ namely,
$\log\left({N_{b}}\right)\geq\log\left({N_{b-1}}\right)+n_{b}/N_{b}.$ Repeat
the above procedure by letting $t=n_{j}/N_{j-1}$ for $j=b-1,\ldots,2.$ It then
follows that

$\displaystyle\log\left({N_{b}}\right)$ $\displaystyle\geq$
$\displaystyle\log\left({N_{b-1}}\right)+\frac{n_{b}}{N_{b}}$
$\displaystyle\geq$
$\displaystyle\log\left({N_{b-2}}\right)+\frac{n_{b-1}}{N_{b-1}}+\frac{n_{b}}{N_{b}}$
$\displaystyle\cdots$ $\displaystyle\geq$
$\displaystyle\log(N_{1})+\sum_{j=2}^{b}\frac{n_{j}}{N_{j}}.$

Then, in view of $n_{1}=N_{1}$, (20) holds. The remaining step is to prove
(21). We claim that

$2(\sqrt{a+b}-\sqrt{a})\geq\frac{b}{\sqrt{a+b}},\ \ \text{for}\ a,b>0.$

Let $a=N_{b-1}$ and $b=n_{b}$. We have $2\sqrt{N_{b}}\geq
2\sqrt{N_{b-1}}+n_{b}/\sqrt{N_{b}}.$ Similarly, we repeat this procedure by
choosing $a=N_{j-1},b=n_{j}$ for $j=b-1,\ldots,2.$ It follows that

$\displaystyle 2\sqrt{N_{b}}$ $\displaystyle\geq$ $\displaystyle
2\sqrt{N_{b-1}}+\frac{n_{b}}{\sqrt{N_{b}}}$ $\displaystyle\geq$ $\displaystyle
2\sqrt{N_{b-2}}+\frac{n_{b-1}}{\sqrt{N_{b-1}}}+\frac{n_{b}}{\sqrt{N_{b}}}$
$\displaystyle\cdots$ $\displaystyle\geq$ $\displaystyle
2\sqrt{N_{1}}+\sum_{j=2}^{b}\frac{n_{j}}{\sqrt{N_{j}}}.$

Given $n_{1}=N_{1}$, (21) holds. $\Box$

For any matrix $\bm{A}=(a_{ij})$, let $\|\bm{A}\|_{\infty}$ be the largest
absolute value of its elements, that is,
$\|\bm{A}\|_{\infty}=\max_{i,j}|a_{ij}|$. The next two lemmas give the bound
for the error term $\Delta_{j}$ in Theorem 1.

######  Lemma 4.

Suppose that the conditions in Lemma 1 holds and the subsequent batch size
$n_{j}\geq c\log p,j=2,\ldots,b$, for some constants $c$. If

$s_{r}^{2}\frac{\log p}{N_{b}}\log{\frac{N_{b}}{n_{1}}}=o(1),$ (22)

then, $\lVert\bm{\widehat{z}}_{r}\lVert_{2}=\Omega(\sqrt{N_{b}}).$

###### Proof of Lemma 4.

By the triangle inequality,

$\displaystyle\lVert\bm{\widehat{z}}_{r}\lVert_{2}$ $\displaystyle\leq$
$\displaystyle\lVert\bm{{z}}_{r}\lVert_{2}+\lVert\bm{\widehat{z}}_{r}-\bm{{z}}_{r}\lVert_{2}=\lVert\bm{{z}}_{r}\lVert_{2}+\sqrt{\sum_{j=1}^{b}\lVert\bm{\widehat{z}}^{(j)}_{r}-\bm{{z}}^{(j)}_{r}\lVert_{2}^{2}},$

where
$\bm{z}_{r}=((\bm{z}_{r}^{(1)})^{\top},\ldots,(\bm{z}_{r}^{(b)})^{\top})^{\top}$.

First, we intend to show that
$\lVert\bm{{z}}_{r}\lVert_{2}^{2}=\Omega(N_{b}).$ Recall that
${z}_{r,1}^{(1)}$, ${x}_{r,1}^{(1)}$ and $\bm{X}_{-r,1}^{(1)}$ denote the
first element of ${\bm{z}}_{r}^{(1)}$, ${\bm{x}}_{r}^{(1)}$ and the first row
of $\bm{X}_{-r}^{(1)}$ respectively. Consider
$\zeta_{r}:=\mathbb{E}\\{({z}_{r,1}^{(1)})^{2}\\}=\mathbb{E}\\{({x}_{r,1}^{(1)}-\bm{X}_{-r,1}^{(1)}\bm{{\gamma}}_{r})^{2}\\}.$
Under Assumption 2, $\Lambda_{\min}^{2}\leq\zeta_{r}\leq\Sigma_{j,j}=O(1).$ It
then follows from the law of large numbers that

$\lVert\bm{{z}}_{r}\lVert_{2}^{2}=\mathbb{E}(\lVert\bm{{z}}_{r}\lVert_{2}^{2})+O_{\mathbb{P}}(\sqrt{N_{b}})=\zeta_{r}N_{b}+O_{\mathbb{P}}(\sqrt{N_{b}}).$

Next, we demonstrate that
$\sum_{j=1}^{b}\lVert\bm{\widehat{z}}^{(j)}_{r}-\bm{{z}}^{(j)}_{r}\lVert_{2}^{2}=o_{\mathbb{P}}(N_{b})$.
According to Lemma 1,

$\displaystyle\sum_{j=1}^{b}\lVert\bm{\widehat{z}}^{(j)}_{r}-\bm{{z}}^{(j)}_{r}\lVert_{2}^{2}$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{b}\lVert\bm{X}^{(j)}_{-r}(\bm{\widehat{\gamma}}^{(j)}_{r}-\bm{{\gamma}}_{r})\lVert_{2}^{2}$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{b}n_{j}|(\bm{\widehat{\gamma}}^{(j)}_{r}-\bm{{\gamma}}_{r})^{\top}\widehat{\bm{\Sigma}}^{(j)}_{-r}(\bm{\widehat{\gamma}}^{(j)}_{r}-\bm{{\gamma}}_{r})|$
$\displaystyle\leq$
$\displaystyle\sum_{j=1}^{b}n_{j}\|\widehat{\bm{\Sigma}}^{(j)}_{-r}\|_{\infty}\|\bm{\widehat{\gamma}}^{(j)}_{r}-\bm{{\gamma}}_{r}\|_{1}^{2},$

where
$\widehat{\bm{\Sigma}}^{(j)}_{-r}=(\bm{X}^{(j)}_{-r})^{\top}\bm{X}^{(j)}_{-r}/n_{j}.$
Recall that ${\bm{\Sigma}}_{-r}$ is the principle submatrix of ${\bm{\Sigma}}$
by removing the $r$-th row and the $r$-th column. It can be shown along
similar lines of the proof of Lemma 1 that, with probability at least
$1-p^{-4}$,

$\|\widehat{\bm{\Sigma}}^{(j)}_{-r}-\bm{\Sigma}_{-r}\|_{\infty}\leq
c_{4}\sqrt{\frac{\log p}{n_{j}}},\ \ \text{for }j=1,\ldots,b.$

Since $n_{j}\geq c\log p,j=1,\ldots,b$ and $\|\bm{\Sigma}_{-r}\|_{\infty}$ is
bounded, it follows that for some constant $c_{5}$,

$\|\widehat{\bm{\Sigma}}^{(j)}_{-r}\|_{\infty}\leq\|\bm{\Sigma}_{-r}\|_{\infty}+\|\widehat{\bm{\Sigma}}^{(j)}_{-r}-\bm{\Sigma}_{-r}\|_{\infty}\leq
c_{5},\text{for }j=1,\ldots,b.$

Consequently,

$\displaystyle\sum_{j=1}^{b}\|\bm{\widehat{z}}^{(j)}_{r}-\bm{{z}}^{(j)}_{r}\|_{2}^{2}$
$\displaystyle\leq$ $\displaystyle
c_{5}\sum_{j=1}^{b}n_{j}\|\bm{\widehat{\gamma}}^{(j)}_{r}-\bm{{\gamma}}_{r}\|_{1}^{2}$
$\displaystyle\leq$ $\displaystyle c_{5}s_{r}^{2}\log
p\sum_{j=1}^{b}\frac{n_{j}}{N_{j}}$ $\displaystyle\leq$ $\displaystyle
c_{5}{s_{r}^{2}\log p}\left(1+\log{\frac{N_{b}}{n_{1}}}\right),$

where the last inequality is from (20) in Lemma 3. As

$s_{r}^{2}\frac{\log p}{N_{b}}\log{\frac{N_{b}}{n_{1}}}\to 0\text{ as
}N_{b}\to\infty,$

then

$\sum_{j=1}^{b}\lVert\bm{\widehat{z}}^{(j)}_{r}-\bm{{z}}^{(j)}_{r}\lVert_{2}^{2}=o_{\mathbb{P}}(N_{b}).$

As a result, we have shown that $\lVert\bm{\widehat{z}}_{r}\lVert_{2}^{2}\leq
c_{6}N_{b}$ for some constant $c_{6}$ in probability. Similarly, in view of
the fact that

$\lVert\bm{\widehat{z}}_{r}\lVert_{2}\geq\lVert\bm{{z}}_{r}\lVert_{2}-\sqrt{\sum_{j=1}^{b}\lVert\bm{\widehat{z}}^{(j)}_{r}-\bm{{z}}^{(j)}_{r}\lVert_{2}^{2}},$

we can conclude that $\lVert\bm{\widehat{z}}_{r}\lVert_{2}^{2}\geq c_{7}N_{b}$
for some constant $c_{7}.$ We complete the proof of Lemma 4. $\Box$

######  Lemma 5.

Suppose that the conditions in Lemma 1 hold and the subsequent batch size
$n_{j}\geq c\log p,j=2,\ldots,b$, for some constants $c$. If

$s_{0}s_{r}\frac{\log(p)}{\sqrt{N_{b}}}=o(1),$

Then, $\Big{\lvert}\sum_{k\neq
r}\bm{\widehat{z}}_{r}^{\top}\bm{x}_{k}({\widehat{\beta}}^{(b)}_{k}-{\beta}_{0,k})\Big{\lvert}=O_{\mathbb{P}}\left(s_{0}s_{r}\log(p)\right)=o_{\mathbb{P}}\left(\sqrt{N_{b}}\right).$

###### Proof of Lemma 5.

As mentioned earlier in Remark 2, due to
$\bm{\widetilde{z}}^{(j)}_{r}=\bm{x}_{r}^{(j)}-\bm{X}_{-r}^{(j)}\bm{\widehat{\gamma}}^{(b)}_{r}$,
we cannot directly apply KKT condition here. To see the difference with the
proof in the offline debiased lasso, we first separate $\sum_{k\neq
r}\bm{\widehat{z}}_{r}^{\top}\bm{x}_{k}({\widehat{\beta}}^{(b)}_{k}-{\beta}_{0,k})$
into two parts: the offline term and one additional term from the online
algorithm. The upper bound of the former is derived from KKT condition while
the latter is tackled differently. Consider
$\bm{\widetilde{z}}_{r}=((\bm{\widetilde{z}}^{(1)}_{r})^{\top},\ldots,(\bm{\widetilde{z}}^{(b)}_{r})^{\top})^{\top}\in\mathbb{R}^{N_{b}}$
where
$\bm{\widetilde{z}}^{(j)}_{r}=\bm{x}^{(j)}_{r}-\bm{X}^{(j)}_{-r}\bm{\widehat{\gamma}}^{(b)}_{r}.$
Write

$\displaystyle\sum_{k\neq
r}\bm{\widehat{z}}_{r}^{\top}\bm{x}_{k}({\widehat{\beta}_{k}}^{(b)}-{\beta}_{0,k})$
$\displaystyle=$ $\displaystyle\sum_{k\neq
r}\bm{\widetilde{z}}_{r}^{\top}\bm{x}_{k}({\widehat{\beta}}^{(b)}_{k}-{\beta}_{0,k})+\sum_{k\neq
r}(\bm{\widehat{z}}_{r}-\bm{\widetilde{z}}_{r})^{\top}\bm{x}_{k}({\widehat{\beta}}^{(b)}_{k}-{\beta}_{0,k})$
$\displaystyle:{=}$ $\displaystyle\Pi_{\text{off}}+\Pi_{\text{on}},$

where $\Pi_{\text{off}}$ pertains to an offline term and $\Pi_{\text{on}}$
pertains to an online error. First,

$\Pi_{\text{off}}\leq\|\bm{\widehat{\beta}}^{(b)}-\bm{\beta}_{0}\|_{1}\max_{k\neq
r}|\bm{\widetilde{z}}_{r}^{\top}\bm{x}_{k}|.$

By the Karush-Kuhn-Tucker (KKT) condition and Lemma 2,

$\max_{k\neq r}|\bm{\widetilde{z}}_{r}^{\top}\bm{x}_{k}|\leq
N_{b}\lambda_{b}=c_{2}\sqrt{N_{b}{\log(p)}},\
\|\bm{\widehat{\beta}}^{(b)}-\bm{\beta}_{0}\|_{1}\leq
c_{3}s_{0}\sqrt{\frac{\log(p)}{N_{b}}}.$

Then,

$\Pi_{\text{off}}=O_{\mathbb{P}}\left(s_{0}\sqrt{\frac{\log(p)}{N_{b}}}\left\\{\sqrt{N_{b}{\log(p)}}\right\\}\right)=O_{\mathbb{P}}(s_{0}{\log(p)}).$

Second, for the online error,

$\displaystyle\Pi_{\text{on}}$ $\displaystyle\leq$
$\displaystyle\|\bm{\widehat{\beta}}^{(b)}-\bm{\beta}_{0}\|_{1}\max_{k\neq
r}|(\bm{\widehat{z}}_{r}-\bm{\widetilde{z}}_{r})^{\top}\bm{x}_{k}|$
$\displaystyle=$
$\displaystyle\|\bm{\widehat{\beta}}^{(b)}-\bm{\beta}_{0}\|_{1}\max_{k\neq
r}\Big{\lvert}\sum_{j=1}^{b}(\bm{\widehat{z}}^{(j)}_{r}-\bm{\widetilde{z}}^{(j)}_{r})^{\top}\bm{x}^{(j)}_{k}\Big{\lvert}.$

By Lemma 1, we obtain

$\displaystyle|(\bm{\widehat{z}}^{(j)}_{r}-\bm{\widetilde{z}}^{(j)}_{r})^{\top}\bm{x}^{(j)}_{k}|$
$\displaystyle\leq$
$\displaystyle\lVert\bm{\widehat{\gamma}}^{(j)}_{r}-\bm{\widehat{\gamma}}^{(b)}_{r}\lVert_{1}\max_{m\neq
r}|[\bm{x}^{(j)}_{m}]^{\top}\bm{x}^{(j)}_{k}|$ $\displaystyle=$ $\displaystyle
O_{\mathbb{P}}\left(s_{r}\sqrt{\frac{\log(p)}{N_{j}}}\times
n_{j}\|\widehat{\bm{\Sigma}}^{(j)}_{-r}\|_{\infty}\right).$

In the proof of Lemma 4, we have shown that
$\|\widehat{\bm{\Sigma}}^{(j)}_{-r}\|_{\infty}$ is bounded with probability
tending to 1. As a result,

$\displaystyle\Big{\lvert}\sum_{j=1}^{b}(\bm{\widehat{z}}^{(j)}_{r}-\bm{\widetilde{z}}^{(j)}_{r})^{\top}\bm{x}^{(j)}_{k}\Big{\lvert}$
$\displaystyle=$ $\displaystyle
O_{\mathbb{P}}\left(\sum_{j=1}^{b}s_{r}\sqrt{\frac{n_{j}^{2}\log
p}{N_{j}}}\right)$ $\displaystyle=$ $\displaystyle
O_{\mathbb{P}}\left(s_{r}\sqrt{N_{b}\log(p)}\right),$

where the last equation is from (21) in Lemma 3.

Then,

$\Pi_{\text{on}}=O_{\mathbb{P}}\left(s_{0}\sqrt{\frac{\log(p)}{N_{b}}}\left\\{s_{r}\sqrt{N_{b}\log(p)}\right\\}\right)=O_{\mathbb{P}}\left(s_{0}s_{r}\log(p)\right).$

Since $s_{0}s_{r}{\log(p)}/{\sqrt{N_{b}}}=o(1)$, the statement of the lemma
follows.

$\Box$

###### Proof of Theorem 1.

Recall that
$\bm{X}=((\bm{X}^{(1)})^{\top},\ldots,(\bm{X}^{(b)})^{\top})^{\top},\bm{y}=((\bm{y}^{(1)})^{\top},\ldots,(\bm{y}^{(b)})^{\top})^{\top}$,
$\bm{x}_{r}=((\bm{x}^{(1)}_{r})^{\top},\ldots,(\bm{x}_{r}^{(b)})^{\top})^{\top}$
and
$\bm{\widehat{z}}_{r}=((\bm{\widehat{z}}^{(1)}_{r})^{\top},\ldots,(\bm{\widehat{z}}_{r}^{(b)})^{\top})^{\top}$.
Then, we write the online debiased estimator in (11) into the following vector
form:

${\widehat{\beta}}^{(b)}_{\text{on},r}={\widehat{\beta}}^{(b)}_{r}-\frac{\bm{\widehat{z}}_{r}^{\top}(\bm{y}-\bm{X}\bm{\widehat{\beta}}^{(b)})}{\bm{\widehat{z}}_{r}^{\top}\bm{x}_{r}}.$

Subtract the true parameter ${\beta}_{0,r}$ and obtain

${\widehat{\beta}}^{(b)}_{\text{on},r}-{\beta}_{0,r}=\frac{\lVert\bm{\widehat{z}}_{r}\lVert_{2}}{\bm{\widehat{z}}_{r}^{\top}\bm{x}_{r}}\left\\{\frac{\bm{\widehat{z}}_{r}^{\top}\bm{\epsilon}}{\lVert\bm{\widehat{z}}_{r}\lVert_{2}}-\frac{\sum_{k\neq
r}\bm{\widehat{z}}_{r}^{\top}\bm{x}_{k}({\widehat{\beta}}^{(b)}_{k}-{\beta}_{0,k})}{\lVert\bm{\widehat{z}}_{r}\lVert_{2}}\right\\}.$

What remains to be shown is that,

$\displaystyle\lVert\bm{\widehat{z}}_{r}\lVert_{2}={\Omega(\sqrt{N_{b}})},$
$\displaystyle\sum_{k\neq
r}\bm{\widehat{z}}_{r}^{\top}\bm{x}_{k}({\widehat{\beta}}^{(b)}_{k}-{\beta}_{0,k})=o_{\mathbb{P}}(\sqrt{N_{b}}).$

as detailed by Lemma 4 and Lemma 5 respectively. $\Box$

###### Proof of Theorem 2 and Theorem 3.

Theorem 2 and Theorem 3 can be proved in the same fashion as Theorem 1. Due to
limited space, we only point out the main difference. In the proof of Theorem
2, we will show the upper bound of
$\|\widehat{\bm{\Sigma}}^{(j)}_{-r}\|_{\infty}$ is $O_{\mathbb{P}}(\log p)$,
by replacing $n_{j},j=2,\ldots,b$ with $1$ in Lemma 4 and Lemma 5. For Theorem
3, the major difference is to establish a similar lemma to Lemma 1 with
conclusion
$\lVert\bm{\widehat{\gamma}}^{(j)}_{r}-\bm{{\gamma}}_{r}\lVert_{1}\leq
cs_{r}^{2}\lambda_{j}$ under Assumption 2.

$\Box$

Figure 8: QQ plots of standardized $\widehat{\beta}_{\text{on},r}^{(b)}$ with
total sample size $N_{b}=1200$, $p=1000$ and
$\bm{\Sigma}=\\{0.5^{|i-j|}\\}_{i,j=1,\dots,p}$. Each column represents the
estimated parameter $\widehat{\beta}_{\text{on},r}^{(b)}$ at data batches
$b=2,6,10$. Each row corresponds to a true value of parameter
$\bm{\beta}_{0}$, i.e. $\beta_{0,r}=0,0.01,1$.