\section{Introduction}

Spin and orbital angular momentum  (SAM and OAM, respectively) are nearly independent quantities when light beams propagate paraxially in vacuum or through  isotropic and homogeneous media \cite{Int_1,Int_2,Int_3}, and are independently conserved during both paraxial and non-paraxial propagation \cite{Int_4}. However, propagation through inhomogeneous or anisotropic media \cite{Int_non_parax}, scattering processes \cite{Int_5}, as well as tight focusing in isotropic inhomogeneous media lead to interactions between SAM and OAM, giving rise to tangible effects in the mesoscopic scale with applications in nanosensing \cite{Int_5}, or particle manipulation \cite{the_3}.  In the case of tight focusing in optical tweezers, even a fundamental Gaussian beam evolves non-paraxially due to the presence of a spin-orbit and orbit-oribit interaction term in the expression of total angular momentum \cite{Int_5,Int_6,Int_7,Int_8,Int_9}. Now, both the spin-orbit interaction term (SOI) and orbit-orbit interaction terms are enhanced as a result of the generation of a large longitudinal component of the electric field due to tight focusing, and further accentuated by inserting a refractive index contrast in the path of the light beam that increases the geometric phase gradient of the focused light. On a different note, another interesting manifestation of the interplay between spin and orbital angular momentum arises in the spin and orbital Hall effect, which is essentially the transverse spatial separation of opposite angular momentum components induced by the spin-orbit and orbit-orbit interaction, preserving the angular momentum conservation law. This angular momentum dependent separation may be represented in terms of spin-spin separation called spin Hall effect \cite{Int_7,Int_10,Int_11,Int_12,Int_13,Int_14,Int_15}, orbit-orbit separation called orbital Hall effect \cite{Int_16,Int_17,Int_18}, and spin-orbit separation called spin-orbit Hall effect \cite{Int_19} - where the transverse spatial separation of the respective opposite spin and/or orbital angular momentum components is observed.

The spin and orbital Hall effect have been extensively studied in both fundamental Gaussian beams \cite{the_3}, and in orbital angular momentum carrying Laguerre-Gaussian (LG) beams \cite{the_9,Int_1}, with interesting effects in spin-orbit and orbit-orbit conversions in tightly focused LG beams leading to intriguing orbital motion of single optically trapped birefringent particles \cite{the_8}. However, effects of spin-orbit interaction in structured Gaussian beams carrying no intrinsic orbital angular momentum (typically called Hermite-Gaussian or HG beams) have largely been ignored in the literature.  A study on such beams, however, merits attention, since any HG mode may be written as a superposition of two LG modes having opposite topological charge. This may lead to very interesting effects of spin-orbit interaction, especially in the case of input spin-polarized HG beams. In addition, the breaking of radial symmetry in the intensity profile distribution of such beams may lead to direct visualization of SOI effects in the beam structure itself.

In this paper, we study this problem in detail, and demonstrate clear manifestations of azimuthal rotation of the transverse intensity profile of  a circularly polarized $HG_{10}$ beam after it is tightly focused using optical tweezers. We describe this as a rotational spin-Hall effect. Thus, we observe that for input right or left-circularly polarized light (RCP and LCP, respectively), tight focusing using a high numerical aperture (NA) objective lens in an optical tweezers setup, leads to the generation of the opposite spin component that is also coupled with a corresponding orbital angular momentum mode - with the final effect being the rotation of the transverse intensity profile of the output.  For the same reason, the transverse component of the longitudinal field intensity generated due to non-paraxial propagation of the light, also displays a rotation.  The magnitude of rotation for the respective components is determined by the diffraction integrals (or Debye-Wolf integrals) $ I_{11}, I_{12}$ and $I_{14}$ (for transverse field components), and $I_{10}$ and $I_{13}$ (for longitudinal field components) - which provide the extent of the spin to orbit conversion of angular momentum \cite{Int_5, the_4}. We also study the dependence of the rotation on the refractive index (RI) contrast in the path of the beam after it is tightly focused, and observe that the rotational spin Hall effect in general increases monotonically with increasing RI contrast. We verify this experimentally by projecting out the opposite spin polarized transverse intensity component from the input, and the longitudinal component in an optical tweezers configuration.

 We employ the Debye-Wolf theory or angular spectrum method \cite{the_2, the_4} to determine the electric field at the output of the high NA objective lens for an input spin-polarized HG beam, also considering RI stratification of the media through which the beam travels after focusing \cite{the_1,the_2,the_3,the_4,the_5}. The expression of the output electric field from the input electric field may be written as:
\onecolumngrid
\begin{equation}
\left[\begin{array}{c}
E_{x}^{o} \\
E_{y}^{o} \\
E_{z}^{o}
\end{array}\right]=\left[\begin{array}{ccc}
i I_{11} \cos \psi+i I_{14} \cos 3 \psi & -i I_{12} \sin \psi+i I_{14} \sin 3 \psi & 2 I_{10}-2 I_{13} \cos 2 \psi \\
-i I_{12} \sin \psi+i I_{14} \sin 3 \psi & i\left(I_{11}+2 I_{12}\right) \cos \psi-i I_{14} \cos 3 \psi & -2 I_{13} \sin 2 \psi \\
-2 I_{10}+2 I_{13} \cos 2 \psi & 2 I_{13} \sin 2 \psi & i\left(I_{11}-I_{12}\right) \cos \psi
\end{array}\right] \times\left[\begin{array}{c}
E_{x}^{i} \\
E_{y}^{i} \\
E_{z}^{i}
\end{array}\right]
\label{eq1}
\end{equation}
\twocolumngrid
Where $E_{o}$ and $E_{i}$ denote the output and input electric fields related through the 3x3 Jones matrix, respectively, and $I_{10}, I_{11}, I_{12},  I_{14} $ and $I_{13} $  are the Debye-Wolf integrals. Now, given that the Jones vectors for $x$ and $y$-polarized input light are $\mathrm{E}_{\mathrm{x}}^i=\left[\begin{array}{lll}
1 & 0 & 0
\end{array}\right]^{T} \text { and } \mathrm{E}_{\mathrm{y}}^i=\left[\begin{array}{lll}
0 & 1 & 0
\end{array}\right]^{T}$, \\
we have from Eq.~\ref{eq1}, the output electric field for $HG_{10}$ input $x$-polarized light as
\begin{equation}
\left[\begin{array}{c}
E_{x}^{o} \\
E_{y}^{o} \\
E_{z}^{o}
\end{array}\right]_{x-pol}=\left[\begin{array}{c}
i I_{11} \cos \psi+i I_{14} \cos 3 \psi \\
-i I_{12} \sin \psi+i I_{14} \sin 3 \psi \\
-2 I_{10}+2 I_{13} \cos 2 \psi
\end{array}\right]
\label{eq2}
\end{equation}. It is important to note that the  circular polarized $HG_{10}$ beam does not carry OAM, but does possess an SAM of magnitude ±$\hslash$, so that the total angular momentum equals the SAM. Now, noting that the Jones vectors for input RCP and LCP light are $E_{r c p / l c p}=\left[\begin{array}{lll}
1 & \pm i & 0
\end{array}\right]^{T}
$, we use Eq.~\ref{eq1}, to determine the output electric field in both cases as
\begin{equation}
{\left[E_{x}^{0}~E_{y}^{0}~E_{z}^{0}\right]_{R C P / L C P}^{T}=\left[\begin{array}{lll}
a & b & c
\end{array}\right]^{\top}},
\label{eq3}
\end{equation}
where, $a=i I_{11} \cos \psi\pm I_{12} \sin \psi + I_{14} (i \cos 3\psi \mp \sin 3 \psi)$, $b=-i I_{12} \sin \psi+i I_{14} \sin 3 \psi \mp\left(I_{11}+2 I_{12}\right) \cos \psi \pm I_{14} \cos 3 \psi$, and $c=-2 I_{0}+2 I_{13} \cos 2 \psi \pm 2 i I_{13} \sin 2 \psi$.

\begin{center}
\begin{figure}[H]
\includegraphics[width=0.47\textwidth]{Fig_1.png}
\caption{Simulation of rotation of the total intensity profile at $z=2~\mu~m$ away the focus of the high NA objective (trap focus) for an input $HG_{10}$ beam mode for matched RI ((a), (b), (c)) and RI 1.814 ((d), (e) and  (f)). (a) and (d) are for input $x$-polarized light, while  (b) and (e) are for input RCP, (c) and (f) are for input LCP light, (g) rotation due to transverse components for RI 1.814, (h) rotation due to transverse component of longitudinal intensity profile for input RCP}
\label{ssphase0}
\end{figure}
\end{center}
We then decompose the above output electric fields for input RCP (LCP) light in terms of the SAM and OAM components - note that the effects of SOI cause all helicity components to be present along with the corresponding OAM modes in the output electric field.
\onecolumngrid
\begin{equation}
\left[\begin{array}{c}
E_{x}^{0} \\
E_{y}^{0} \\
E_{z}^{0}
\end{array}\right]_{R/L}=\frac{i}{2}\left(I_{11} e^{\pm i \psi}+I_{11} e^{\mp i \psi}+I_{12} e^{\pm i \psi}+I_{12} e^{\mp i \psi}\right)\\ \times \left[\begin{array}{c}
1 \\
\pm i \\
0
\end{array}\right] +i\left(I_{14} e^{\pm 3 i \psi}-I_{12} e^{\pm i \psi}\right)\left[\begin{array}{c}
1 \\
\mp i \\
0
\end{array}\right]\\+\left (2 I_{13}e^{\pm 2 i \psi }-2 I_{10}\right)\left[\begin{array}{c}
0 \\
0 \\
1
\end{array}\right]
\label{eq4}
\end{equation}
\twocolumngrid
It is clear that the first four terms are associated with the same helicity as the input light, coupled correspondingly with appropriate positive and negative OAM, respectively, while SOI leads to the fifth and sixth terms containing opposite helicity with $ l=\pm \hslash$ and $l=\pm 3\hslash $, respectively, satisfying the conservation of total angular momentum. The longitudinal component of the field with $l=\pm 2\hslash $ and $l=0$ appear in the last two terms. Understandably, the effects of SOI would be extracted by projecting out the field components having the opposite helicity as the input field, as well as the longitudinal component. This is manifested in the corresponding output intensity profile obtained from Eq.~\ref{eq4} - where, for an input RCP/LCP $HG_{10}$ mode, both the transverse (opposite spin to the input) and longitudinal intensity profiles display an azimuthal rotation about the input beam axis in a clockwise/anticlockwise direction (Fig.~\ref{ssphase0}). We observe these effects in our numerical simulations that we describe next.

We have simulated our experimental system, which consists of a stratified medium in the path of the tightly focused light after the output of the high NA objective in an optical tweezers configuration, and attempted to observe the rotational spin Hall-effect of light. As is well known, the linearly polarized light can be decomposed into RCP and LCP components in the transverse plane. Thus, the laser beam of wavelength 671 nm is incident on the 100× oil-immersion objective of NA 1.4 followed by (a) an oil layer of thickness around 5 $\mu m$ and refractive index (RI) 1.516, (b) a 160 $\mu m$ thick coverslip having refractive indices 1.516, 1.572, 1.695 and 1.814 (note that the case where the RI = 1.516 is henceforth referred to as the ``matched condition'', which is typically employed in optical tweezers to minimize spherical aberration effects in the focused beam spot.) (c) a water layer having a refractive index of 1.33 with a depth of 35 $\mu m$, and finally (chosen since probe particles in future experiments will be immersed in water) (d) a glass slide of refractive index 1.516 whose thickness we consider to be semi-infinite (~1500 $\mu m$). In the simulation, the origin of coordinates is taken inside the water layer at an axial distance of 5 $\mu m$ from the interface between water and the cover-slip. Thus, the objective-oil interface is at -170 $\mu m$, oil-coverslip interface is at -165 $\mu m$, cover-slip-sample chamber interface is at -5 $\mu m$, and sample chamber-glass slide interface is at +30 $\mu m$. On tight-focusing of the Hermite-Gaussian beam linearly polarized in the transverse plane, we observe that the electric field in the focal plane exhibits a component not only along the incident polarization direction, but also in the orthogonal and longitudinal directions. As we showed previously, the electric field in the transverse plane depends on Debye-Wolf integrals $I_{11}, I_{12},$ and $I_{14}$, while the longitudinal component depends on $I_{10}$ and $I_{13}$.

In order to extract the effects of SOI which lead to the rotational spin Hall effect, we plot the Debye-Wolf components associated with the polarization which is orthogonal to the input state, viz. RCP for input LCP, and vice versa. Thus, in Fig.~\ref{ssphase0}, we plot the transverse components ($I_{12}$ and $I_{14}$) of orthogonal polarization to the input state and longitudinal components of the total intensity (from Eq.~\ref{eq2}, \ref{eq4}) at $z=2~\mu$m away from the beam focus, and clearly observe rotations of the intensity profiles relative to the axis of the input modal distribution (rotational spin-Hall shift). We use both the matched RI for the cover slip (Fig.~\ref{ssphase0}(a)-(c)), and a mismatched RI of 1.814 (Fig.~\ref{ssphase0}(d)-(f)). While the rotational spin Hall shifts for input linear polarization ((Fig.~\ref{ssphase0}(a), (d)) is zero (since the shifts for the constituent RCP and LCP components cancel out), they are significantly high for both RCP ((Fig.~\ref{ssphase0}(b), (e)), and LCP (Fig.~\ref{ssphase0}(c), (f)). Note that the non-paraxial propagation of light also leads to a transverse component of the the longitudinal field intensity that also displays a rotation. We compare the rotational effects of the different field intensity components in Fig.~\ref{ssphase0}(g)-(h), where (g) shows the spin rotational shift due to the transverse components (LCP), and (h) shows the spin rotational shift due to the transverse component of the longitudinal intensity for input RCP light, and (e) the total (which term we now use to denote transverse LCP + transverse component of the longitudinal field intensities) shifts for input RCP light. We measure the rotations with respect to the axis of the linearly polarized beam mode that does not undergo rotation after focusing, and display the measured angles of rotation for an input RCP $HG_{10}$ beam as a function of RI from simulations in Fig.~\ref{rotmeas}. We observe a maximum of ~7.3 degrees azimuthal rotation for the transverse components for an RI of 1.814, while that for the transverse component of the longitudinal intensity profile is ~4.5 degrees. The total intensity - a combination of the rotations of both field intensity components - profile is observed to rotate by ~4.7 degrees azimuthally. We also observe that an increase in RI contrast tends to increase the observed azimuthal rotation. This can be understood from the fact that both the geometric phase acquired by the circularly polarized light during propagation, and its gradient, are enhanced with increasing RI contrast. Note that it is the geometric phase that is responsible for the SOI, while its gradient is at the heart of the observed spin-dependent azimuthal rotation of the intensity profile \cite{the_7}.  In addition, since the contribution of the transverse component of the longitudinal field intensity is ~90\% or even more of the total intensity (see Fig.~2 of Supplementary Information), it is the rotation of this component that finally dominates in the rotation of the total intensity profile, which we finally measure experimentally. We now proceed to describe our experiments to verify these simulations.
\begin{center}
\begin{figure}[!h]
\includegraphics[width=0.35\textwidth]{Fig_2.png}
\caption{ (a) Measurement of the rotational spin Hall effect for simulated transverse  (red solid circles), simulated longitudinal (blue solid circles), simulated total intensity (transverse + longitudinal, green solid circles)), and experimentally measured (transverse + longitudinal, black solid circles)) as a function of RI of the cover slip.}
\label{rotmeas}
\end{figure}
\end{center}

We use a conventional optical tweezers configuration consisting of an inverted microscope (Carl Zeiss Axioert.A1) with an oil-immersion 100× objective (Zeiss, NA 1.4) and a He-Ne laser (671 $nm$, 200 $mW$) coupled into the back port of the microscope. We use a vortex half-wave retarder of zero-order for generating a vector beam (i.e. LG beam). This beam is then passed through a linear polarizer to convert the LG mode into a $HG_{10}$ mode. Finally, we pass the beam through a quarter wave plate (QWP) centred at 671 $nm$ and oriented at 45 degrees to the input beam axis in order to circularly polarize the $HG_{10}$ mode. We then couple the circularly polarized HG beam into the microscope so that it is tightly focused into the stratified medium inserted beyond the objective lens. The stratified medium - as described earlier - consists of immersion oil (RI 1.516), cover slips (1.516, 1.572, 1.695 and 1.814 - we use individual cover slips for individual experiments), water (RI 1.33),  and a top glass slide (RI 1.516).

The back reflected light, from the cover slip passes through another QWP in order to project out the intensity of the spin state we would like to detect. Thus, for  incident RCP light, the output QWP filters out the RCP  component - so we detect the intensity only of the LCP component (superposed with the longitudinal component) on a CCD that we use to image the intensity profile. Importantly, tight focusing of linearly polarized HG beam does not show any rotation of the axis of $HG_{10}$ mode which we verify first. In order to measure rotation of the output mode, we perform the experiments first with linear polarized light that display no rotation, and can thus be employed as reference to measure the rotation for input spin polarized states. To determine the angle of rotation,  we use the CCD image of the mode to find out the distribution of pixels along the mode central region which have  minimum intensity, and fit a straight line connecting these pixels. This then defines the singularity axis of the mode (see Fig.~\ref{expresults}). We perform this analysis for all output modes for the input linear and spin polarized states, and measure the angle of rotation of the axis from the slope with respect to a vertical axis. We display the results of the rotation of the beam axes in Fig.~\ref{expresults}(a)-(f). Fig.~\ref{expresults}(a)-(c) display the rotational spin Hall shift for the matched RI, whereas (d)-(f) are for an RI of 1.814. Also, Fig.~\ref{expresults}(a) and (d) represent the reference linear polarization states with respect to which we measure the rotational spin Hall effect for the matched condition (Fig.~\ref{expresults} (b) and (c) - input RCP and LCP, respectively), and RI of 1.814 (Fig.~\ref{expresults} (e) and (f) - input RCP and LCP, respectively). Clearly, the beam axis appear tilted in opposite directions for opposite input spin states, with the tilt higher for the mismatched condition. The actual measured values of rotation are shown in Fig.~\ref{rotmeas} (black solid dots) - and we obtain a very good match with the values obtained from simulation (green solid dots). The errors in determining the rotation angles are between 5-20\% for individual cases.
\begin{figure}
    \centering
    \includegraphics[width = 0.35\textwidth]{Fig_3.png}
    \caption{Experimental measurements of the rotational Hall shift for RI matched and mismatched cover slips. (a)-(c) are for matched conditions (RI 1.516) - with (a) CCD image of the output intensity profile for input linearly polarized state which we use as a reference to measure  rotations for input spin polarized states, (b) output intensity profile for input RCP, (c) LCP. (d)-(f) are for RI 1.814, with (d) reference linearly polarized ouput state, (e) output intensity profile for input RCP, (f) LCP. The white straight lines through the beam axes show the linear fit for lowest intensity pixels along the beam axis which we determine from a least squared fit in MATLAB.}
    \label{expresults}
\end{figure}

In conclusion, we study the rotational spin Hall effect in a structured Gaussian beam ($HG_{01})$ mode) carrying no intrinsic OAM, when it is tightly focused through a high NA objective lens in an optical tweezers configuration with a stratified medium in the path of the focused light. For input circularly polarized states, the SOI effects due to tight focusing cause the generation of opposite circular polarization coupled with a corresponding OAM state to satisfy conservation of total AM. As a result, the intensity profile corresponding to these states undergoes a rotation with respect to the axis of the input state, allowing us to visualize the effects of SOI in the output intensity profile itself. In addition, the transverse component of the longitudinal field intensity profile also undergoes a rotation due to coupling with OAM states due to SOI. We quantify such rotations using the complete vector diffraction theory described by Debye and Wolf by numerical simulations of the output electric field after tight focusing through a stratified medium, and proceed to verify our simulations via careful experiments, where we project out the opposite circular polarization state along with the longitudinal component for a given input circular polarization. We obtain excellent matches with simulation, and thus clearly demonstrate an interesting manifestation of SOI for higher order Gaussian modes of light. The rotational spin Hall shift increases with RI contrast of the stratified medium, which is expected since the magnitude of the SOI increases correspondingly due to the enhanced geometric phase gradient in these cases. We would also like to point out that, while these studies were carried out on the spatial intensity profiles of structured beams in an optical tweezers configuration, we would like to conduct future experiments on the observation of rotation of optically trapped birefringent particles around the beam axis due to the generation of OAM by SOI in such beams. These would thus lead to interesting routes of complex particle manipulation using optical tweezers, in which direction in our research is currently engaged.

The authors acknowledge the SERB, Department of Science and Technology, Government of India (Project No. EMR/2017/001456) and IISER Kolkata IPh.D fellowship for research. They also acknowledge Anand Dev Ranjan and Sauvik Roy for their help in experiments and simulation.

\bibliographystyle{apsrev4-2}

\providecommand{\noopsort}[1]{}\providecommand{\singleletter}[1]{#1}%

\section{THEORETICAL CALCULATIONS}
The angular spectrum for the focal-field can be expressed as a function of the field at the focus as \cite{the_4}:
\begin{eqnarray}\label{polarint}
\vec{E}(\rho,\psi,z)&=&  i\frac{kfe^{-ikf}}{2\pi}\int_0^{\theta_{max}}\int_0^{2\pi}
\vec{E}_{res}(\theta,\phi)e^{ikz\cos\theta}\nonumber\\& \times & e^{ik\rho\sin\theta\cos(\phi-\psi)}
sin(\theta)\>d\theta d\phi,
\end{eqnarray}
where, $f$ is the focal length of the lens, and $\theta_{max}$ is because of the finite size of the aperture which is decided by the numerical aperture (NA) of the objective lens. Note that here we cannot use the paraxial approximation for calculating the electric field at the focal plane, because we are using a high NA objective lens (NA=1.4) for tight focusing. Thus, we use the angular spectrum method, which are the exact solutions for the non-paraxial regime to determine the nature of electric field distribution of Hermite Gaussian modes propagating in a medium or across stratified medium. This approach operates in the frequency domain, calculates the Fourier transform (FT) of the input field and multiplies the result with a transfer function, so that the desired output field is obtained by the inverse FT. The transfer function is given by $A=R_{z}(\phi) T R_{y}(\theta) R_{z}(-\phi)$ where $R_{z}$ and $R_{y}$ are $\mathrm{SO}(3)$ rotation matrices \cite{the_5}. Since the stratification of any medium makes the field propagating in the medium dependent on the input polarization, we incorporate the Fresnel transmission coefficients $\left(T_{s}\right.$
and $T_{p}$ ) as well as the Fresnel reflection coefficients $\left(R_{s}\right.$ and $\left.R_{p}\right)$ considering the contributions from both the $s$ and $p$ polarization. The output
and the incident field are related through a transfer function, A as:
$ \mathbf{E}_{r e s}(\theta, \phi)=A \mathbf{E}_{i n c}(\theta, \phi),$ where, the $T$ and $R$ matrices are given by:
\begin{equation}
T=\left(\begin{array}{ccc}
T_{p} & 0 & 0 \\
0 & T_{s} & 0 \\
0 & 0 & T_{p}
\end{array}\right) ; R=\left(\begin{array}{ccc}
-R_{p} & 0 & 0 \\
0 & R_{s} & 0 \\
0 & 0 & -R_{p}
\end{array}\right),
\end{equation}

Note that we have $ T_{i}^{(1, j)}=\dfrac{E_{i+}^{j}}{E_{i+}^{1}} ; R_{i}^{(1, j)}=\dfrac{E_{i-}^{j}}{E_{i+}^{1}};$. Here, $i$ specifies the polarization $(s$ and $p),~ +/-$ signifies a wave propagating forward and back-ward, respectively, and $j$ in the superscript specifies the layer of the stratified medium in which the optical tweezers (trapping laser) focus lies.

For an input $HG_{10}$ beam \cite{the_4}, we have
\begin{equation*}
\begin{array}{c}
E_{\text {inc }}(\theta, \phi)=E_{0}\left(2 x_{\infty} / w_{0}\right) \exp \left\{-\frac{x_{\infty}^{2}+y_{\infty}^{2}}{w_{0}^{2}}\right\}~~~~~~~~~~~~~~\\=\left(2 E_{0} f / w_{0}\right) \sin \theta \cos \phi e^{-f^{2} \sin ^{2} \theta / w_{0}^{2}}
\end{array}
\end{equation*}
From these equations, we can write the output electric field from the input electric field in the form of a matrix equation as:
\vspace*{1 cm}
\onecolumngrid
\begin{equation}
\left[\begin{array}{c}
E_{x}^{o} \\
E_{y}^{o} \\
E_{z}^{o}
\end{array}\right]=\left[\begin{array}{ccc}
i I_{11} \cos \psi+i I_{14} \cos 3 \psi & -i I_{12} \sin \psi+i I_{14} \sin 3 \psi & 2 I_{10}-2 I_{13} \cos 2 \psi \\
-i I_{12} \sin \psi+i I_{14} \sin 3 \psi & i\left(I_{11}+2 I_{12}\right) \cos \psi-i I_{14} \cos 3 \psi & -2 I_{13} \sin 2 \psi \\
-2 I_{10}+2 I_{13} \cos 2 \psi & 2 I_{13} \sin 2 \psi & i\left(I_{11}-I_{12}\right) \cos \psi
\end{array}\right] \times\left[\begin{array}{c}
E_{x}^{i} \\
E_{y}^{i} \\
E_{z}^{i}
\end{array}\right],
\end{equation}
where $E_{o}$ and $E_{i}$ denote the output and input Jones polarization vectors, respectively, and $I_{10}, I_{11}, I_{12},  I_{14} $ and $I_{13} $  are the Debye–Wolf integrals for the transmitted and reflected waves, which are given as \cite{the_1, the_2}
$$
\begin{array}{c}
I_{10}^{t}(\rho)=\int_{0}^{\theta_{\max }} E_{i n c}(\theta) \sqrt{\cos \theta} \sin ^{2} \theta T_{p} \sin \theta_{j} J_{0}\left(k_{1} \rho \sin \theta\right) e^{i k_{j} z \cos \theta_{j}} d \theta ,\\\\
I_{11}^{t}(\rho)=\int_{0}^{\theta_{\max }} E_{i n c}(\theta) \sqrt{\cos \theta} \sin ^{2} \theta\left(T_{s}+3 T_{p} \cos \theta_{j}\right) J_{1}\left(k_{1} \rho \sin \theta\right) e^{i k_{j} z \cos \theta_{j}} d \theta,\\\\
I_{12}^{t}(\rho)=\int_{0}^{\theta_{\max }} E_{i n c}(\theta) \sqrt{\cos \theta} \sin ^{2} \theta\left(T_{s}-T_{p} \cos \theta_{j}\right) J_{1}\left(k_{1} \rho \sin \theta\right) e^{i k_{j} z \cos \theta_{j}} d \theta,\\\\
I_{13}^{t}(\rho)=\int_{0}^{\theta_{\max }} E_{i n c}(\theta) \sqrt{\cos \theta} \sin ^{2} \theta T_{p} \sin \theta_{j} J_{2}\left(k_{1} \rho \sin \theta\right) e^{i k_{j} z \cos \theta_{j}} d \theta,\\\\
I_{14}^{t}(\rho)=\int_{0}^{\theta_{\max }} E_{i n c}(\theta) \sqrt{\cos \theta} \sin ^{2} \theta\left(T_{s}-T_{p} \cos \theta_{j}\right) J_{3}\left(k_{1} \rho \sin \theta\right) e^{i k_{j} z \cos \theta_{j}} d \theta,\\\\
I_{11}^{r}(\rho)=\int_{0}^{\theta_{\max }} E_{\text {inc }}(\theta) \sqrt{\cos \theta} \sin ^{2} \theta\left(R_{s}-3 R_{p} \cos \theta_{j}\right) J_{1}\left(k_{1} \rho \sin \theta\right) e^{-i k_{j} z \cos \theta_{j}} d \theta, \\\\
I_{12}^{r}(\rho)=\int_{0}^{\theta_{\max }} E_{\text {inc }}(\theta) \sqrt{\cos \theta} \sin ^{2} \theta\left(R_{s}-R_{p} \cos \theta_{j}\right) J_{1}\left(k_{1} \rho \sin \theta\right) e^{-i k_{j} z \cos \theta_{j}} d \theta, \\\\
I_{13}^{r}(\rho)=-\int_{0}^{\theta_{\max }} E_{\text {inc }}(\theta) \sqrt{\cos \theta} \sin ^{2} \theta R_{p} \sin \theta_{j} J_{2}\left(k_{1} \rho \sin \theta\right) e^{-i k_{j} z \cos \theta_{j}} d \theta, \\\\
I_{14}^{r}(\rho)=\int_{0}^{\theta_{\max }} E_{\text {inc }}(\theta) \sqrt{\cos \theta} \sin ^{2} \theta\left(R_{s}+R_{p} \cos \theta_{j}\right) J_{3}\left(k_{1} \rho \sin \theta\right) e^{-i k_{j} z \cos \theta_{j}} d \theta,
\end{array}
$$
Here, $t$ and $r$ superscripts denote the transmitted and reflected components, respectively, and $J_{n}$ is the Bessel function of the first kind of order $n$.\\

The radial intensity distribution of the output electric field with  input RCP and LCP HG beam is given by
\onecolumngrid
\begin{equation}
\begin{array}{c}
I_{\text {RCP/LCP }}(\rho)=2\left(\left|I_{12}\right|^{2}+\left|I_{14}\right|^{2}\right)+2\left(\left|I_{11}\right|^{2}+\left|I_{12}\right|^{2}\right) \cos ^{2} \psi+2\left(I_{11} I_{12}^{*}+I_{11}^{*} I_{12}\right) \cos ^{2} \psi \\
-2\left(I_{12} I_{14}^{*}+I_{12}^{*} I_{14}\right) \cos 2 \psi \pm 2 i\left(I_{12} I_{14}^{*}-I_{12}^{*} I_{14}\right) \sin 2 \psi
\end{array}
\label{radialeqn}
\end{equation}
The equations for the intensity profile are different for input RCP and LCP $HG_{10}$ beams. Therefore, the rotation of the intensity profile for input RCP and LCP beams are in clockwise and anti-clockwise directions, respectively.
\section{NUMERICAL SIMULATION}
We now run simulations on our experimental system (stratified medium in the path of the optical tweezers light beam) as described in the primary manuscript. According to Eq.~\ref{radialeqn}, the transverse intensity of RCP and LCP components are different, therefore intensity profile for RCP and LCP input $HG_{10}$ mode show clockwise and anti-clockwise rotation with respect to x-axis in transverse plane. We had provided results for the azimuthal rotation for RI 1.516 and 1.814 (Fig.~1 of primary manuscript) - here we provide those for RI 1.572 and 1.659. Once again, we keep only the transverse field intensities corresponding to the circular polarization opposite to the input, and the transverse component of the longitudinal field intensity profile. The simulation results are displayed in Fig.~\ref{ssphase0}(a)-(f). The measurement of the rotation angles are performed using the method described in the primary manuscript.
\begin{center}
\begin{figure}[]
\includegraphics[width=0.42\textwidth]{fig_4.png}
\caption{Simulation of rotation of the total intensity profile at $z=2~\mu~m$ away the focus of the high NA objective (trap focus) for an input $HG_{10}$ beam mode for mismatched RI 1.572 ((a), (b), (c)) and RI 1.695 ((d), (e) and  (f)). (a) and (d) are for input $x$-polarized light, while  (b) and (e) are for input RCP, (c) and (f) are for input LCP light.}
\label{ssphase0}
\end{figure}
\end{center}

Another useful exercise is to compare the transverse component of the longitudinal intensity profile ($E_{zt}$)to that of the transverse intensity component of the SOI-generated circular polarization that is opposite to the input polarization. This is displayed in Fig.~\ref{intensitycomp}. We observe that field corresponding to $E_{zt}$ (normalized) dominates over the transverse component of the radial intensity - with the former being almost 90\% of the sum of these two intensity components (refer to the total intensity defined in the text and in Fig.~1(g) of the primary manuscript) in the radial direction.  We also observe the contribution of the intensity profile of $E_{zt}$ to be increasing with increasing RI contrast of the stratified medium.  However, when we consider the total intensity in the radial direction, which also contains the intensity for the circular polarization in the direction of the input polarization, the $E_{zt}$ field intensity is only 20 to 30 \% of the total radial intensity. We also observe that this ratio reduces with increasing RI contrast, since the spread in the $k$ vector is reduced correspondingly, resulting in a lower z-component of the field compared to the total transverse component.
\begin{center}
\begin{figure}[]
\includegraphics[width=0.42\textwidth]{fig_5.png}
\caption{ (b) Percentage of longitudinal component of light compared to  spin-polarized transverse component induced by SOI (red solid circles), and total transverse intensity (blue solid circles).}
\label{intensitycomp}
\end{figure}
\end{center}
\begin{center}
\begin{figure}[!h!t]
\includegraphics[width=0.42\textwidth]{fig_6.png}
\caption{Experimental measurements of the rotational spin Hall shift for RI mismatched cover slips. (a)-(c) are RI 1.572 - with (a) CCD image of the output intensity profile for input linearly polarized state which we use as a reference to measure rotations for input spin polarized states,(b) output intensity profile for input RCP, (c) LCP. (d)-(f) are for R.I 1.695 with (d) reference linearly polarized output state, (e) output intensity profile for input RCP, (f) LCP. The white straight lines through the beam axes show the linear fit for lowest intensity pixels along the beam axis.}
\label{exprotri1572}
\end{figure}
\end{center}
\section{Experiment}
We have also performed experiments to determine the azimuthal rotation for an input RCP $HG_{10}$ beam for RI 1.572 and 1.695 of the cover slip in the stratified medium. The beam profiles are displayed in Fig.~\ref{exprotri1572}(a)-(c) for RI 1.572 and (d)-(f) for RI 1.695. The measured rotations agree well with the theoretically simulated values as shown in Fig. 2 of the primary manuscript.

\bibliographystyle{apsrev4-2}

\providecommand{\noopsort}[1]{}\providecommand{\singleletter}[1]{#1}%


\section{Introduction}
The advent of distributed online learning systems such as Apache Flink~\citep{Carbone2015ApacheFS} 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~\citep{Mobile_users_2018}, e-commerce purchases~\citep{ecommerce_2016}, infectious disease surveillance programs~\citep{disease2016,disease2020}, mobile health consortia~\citep{mHealth2017,mHealth2020}, and financial trading floors~\citep{finance2018}. 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 \citep
	truncated_SGD_Langford_2009,

	AdaGrad2011,
	Tarres2014,

	Sun2018ANF}.
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$~\citep{Schifano2016CUEE, Luo2020}. The goal of this paper is to develop an online algorithm and statistical inference procedure for analyzing high-dimensional streaming data.

\subsection{Related work}
The last decade has witnessed enormous progress on statistical inference in high-dimensional models, see, for example,
\citet{Zhang_delasso_2014,vandegeer2014,javanmard2014confidence}, as well as the review paper~\citet{hdi2015} 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~\citep{tibshirani1996regression}.
For example,~\citet{truncated_SGD_Langford_2009} proposed an online $\ell_1$-regularized method via a variant of the truncated SGD.~\citet{online_sparse_Fan_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 \citep{zou2005regularization}, Minimax Convex Penalty (MCP) \citep{zhang2010nearly}, and Feature Selection with Annealing (FSA) \citep{AdaGrad2011} has been studied by \citet{Sun2018ANF}.
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~\citep{Sun2018ANF}.
Recently, \citet{Deshpande2019OnlineDF} considered a class of online estimators in a high-dimensional auto-regressive model and studied the asymptotic properties via martingale theories.~\citet{shi2020statistical} 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.

\subsection{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.
\begin{itemize}
	\setlength\itemsep{-0.10 cm}
	\item 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.
	\item  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.
	\item 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.

	\item
	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.
\end{itemize}

The rest of the paper is organized as follows. Section~\ref{sec:online_debiased_lasso} presents the model formulation and our proposed ODL procedure. Section~\ref{sec3} includes the theoretical properties of the proposed ODL estimator. Simulation experiments are given in Section~\ref{sec:sim} to evaluate the performance of our proposed ODL in comparison to the offline ordinary least square estimator. In Section~\ref{sec:data} we demonstrate the application of the proposed method on an air pollution dataset. Concluding remarks are given in Section~\ref{sec:discussion}. Detailed proof of the theoretical properties is included in the appendix.

\section{Online debiased lasso}\label{sec:online_debiased_lasso}
Consider a time point $b\geq2$ 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:
\begin{equation}\label{lm1}
	\bm{y}^{(j)} = \bm{X}^{(j)}\bm{\beta}_0 + \bm{\epsilon}^{(j)}, \ j=1,\dots,b
\end{equation}
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 \ge 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 \citep{ Zhang_delasso_2014, vandegeer2014, javanmard2014confidence} 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~\eqref{lm1}.
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
\begin{equation}\label{lasso_step_1}
	\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\},
\end{equation}
where  $\lambda_1 \ge 0 $ is a regularization parameter.  Following \citet{Zhang_delasso_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
\begin{equation}\label{projection_step_1}
	\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\},
\end{equation}
and $\lambda_1$ is taken to be the same as in (\ref{lasso_step_1}) for simplicity.
Then,  the offline debiased lasso estimator of $ \bm{\beta}_{0}$, $r=1, \ldots,p$,  is defined as
\begin{equation}\label{debiased_step_1}
	{\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)}}.
\end{equation}
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~\eqref{lasso_step_1}-\eqref{debiased_step_1}.
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~\ref{sec: online lasso}-\ref{sec: online debias} below.

\subsection{Online lasso}\label{sec: online lasso}
Upon the arrival of data batch $\mathcal{D}_b=\{\bm{y}^{(b)},\bm{X}^{(b)} \}$ with $b\geq2$, 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:
\begin{equation}\label{lasso_step_j}
	\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\},
\end{equation}
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~\ref{sec: tuning parameter}.
However, since we only assume the availability of summary statistics of historical data, we cannot use the algorithms such as coordinate descent \citep{fhht2007} that requires the availability of the whole dataset.

Note that the objective function in~\eqref{lasso_step_j} depends on the data only through the following summary statistics:
\begin{equation}\label{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)}.
\end{equation}
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 \eqref{lasso_step_j} by the gradient descent algorithm.

Let $ \mathcal{L}(\boldsymbol{\beta}) = \sum_{j =1}^{b}\|\bm{y}^{(j)} - \bm{X}^{(j)}\bm{\beta} \|_2^2/(2N_b) $, and the gradient of $\mathcal{L}(\boldsymbol{\beta})$  is given by
\begin{equation}\label{gradient_descent}
	\frac{\partial \mathcal{L}(\boldsymbol{\beta}) }{\partial \boldsymbol{\beta}}  = \frac{1}{N_b}(\bm{{S}}^{(b)} \boldsymbol{\beta} - \bm{U}^{(b)}).
\end{equation}
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 \citep{daubechies2004iterative, donoho1994}. Specifically,
each iteration consists of the following two steps:
\begin{itemize}
	\item Step 1 (Gradient descent): update $\widehat{\boldsymbol{\beta}}^{(b)}$ through
	\begin{equation
		\widehat{\boldsymbol{\beta}}^{(b)} \leftarrow \widehat{\boldsymbol{\beta}}^{(b)} - \eta \frac{\partial \mathcal{L}(\boldsymbol{\beta})}{\partial \boldsymbol{\beta}}
		=\widehat{\boldsymbol{\beta}}^{(b)} - \frac{\eta}{N_b} (\boldsymbol{S}^{(b)}\widehat{\boldsymbol{\beta}}^{(b)} - \boldsymbol{U}^{(b)}),
	\end{equation}
	where $ \eta $ is the learning rate in the gradient descent;
	\item 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{\boldsymbol{\beta}}^{(b)}$  in step 1, for $r=1,\dots,p$, where
	$\mathcal{S}(x, \lambda)=\text{sgn}(x)(|x|-\lambda)_+$.

\end{itemize}
These two steps are carried out iteratively till convergence. In the implementation, the stopping criterion is set as $\|{\partial \mathcal{L}(\boldsymbol{\beta})}/{\partial\boldsymbol{\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
\begin{equation*}
	\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)}.
\end{equation*}
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
\begin{equation}\label{eq:sigma_sq}
	(\widehat{\sigma}^2_{\epsilon})^{(b)} \equiv \frac{N_{b-1}}{N_b}(\widehat{\sigma}_{\epsilon}^2)^{(b-1)} + \frac{n_b}{N_b} (\boldsymbol{y}^{(b)} - \boldsymbol{X}^{(b)}\widehat{\boldsymbol{\beta}}^{(b)})^\top (\boldsymbol{y}^{(b)} - \boldsymbol{X}^{(b)}\widehat{\boldsymbol{\beta}}^{(b)}),
\end{equation}
which will be used in constructing the confidence intervals.

\subsection{Online low-dimensional projection}\label{sec: online projection}

Next, we obtain an online estimator for the low-dimensional projection.  Let
\begin{equation}\label{projection_step_j}
	\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\},
\end{equation}
where $ N_b$ and $\lambda_b$ are the same as in (\ref{lasso_step_j}).
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 $ \boldsymbol{S}^{(b)} $ directly with $\bm{R}^{(b)}=\boldsymbol{S}^{(b)}_{-r,-r}$ and $\bm{T}^{(b)}=\boldsymbol{S}^{(b)}_{-r,r} $, where $\boldsymbol{S}^{(b)}_{-r,-r}$ is a sub-matrix of $\boldsymbol{S}^{(b)}$ excluding the $r$-th row and the $r$-th column, and $\boldsymbol{S}^{(b)}_{-r,r} $ is a sub-matrix of $\boldsymbol{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~\ref{sec: online debias}

\subsection{Online debiased lasso estimator}\label{sec: online debias}
When data batch  $\mathcal{D}_b$ arrives,  the ODL
estimator  for $\beta_{0, r}$, $r=1,\ldots, p$, is defined as
\begin{equation}\label{online_debiased_algorithm}
	\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}
\end{equation}
Although  the historical data are  still involved in \eqref{online_debiased_algorithm}, 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 \eqref{online_debiased_algorithm}, we obtain the ODL estimator $ {\widehat{\beta}}^{(b)}_{\text{on}, r}. $ Meanwhile, as discussed in Section \ref{sec3}, the estimated standard error is given by $ \widehat{\sigma}^{(b)}_\epsilon\widehat{\tau}^{(b)}_{r} $,
where
\begin{equation}\label{tau_b}
	\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}},
\end{equation}
and $(\widehat{\sigma}^{2}_\epsilon)^{(b)}$  is given in \eqref{eq:sigma_sq} in Subsection \ref{sec: online lasso}.

\subsection{Tuning parameter selection}\label{sec: tuning parameter}
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~\ref{fig:onlineCV}. 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"~\citep{Tashman_recal_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}\|(\boldsymbol{y}^{(b)} - \boldsymbol{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{\boldsymbol{\beta}}^{(b)}(\lambda_b)$ as the starting point for ODL.
\begin{figure}[h]
	\centering
	\includegraphics[width=3.2 in, height=1.6 in]{figures/onlineCV}
	\caption{\label{fig:onlineCV}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{\boldsymbol{\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=\{\boldsymbol{y}^{(b)},\boldsymbol{X}^{(b)}\}$ is the testing set.}
\end{figure}

\subsection{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~\ref{fig:algorithm}. It consists of two main blocks: one is~\textit{online lasso estimation} and the other is \textit{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 \textit{online lasso estimation} block, where the summary statistics $\left\{\boldsymbol{S}^{(b-1)},\boldsymbol{U}^{(b-1)} \right\}$ are updated to $\left\{\boldsymbol{S}^{(b)},\boldsymbol{U}^{(b)} \right\}$. These summary statistics facilitate the updating of the lasso estimator $\widehat{\boldsymbol{\beta}}^{(b-1)}$ to $\widehat{\boldsymbol{\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{\boldsymbol{\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 \textit{low-dimensional projection} block for the calculation of  $\widehat{\boldsymbol{\gamma}}_r^{(b)}(\lambda_b)$. The idea of online updating is the same as in \textit{lasso estimation}, except the relevant summary statistics are the sub-matrices of $\boldsymbol{S}^{(b)}$. The resulting projection $\widehat{\boldsymbol{z}}_r^{(b)}$ output from the \textit{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.

\begin{figure}[h]
	\centering
	\includegraphics[width=\linewidth]{figures/algorithm}
	\caption{\label{fig:algorithm}Flowchart of the online debiasing algorithm. When a new data batch $\mathcal{D}_b$ arrives, it is sent to the  \textit{lasso estimation} block for updating $\widehat{\boldsymbol{\beta}}^{(b-1)}$ to $\widehat{\boldsymbol{\beta}}^{(b)}$. At the same time, it is also viewed as test set for adaptively choosing tuning parameter $\lambda_b$. In the \textit{low-dim projection}  block, we extract sub-matrices from the updated summary statistic $\boldsymbol{S}^{(b)}$ to compute $\widehat{\boldsymbol{\gamma}}_r^{(b)}(\lambda_b)$ and the corresponding low-dimensional projection $\widehat{\boldsymbol{z}}_r^{(b)}$. Outputs $\widehat{\beta}_r^{(b)}(\lambda_b)$ and $\widehat{\boldsymbol{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)})$.}
\end{figure}

\begin{remark}
When $p$ is large,
the online algorithm presented in Figure \ref{fig:algorithm} requires a sufficient large storage capacity, since sample covariance matrix $\boldsymbol{S}^{(b)}$ requires $\mathcal{O}(p^2)$ space complexity. To reduce memory usage, we can apply the eigenvalue decomposition (EVD) of $\boldsymbol{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 $\boldsymbol{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~\citep{onlinePCA2018} 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}^2p) $ 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.
\end{remark}

\section{Theoretical properties}\label{sec3}
To establish the asymptotic properties of the ODL estimator proposed in Section \ref{sec:online_debiased_lasso}, 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 $ \boldsymbol{X} $ and the error terms are imposed to establish the asymptotic results. Specifically, we assume that $ \boldsymbol{X} $ has either i.i.d sub-Gaussian or bounded rows. We first consider the sub-Gaussian case.

\begin{assumption}\label{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.
\end{assumption}

\begin{theorem}\label{thm_online}
	Assume Assumption \ref{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
	\begin{eqnarray}\label{online_additional_assumption}
		s_0s_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),
	\end{eqnarray}
	then, for any $ r = 1, \ldots, p $ and sufficiently large $ N_b $,
	\begin{eqnarray*}
		&&({\widehat{\beta}}^{(b)}_{\text{on},r} - { \beta}_{0,r})/\widehat{\tau}^{(b)}_{r}= W_r + \Delta_r,\\
		&&W_r=\frac{\bm{\widehat{z}}_{r}^\top \bm \epsilon}{\lVert\bm{\widehat z}_{r}\lVert_2},\ |\Delta_r| = o_{\mathbb{P}}(1),
	\end{eqnarray*}
	where $ \widehat{\tau}^{(b)}_{r} $ is defined in \eqref{tau_b}.
\end{theorem}

\begin{remark}\label{remark 2}
	\textit{
		Similar to the  offline debiased lasso estimator \citep{Zhang_delasso_2014, vandegeer2014},
		Theorem \ref{thm_online} 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  \citet{Zhang_delasso_2014}
		and \citet{vandegeer2014} heavily use the KKT condition. Therefore, different arguments are needed to
		establish Theorem \ref{thm_online}.
	}
\end{remark}

\begin{remark}
	\textit{Theorem \ref{thm_online} is established for the proposed online debiased lasso estimators based on the algorithm described in Section \ref{sec:online_debiased_lasso}. Indeed, the proof of Theorem  \ref{thm_online}  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
		~\citet{truncated_SGD_Langford_2009} and \citet{online_sparse_Fan_2018} will have similar asymptotic distributional properties.
	}
\end{remark}

\begin{remark}
	\textit{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 \citet{vandegeer2014}. }
\end{remark}

The requirement on the minimum batch size in Theorem \ref{thm_online} 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$.

\begin{theorem}\label{thm_online_2}
	Assume Assumption  \ref{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
	\begin{eqnarray*}\label{online_additional_assumption_2}
		s_0s_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),
	\end{eqnarray*}
	then, for any $ r = 1, \ldots, p $ and sufficiently large $ N_b $,
	\begin{eqnarray*}
		&&({\widehat{\beta}}^{(b)}_{\text{on},r} - { \beta}_{0,r})/\widehat{\tau}^{(b)}_{r}= W_r + \Delta_r,\\
		&&W_r=\frac{\bm{\widehat{z}}_{r}^\top \bm \epsilon}{\lVert\bm{\widehat z}_{r}\lVert_2},\ |\Delta_r| = o_{\mathbb{P}}(1),
	\end{eqnarray*}
	where $ \widehat{\tau}^{(b)}_{r} $ is defined in \eqref{tau_b}.
\end{theorem}

\begin{remark}
	\textit{The requirement of the first batch size $ n_1 \geq c_1s_r\log p $  in Theorems \ref{thm_online} and  \ref{thm_online_2}
		is needed to establish the consistency of the lasso-typed estimator \citep{buhlmann2011statistics}; otherwise, the error bound of
		$ \bm{\widehat{\gamma}}^{(1)}_{r}$ defined in \eqref{projection_step_j} 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}  \} $.
	}
\end{remark}

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}|$.

\begin{assumption}\label{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) $.
\end{assumption}

\begin{theorem}\label{thm_online_3}
	Assume Assumption \ref{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
	\begin{eqnarray*}\label{online_additional_assumption_3}
		s_0s_r \frac{\log (p)}{\sqrt{N_b}} = o(1), \ s_r^2{\frac{\log (p)}{N_b}\log{{N_b}}} = o(1),
	\end{eqnarray*}
	then, for any $ r = 1, \ldots, p $ and sufficiently large $ N_b $,
	\begin{eqnarray*}
		&&({\widehat{\beta}}^{(b)}_{\text{on},r} - { \beta}_{0,r})/\widehat{\tau}^{(b)}_{r}= W_r + \Delta_r,\\
		&&W_r=\frac{\bm{\widehat{z}}_{r}^\top \bm \epsilon}{\lVert\bm{\widehat z}_{r}\lVert_2},\ |\Delta_r| = o_{\mathbb{P}}(1),
	\end{eqnarray*}
	where $ \widehat{\tau}^{(b)}_{r} $ is defined in \eqref{tau_b}.
\end{theorem}

Theorems \ref{thm_online_2} and \ref{thm_online_3} are established without specific assumptions on data batch sizes except for the first batch. Comparing to Theorem~\ref{thm_online_2},
Theorem \ref{thm_online_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} (\boldsymbol{y}^{(j)} - \boldsymbol{X}^{(j)}\widehat{\boldsymbol{\beta}}^{(j)})^\top (\boldsymbol{y}^{(j)} - \boldsymbol{X}^{(j)}\widehat{\boldsymbol{\beta}}^{(j)}), $ we can see that
$ \widehat{\sigma}^{(b)}_\epsilon $ is consistent for $ \sigma_\epsilon $  in view of  the consistency of $ \widehat{\boldsymbol{\beta}}^{(b)} $ in Lemma \ref{lemma_1_2}.

The proofs of Theorems  \ref{thm_online}--\ref{thm_online_3} are  included in the  appendix.
According to Theorems \ref{thm_online} - \ref{thm_online_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
\begin{equation*}
	{\widehat{\beta}}^{(b)}_{\text{on},r} \pm  \Phi^{-1}(1 - \frac{\alpha}{2}) (\widehat{\sigma}^{(b)}_\epsilon \widehat{\tau}^{(b)}_{r}),
\end{equation*}
where $ \widehat{\sigma}^{(b)}_\epsilon $ is defined in \eqref{eq:sigma_sq},   $\Phi(\cdot)$ is the cumulative distribution function of the standard normal distribution and $\Phi^{-1}$ is its inverse function.

\section{Simulation experiments}\label{sec:sim}
\subsection{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 = \boldsymbol{\beta}_0^\top\boldsymbol{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})$, $\boldsymbol{x}_i^{(j)}\sim \mathcal{N}(\bm{0},\boldsymbol{\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 $\boldsymbol{\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  $\boldsymbol{\Sigma}$ are considered:
(a) $\boldsymbol{\Sigma}=\boldsymbol{I}_p$; (b) $\boldsymbol{\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 $\boldsymbol{\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~\texttt{lm} as a benchmark for comparison.
The results are reported in Tables~\ref{tab:p_400_ind}-\ref{tab:p_1000_AR1}.

\subsection{Bias and coverage probability}
It can be seen from Tables~\ref{tab:p_400_ind}-\ref{tab:p_1000_AR1} 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~\ref{tab:p_400_ind} and~\ref{tab:p_400_AR1}. 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~\ref{fig:QQ_400}. 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~\ref{tab:p_1000_ind}-\ref{tab:p_1000_AR1}. We believe this is reasonable, as each column in the design matrix, denoted by $\boldsymbol{x}_r\in\mathbb{R}^{N_b\times 1}$, is of the same marginal distribution, and thus the estimated standard errors computed according to \eqref{tau_b}  in the simulations are identical up to a certain decimal for every component in $\boldsymbol{\beta}$,  regardless of signal strength.

\begin{table}
	\caption{\label{tab:p_400_ind} $N_b=420$, $b=12$, $p=400$, $s_0=6$, $\boldsymbol{\Sigma}= \boldsymbol{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~\ref{sec: tuning parameter}. Simulation results are summarized over 200 replications. In the table, we report the $\lambda$ selected with highest frequency among 200 replications.}
	\centering
	\begin{tabular}{l l| c | cccccc}

		&$\beta_{0,r}$ &OLS  &\multicolumn{6}{c}{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 \\

		\hline
		\hline
		\multirow{3}{*}{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\\
		\hline
		\multirow{3}{*}{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\\
		\hline
		\multirow{3}{*}{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\\
		\hline
		\multirow{3}{*}{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\\
		\hline
		\multirow{3}{*}{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\\
		\hline
	\end{tabular}
\end{table}

\begin{table}
	\caption{\label{tab:p_400_AR1} $N_b=420$, $b=12$, $p=400$, $s_0=6$, $\boldsymbol{\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~\ref{sec: tuning parameter}. Simulation results are summarized over 200 replications. In the table, we report the $\lambda$ selected with highest frequency among 200 replications.}
	\centering
	\begin{tabular}{l l| c| ccc ccc}
		&$\beta_{0,r}$ &OLS  &\multicolumn{6}{c}{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\\
		\hline
		\hline
		\multirow{3}{*}{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\\
		\hline
		\multirow{3}{*}{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\\
		\hline
		\multirow{3}{*}{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\\
		\hline
		\multirow{3}{*}{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\\
		\hline
		\multirow{3}{*}{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\\
		\hline
	\end{tabular}
\end{table}

\begin{figure}
	\centering
	\includegraphics[width=\linewidth]{figures/p_400_ind}
	\caption{\label{fig:QQ_400}QQ plots of standardized $\widehat{\beta}_{\text{on},r}^{(b)}$ with total sample size $N_b=420$, $p=400$ and $\boldsymbol{\Sigma}= \boldsymbol{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 $\boldsymbol{\beta}_0$, i.e. $\beta_{0,r}=0, 0.01, 1$. }
\end{figure}

\begin{table}
	\caption{\label{tab:p_1000_ind} $N_b=1200$, $b=12$, $p=1000$, $s_0=20$, $\boldsymbol{\Sigma}= \boldsymbol{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~\ref{sec: tuning parameter}. Simulation results are summarized over 200 replications. In the table, we report the $\lambda$ selected with highest frequency among 200 replications.}
	\centering
	\begin{tabular}{l l|  c | cccccc}
		&$\beta_{0,r}$ &OLS  &\multicolumn{6}{c}{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\\

		\hline
		\hline
		\multirow{3}{*}{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\\
		\hline
		\multirow{3}{*}{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\\
		\hline
		\multirow{3}{*}{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\\
		\hline
		\multirow{3}{*}{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 \\
		\hline
		\multirow{3}{*}{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\\
		\hline
	\end{tabular}
\end{table}

\begin{table}
	\caption{\label{tab:p_1000_AR1} $N_b=1200$, $b=12$, $p=1000$, $s_0=20$, $\boldsymbol{\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~\ref{sec: tuning parameter}. Simulation results are summarized over 200 replications. In the table, we report the $\lambda$ selected with highest frequency among 200 replications.}
	\centering
	\begin{tabular}{l l|c | cccccc}
		&$\beta_{0,r}$ &OLS  &\multicolumn{6}{c}{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 \\

		\hline
		\hline
		\multirow{3}{*}{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 \\
		\hline
		\multirow{3}{*}{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\\
		\hline
		\multirow{3}{*}{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\\
		\hline
		\multirow{3}{*}{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\\
		\hline
		\multirow{3}{*}{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\\
		\hline

	\end{tabular}
\end{table}

\section{Applications}
\subsection{Analysis of  Beijing PM 2.5 data}\label{sec:data}
We apply the proposed ODL to analyze the Beijing PM2.5 Data by \cite{liang2015assessing}, 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~\ref{fig:real_wind_season}. 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~\cite{liang2015assessing}. 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.

\begin{figure}
	\centering
	\includegraphics[width=0.31\textwidth]{figures/p1}
	\includegraphics[width=0.31\textwidth]{figures/p2}
	\includegraphics[width=0.31\textwidth]{figures/p7}
	\caption{\label{fig:real_wind_season} 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.}
\end{figure}

Next, we focus on another two variables: pressure and dew point. The results are shown in Figure \ref{fig:real_PD_season}. 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 \cite{liang2015assessing}. 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 \ref{fig:real_wind_season}. 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 \ref{fig:trace}. 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.

\begin{figure}
	\centering
	\includegraphics[width=0.31\textwidth]{figures/p4}
	\includegraphics[width=0.31\textwidth]{figures/p5}
	\includegraphics[width=0.31\textwidth]{figures/p8}
	\caption{\label{fig:real_PD_season} 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.}
\end{figure}

\begin{figure}
	\centering
	\includegraphics[width=0.85\linewidth]{figures/full_trace}
	\caption{\label{fig:trace} 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.}
\end{figure}

\subsection{Hang Seng Index fund data}
{\color{black} 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 \cite{lan2016testing},
 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 \ref{fig:HSI_stock}.
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 \ref{fig:HSI_stock},  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.

\begin{figure}
	\centering
	\includegraphics[width=0.48\textwidth]{figures/stock_coeff}
	\includegraphics[width=0.48\textwidth]{figures/stock_std}
	\caption{\label{fig:HSI_stock} 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.}
\end{figure}

}

\section{Discussion}\label{sec: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.

\section*{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.

\bibliographystyle{authordate1}