\section{Executive Summary}

The Laser Interferometer Space Antenna (LISA) will open three decades
of gravitational wave (GW) spectrum between 0.1 and 100\,mHz, the
mHz band \citep{LISAProposal2017}. This band is expected to be the
richest part of the GW spectrum, in types of sources, numbers of sources,
signal-to-noise ratios and discovery potential. When LISA opens the
low-frequency window of the gravitational wave spectrum, around 2034,
the surge of gravitational-wave astronomy will strongly compel a subsequent
mission to further explore the frequency bands of the GW spectrum
that can only be accessed from space. The 2020's is the time to start
developing technology and studying mission concepts for a large-scale
mission to be launched in the 2040's. The mission concept would then
be proposed to Astro2030.

Only space-based missions can access the GW spectrum between $10^{-8}$
and 1\,Hz because of the Earth's seismic noise. This white paper
surveys the science in this band and mission concepts that could accomplish
that science. The proposed small scale activity is a technology development
program that would support a range of concepts and a mission concept
study to choose a specific mission concept for Astro2030. In this
white paper, we will refer to a generic GW mission beyond LISA as
\emph{bLISA.}

\section{Advancing mHz Gravitational Wave Astronomy Beyond LISA}

Gravitational Wave astronomy is a new and promising field. LIGO \citep{aLIGO15}
has shown that gravitational-wave observatories can make routine observations
of sources that are invisible to electromagnetic (EM) observations,
and of sources for which complementary EM and GW information is extraordinarily
powerful. In the first ever observation of GWs, GW150914 \citep{GW150914},
LIGO demonstrated the discovery potential of GW observations, by detecting
merging stellar black holes with masses substantially higher than
expected. The subsequent observation of a merging neutron star binary,
GW170817 \citep{GW170817}, that was also widely observed across the
EM spectrum simultaneously showed the power of coordinated observations
and the power of multiple detectors (i.e., LIGO Hanford, LIGO Livingston,
and Virgo). The extraordinary campaign of EM observations that followed
GW170817 provided a tremendous impetus to multi-messenger and time
domain astronomy.

\begin{figure}
\noindent \begin{centering}
\includegraphics[width=0.9\textwidth]{GWplotter.png}
\par\end{centering}
\caption{GW Spectrum, showing frequency bands where PTAs, LISA and ground-based
GW detectors operate (generated by http://gwplotter.com/, see \citep{GW_Plotter}
for details).}

\noindent \centering{}\label{fig:GWspectrum}
\end{figure}

The GW spectrum in Figure~\ref{fig:GWspectrum} shows the detection
strategies for GW astronomy. Roughly speaking, GW sources are inspiraling
binaries (or triples) of compact objects whose GW frequency chirps
up to a final merger frequency, followed by ring-down of the final
object. The gravitational wave frequency roughly scales inversely
with system mass. Pulsar timing arrays (PTAs) detect the largest SMBHs
over a decade of frequency in the nHz region, and ground-based interferometers
detect stellar-mass systems between 40 and 1000\,Hz. LISA spans 0.1
to 100\,mHz. All three strategies are amplitude – not power – detectors,
meaning signal falls off as $1/\text{R}$, rather than $1/\text{R}^{2}$.

Based on the success of LISA Pathfinder \citep{LPF18a,LPF18,LPF18b,LPF19a,LPF19b,LPF19c}
ESA is leading the LISA mission to open up the mHz GW band, with NASA
as a junior partner. This band is expected to have many more detectable
types of sources than any other in the GW spectrum: $10^{4}-10^{7}{\rm M}_\odot$
massive black hole binaries (MBHBs), intermediate mass black hole
binaries (IMBHBs), extreme-mass-ratio inspirals (EMRIs) with mass
ratios of $10^{4}$ to $10^{5}$ with a system mass $<10^{7}{\rm M}_\odot$,
intermediate-mass-ratio inspirals (IMRIs) with smaller system masses
of $10^{3}-10^{4}{\rm M}_\odot$, close compact binaries in the Milky Way,
and the heavy stellar binaries (10's of ${\rm M}_\odot$) seen by LIGO and
Virgo. Larger sources will be detectable back into the re-ionization
era. The number of sources is expected to be 10's of thousands. Signal-to-noise
ratios (SNRs) can range into the thousands for the strongest signals.
Further, the detectable mass of sources at cosmological distances
is augmented by their redshift, a substantial boost at, say, $z\sim10$.
LISA's discovery potential includes cosmic strings, cosmological phase
transitions, unexpected bursts and a stronger cosmological GW background
than standard inflationary models predict. LISA's impact on GW astronomy
will parallel LIGO's impact and a space-based GW detector to follow
LISA will be critical to NASA's portfolio.

Below 10\,nHz, pulsar timing arrays (PTAs) observe through long-term,
precisely timed observations of rapidly rotating radio pulsars. It
is widely anticipated that they will measure signals from the largest
MBH binaries $(\sim10^{9}{\rm M}_\odot)$ at low redshifts in the 2020's.
PTA sensitivity improves with the discovery and long-term observation
of highly stable pulsars \citep{IPTA16}.

At the high frequency end of the spectrum (40-1000 Hz), ground-based
GW detectors like LIGO and Virgo are currently alternating between
observations that are producing ever more sources, and upgrades leading
to increased sensitivity and duty cycle. Increased sensitivity improves
SNR and enlarges the observable volume as the cube because these are
amplitude detectors. Improved seismic isolation and suspensions are
extending the useful observing band down towards 10\,Hz. During the
current observing run, slated to end in 2020 after a year of operation,
LIGO and Virgo are releasing detection alerts about once per week.
Upgrades are scheduled to continue for several years. Sky localization
will continue to improve as new detectors, KAGRA \citep{Kagra} in
Japan and LIGO-India, come online over the same time scale. So-called
3G detectors \citep{EinsteinTelescope,3G}, with longer armlengths,
improved suspensions and other technology improvements, are in the
early planning stages, with implementation notionally in the 2030s.
These detectors would further improve sensitivity, extend their operating
band down towards 1\,Hz and capitalize on the benefits of a global
network for the best sky localization.

To advance beyond LISA, we survey the science goals and then consider
two illustrative examples of mission concepts. So far in advance of
LISA, it is unwise to focus on a single design. It is more prudent
to examine a range of designs, evaluate the technology needs, some
of which are shared, and down-select to a concept when the technology
challenges are better understood and science priorities are more refined.
The schedule advanced here is for technology development and a mission
concept study in the 2020s leading to a concept selection in time
to propose to Astro2030. Ideally, bLISA would start implementation
when LISA is in its extended mission.

There have been previous efforts to examine possible missions succeeding
LISA. One of the earliest was the Big Bang Observer \citep{BBO2003},
a visionary concept to detect the cosmological GW background. Crowder
and Cornish \citep{CrowderCornish05} compared LISA, BBO, ALIA \citep{ALIABender13}
(to be described below) and stereo versions of LISA and ALIA.

In 2012, a NASA study team produced a wide-ranging study of alternative
designs to LISA \citep{NASASGO12}. While this study did not address
science beyond LISA science, it did comprehensively survey the architectural
choices for space-based GW detectors. The low-frequency Folkner design
described below was derived from this study. The GADFLI \citep{GADFLI}
and GEOGRAWI (now published under the acronym gLISA) \citep{gLISA}
concepts, mentioned below, also originated from this study.

\section{Key Science Goals}

Similar to the EM-spectrum, GW sources are present throughout the
universe and emit over many decades in frequency space. GW emission
is directly linked to the mass of the emitting system; heavier binary
systems merge at lower frequencies while lighter masses pass through
the low frequency band, often with detectable amplitudes, and merge
then at higher frequencies. Closing frequency gaps in the observed
spectrum has two distinct scientific ramifications: (1) mergers can
be measured across (nearly) all mass ranges and out to very large
redshifts; and (2) the inspiral phase can be tracked over many months
and often over many years prior to the merger. Both improvements provide
unique opportunities to test GR. At least as important are improvements
in the 3D localization of these sources to enable coordinated EM observations.

Improving angular resolution is an age-old goal for many areas of
astronomy, and this may be especially important for gravitational
wave astronomy, where the GW signal unequivocally identifies the source
while EM and particle observatories identify the EM/particle signatures
of these sources. LISA localizations may typically be on the order
of a few square degrees (varying significantly for different sources).
With incident wavelengths measured in millions of km, the Rayleigh
criterion $\delta\theta\leq\lambda/D$ indicates the challenge of
directly resolving gravitational wave sources. A large fraction of
GW science, however, relies primarily on astrometric location of point
sources, enhanced beyond the resolving power by the SNR ($\rho$).
So increases in sensitivity are less important for a census of the
black hole population but are crucial to improve the angular resolution
and enable multi-messenger observations.

A key detection milestone will be to achieve angular (and distance)
location with sufficient precision to localize extragalactic events
such as EMRI's and merged massive or super-massive BHB systems to
a single galaxy (or cluster). Approximate localization of these BHB
systems to within arcminutes in advance of merger would enable deep
coincident observations to identify an EM counterpart and thus locate
the galaxy. These multi-messenger observations would then also enable
studies of accretion flows and physics along with details of galaxy
properties, all with precise knowledge of the black hole masses, spins
and recent merger history. Some MBHB galaxy identifications may be
possible with LISA, but future GW missions, in concert with advanced
EM facilities should be able to dramatically expand the rich multi-messenger
data set to enable a robust understanding of the roles black holes
and galaxies play in shaping each other through accretion and mergers.

\subsection{The mHz to Hz frequency band}

The mHz to Hz frequency band is the band in which intermediate mass
black hole binaries, even from lower-mass MBH seeds from Pop-III stars,
would merge. A future observatory with significantly improved sensitivity
compared to LISA provides clear statistics on these still uncertain
objects. Beyond proving their existence and identifying their role
between stellar mass BHs and massive BHs, precise measurements of
the phasing of these mergers would probe for extensions of GR, testing
the presence of additional physical fields with (effective) mass (e.g.
ultralight fields/dark matter candidates).

Another unique source for this frequency range are (typically extra-galactic)
double white dwarf (DWD) systems which could be observed at periods
around their point of initial contact, allowing for example the observation
of a gravitational-wave signal accompanying a type Ia supernova. Based
on the known rate of type Ia supernovae, this also requires significant
improvements in the high frequency sensitivity of LISA. Such a DWD-merger
signal would provide unprecedented early warning for a supernovae.

Observations in this band also provide improved and advanced localization
of the stellar-mass NS and BH chirping binaries with mergers observed
by 3G ground-based GW observatories. Mergers would be precisely forecast
enabling unique options for multi-messenger observations. If intermediate
mass BHs exist in the universe, joint observations of these systems
on the ground and in space will allow us to localize these sources
well before they merge and to break parameter estimation degeneracies.
At sufficiently high sensitivities it could be possible to independently
resolve all such systems to high redshift, providing a complete census
of these systems, while clearly exposing any primordial background
like that predicted by some models of inflation.

LISA will be able to detect extreme mass-ratio inspirals (EMRIs) composed
of a $10\,{\rm M}_\odot$ black hole falling into a $3\times10^{5}{\rm M}_\odot$
BH to as far as $z=3$, but this horizon decreases rapidly for different
BH mass, or for smaller secondary masses. While the rates of these
remain quite uncertain, we might expect LISA to detect a few per year
based on current models, but increased sensitivities just above a
mHz would enhance the rate roughly as $\rho^{3}$. Plausible sensitivity
improvements would provide a census of the EMRI population to moderate
redshift, providing insights into the dynamics in stellar clusters.
More sensitive observations of EMRIs would enable high precision spacetime
mapping to verify the Kerr nature of astrophysical black holes. For
the strongest EMRI signals we might even gain access to the overtones
in the ring-down radiation.

\subsection{The sub-$\mu\text{Hz}$ to mHz band}

\begin{wrapfigure}{l}{0.6\textwidth}%
\noindent \centering{} \includegraphics[width=0.6\textwidth]{charstrain}
\caption{Two potential future mission designs compared to the current LISA
design. The low frequency, or Folkner design, assumes that we can
extend LISA's acceleration noise into the low frequency band while
the high frequency design (ALIA) assumes a frequency independent order
of magnitude improvement in the current acceleration noise. The graph
also shows the traces of a few equal-mass black hole mergers as examples
of the science return of each of these missions. The grey shaded area
is a potential representation of the still unknown stochastic gravitational
wave emission from the millions of galactic binaries.}
\label{fig:sensitivity} \end{wrapfigure}%

A similar gain in sensitivity and therefore angular resolution in
the sub-mHz frequency band is frustrated by the stochastic background
of the myriad of GW signals from galactic binary systems. The exact
shape of this spectrum is subject to debate; current and future survey
missions as well as LISA will provide some answers over the next 20
years. However, it is expected that the sensitivity of LISA might
be limited by this background below about one mHz. The grey area in
Figure \ref{fig:sensitivity} shows one of many possibilities for
this background. \emph{One person's noise is another person's signal.}
Any future space mission probing below the LISA band will be able
to measure the spectral and spatial distribution of these binaries
over more than two decades in frequency space and will be able to
identify and isolate probably hundreds of thousands of the stronger
binary systems in our galaxy.

The different models of the galactic binary background still allow
detection of gravitational waves from merging super-massive black
holes in the $10^{7-9}{\rm M}_\odot$ mass range out to large redshifts. These
signals will not only map out the Kerr nature of super-massive black
holes with high SNR but also shed light on many of the most fundamental
questions in cosmology and galaxy evolution.

One example is the still unsolved puzzle of the rapid emergence of
the high-z QSOs with big implications not only on the cosmic evolution
of super-massive black holes but also on their impact on early galaxy
formation. Do they regulate the entropy of the intergalactic medium,
hence the fuel of galaxy formation? LISA will probe the $10^{4-6}{\rm M}_\odot$
seeds of the early QSOs but will not tell us how and how fast they
reached $10^{8-9}{\rm M}_\odot$ (e.g., by Super-Eddington accretion). At
high z it is expected that merger rates of galaxies and subsequent
mergers between their central MBHs is very high. The observation of
the generated GWs is a good probe of the growth function and is probably
the only tool for a reliable census of MBHs during this epoch. In
the mid range of this low frequency detector, we could possibly see
thousands of inspiraling MBHBs with periodicities of days accessible
with future (post LSST) time domain surveys. GWs would provide the
luminosity distance and the EM-observatories the redshift for these
standard sirens. These multi-messenger observations would also provide
powerful insight onto accretion on binaries, the physics of disk-binary
interaction, and many other aspects which need the masses and spins
of the central engines behind these sources. Sensitivity at the high
frequency end of this detector enables earlier detection and localization
of LISA's MBHBs which will greatly expand opportunities for multi-messenger
observations before and during merger.

\section{Example Mission Concepts}

Looking to the future, we need to explore the science opportunities
beyond LISA, and also the technological hurdles necessary to realize
a beyond-LISA mission. Figure\,\ref{fig:sensitivity} compares the
sensitivity curves of two example mission designs to LISA. These designs
were chosen to illustrate the range of mission concepts, and to show
what is possible if a healthy technology development program for future
gravitational wave missions is established parallel to the LISA project
itself.

\begin{wrapfigure}{r}{0.4\textwidth}%
\noindent \centering{} \includegraphics[width=0.4\textwidth]{Folkner_orbits}
\caption{The Folkner low-frequency mission concept placed three spacecraft
in this triangular configuration around the Sun. }
\label{fig:orbits} \end{wrapfigure}%

One mission design explores the frequency band below LISA and is based
on orbits in which the three spacecraft form an equilateral triangle
with the Sun in its center. The arms are about 100 times longer than
the LISA arms. This design keeps one spacecraft in the vicinity of
Earth to simplify communication with the constellation. Communication
to the other spacecraft would probably have to be carried by the laser
links. These orbits were originally suggested by William Folkner during
the GW Mission Concept Study \citep{NASASGO12} prior to the LISA
Pathfinder success. He predicted some sensitivity to large massive
black hole merger systems in a mission without inertial sensors where
each spacecraft would act as a test masses. We added a drag free system
back into the spacecraft and assumed that it maintains a frequency
independent acceleration noise identical to the LISA requirement.

The second mission design, the Advanced Laser Interferometer Antenna
(ALIA) \citep{ALIABender13}, assumes a ten times improved acceleration
noise in the high frequency part of the LISA spectrum. This example
mission is five times shorter than LISA but uses significantly more
laser power and larger telescopes to increase the interferometric
sensitivity. The orbits for ALIA are still heliocentric but other
mission designs with shorter arms explored also geocentric orbits
(GADFLI and GEOGRAWI, a.k.a. gLISA).

\subsection{Low-Frequency Mission Example}

The orbits proposed in the Folkner mission would separate the spacecraft
by around 260\,Gm (see Fig.\,\ref{fig:orbits}). One of them has
to be kept at a reasonable distance to Earth to allow for spacecraft
to ground communication which means that lasercom between the spacecraft
has to be part of the mission design.

Folkner also provided an initial study of the orbital dynamics for
the first five years of the mission. He found that differential velocities
will stay below 3m/s, which affects the Doppler shifts and therefore
the frequencies of the beat signals between the laser beams. Furthermore,
changes in the opening angles between the spacecraft stay below $0.02^{\circ}$
and would likely not require in-flight adjustments of the opening
angle between the two outgoing beams on each spacecraft.

For the long arm interferometer, we used parameters which are very
similar to LISA: 30~cm telescopes and a 3\,W laser system. The received
power at the far end would be approximately 150\,fW and the shot
noise limited displacement sensitivity would be a factor 200 below
LISA leading to the same strain sensitivity when expressed as a linear
spectral density; the sensitivity curves in Figure\,\ref{fig:sensitivity}
are scaled by $\sqrt{f}$ to take into account the longer observation
time for low-f signals. As shown in Figure\,\ref{fig:sensitivity},
the galactic binary background will likely limit the sensitivity well
before shot noise.

\subsection{High-Frequency Mission Example}

We use the ALIA mission design as an example for a high frequency
beyond LISA mission. ALIA's orbits are LISA-like heliocentric orbits
trailing or leading Earth by a few degrees. The distance between the
spacecraft is 500,000\,km; five times shorter than LISA. The resulting
loss in low frequency sensitivity is (over-) compensated by an assumed
ten-fold improvement in acceleration noise. These two parameters can
be fine-tuned once the galactic binary background radiation has been
characterized by LISA. The main improvement comes from a 40,000-fold
increase in the received laser power which ALIA achieves by increasing
the telescope diameter to 1\,m and the laser power to 30\,W. This
leads to a factor 200 improvement in phase/displacement sensitivity
and, due to the factor 5 shorter arms, to a 40-fold improvement in
strain sensitivity. The dynamics of the orbits, relative spacecraft
velocities and angular changes, scale with the arm length (assuming
all other parameters stay the same) and will be reduced by a factor
five compared to LISA.

The spacecraft in the geostationary mission designs were separated
by 73,000\,km. The shorter arm length shifts the sensitivity curve
further to the right. Using ALIA's laser power and telescope diameter
would give the same strain sensitivity at higher frequencies; the
increase in received laser power is compensated by the decrease in
length. However, it requires further improvements in the phase sensing
system. The advantage of geostationary orbits is the reduced launch
costs and simple ground to space communication links. On the other
hand, geostationary constellations likely require station-keeping
to manage the relative velocities (Doppler shifts) and changes in
the opening angles.

\section{Technology}

Space-based GW detectors, based on laser interferometry, generally
have two noise regimes that determine the performance. Residual accelerations
on the inertial reference masses from unwanted disturbances limit
the low frequency performance. Displacement measurement noise from
the interferometry limit the high frequency performance. Together
with the frequency response of a chosen armlength to the GW wavelength,
these two noise types give the bucket-shaped noise curves that characterize
laser-interferometer-based detectors. Any technology improvements
naturally have to address these noise types. Changing a detector's
operating band also requires adjustments to the treatment of these
noise types.

\subsection{Acceleration Noise}

LPF improved on existing drag free systems by several orders of magnitude
and represents the state of the art in force free motion. The LISA
Pathfinder team has done a tremendous job of understanding and characterizing
the limiting noise sources of their gravitational reference sensor
(GRS) which is now baselined for LISA. However, the initial concept
was born in the mid '90s based on a design sensitivity which was defined
prior to understanding the true limitations of such a system. Now
we have a much better understanding of these limitations and, with
that knowledge, should be able to take a fresh look at the design
and find ways to improve it.

\begin{wrapfigure}{r}{0.6\textwidth}%
\noindent \centering{} \includegraphics[width=0.6\textwidth]{figure_SingleTMAccNoise_Alloc_and_CBE_breakdown}
\caption{The CBE of LISA's acceleration noise is based on LPF data and offers
opportunities for future improvements \citep{LISAPerfModel}.}
\label{fig:acceleration_noise} \end{wrapfigure}%

The LPF team has analyzed hundreds of different potential noise sources.
The leading contributors together with the LISA requirement and current
best estimate are shown in Figure \ref{fig:acceleration_noise} \citep{LISAPerfModel}.
At high frequencies, the noise in the actuation forces between the
spacecraft and the test mass limit the performance. These forces are
typically described as a stiffness along each translational and rotational
degree of freedom. In a perfect drag free system, the spacecraft is
tracking the motion of the test mass and the forces would be constant.
A real drag free system is limited by sensor noise; in LISA we assume
$\sim\text{nm}/\sqrt{\text{Hz}}$ and a few hundred $\text{nrad}/\sqrt{\text{Hz}}$
noise in the capacitive sensors, and by the response time of the $\mu\text{N}$-thrusters.
The later is responsible for the increased noise at higher frequencies.
A related noise source is gravitational noise or changes in the local
gravitational field due to spacecraft motion. At high frequencies,
this is again coupled to the non-suppressed spacecraft motion mostly
along the non-sensitive axis due to cross coupling between the different
degrees of freedom. In LISA, these degrees are only sensed using the
capacitive sensors. One possible improvement in future missions would
be to add a much more sensitive interferometric sensing system to
monitor all degrees of freedom of testmass to spacecraft motion and
not only the one along the optical axis. The design of such a system,
the required sensitivity along the other degrees of freedom, added
complexity and all ramifications with respect to the payload design
need to be explored.

At lower frequencies, thermo-elastic deformations of the spacecraft
will also contribute to the overall noise budget. This noise falls
off with distance cubed. Potential needed changes include improved
thermal stability which has been identified as one of the driving
forces behind many of the low frequency noise terms. In addition to
passive shielding, LISA will also actively stabilize the temperature
of the outer thermal shields which encapsulate the sensitive parts
of the payload. The requirements on this shielding and the active
temperature stabilization system has to be extended towards much lower
frequencies; again a potentially solvable problem with significant
design implications. A parallel approach is to use only low-CTE materials
(with low water content) in the vicinity of the test mass to minimize
thermo-elastic deformations. Low-CTE materials in this context might
not be restricted to ULE or Zerodur but might also include support
structures build from negative and positive CTE materials of similar
density.

The limiting force in LPF in the most sensitive frequency range is
Brownian noise caused by residual gas molecules bouncing off the test
mass. This frequency independent Brownian motion scales with the area
to mass ratio of the test mass, the residual gas pressure which in
LPF was estimated to be $\sim2\,\mu\text{Pa}$ at the end of life,
and on the gap size between test mass and housing. LISA carries a
requirement of $1\,\mu\text{Pa}$ to provide additional margin, which
is still well above the ultra-high vacuum pressure that can be achieved
in modern laboratories and many orders of magnitude away from the
residual pressure in orbit. The gap size was optimized to reduce this
gas pressure noise while maintaining enough sensitivity in the capacitive
sensors and control authority in the electrostatic actuators. 1/f-noise
in the electrostatic actuators was one of the dominant frequency dependent
noise sources limiting LPF at low frequencies. Alternative mission
concepts might use a single test mass which is gravitationally balanced
within the center of the spacecraft. If feasible, such a design could
significantly reduce the needs for test mass actuation and allow to
increase the gap size by maybe even an order of magnitude. Obviously,
such a design would again benefit from an all interferometric sensing
system but other aspects such as redundancy are a major concern.

Another noise source, named stray electrostatic noise in Figure\,\ref{fig:acceleration_noise},
is related to the test mass charge and electrostatic fields. LPF used
a discontinuous UV discharging system which was turned on regularly
to keep the charges below a few million electrons. LISA is planning
to use a continuous UV discharging system which is expected to reduce
the test mass charge significantly and reduce this noise as well.
However, further improvements in the charge sensing technique might
be needed to reduce the residual charge on the test mass again. The
current charge sensing system relies on the electrostatic actuators
which would limit the proposed increase in gap size to overcome other
noise sources.

One example for a noise that scales with the volume of the test mass
is magnetic force noise. Changing magnetic fields couple to the non-vanishing
magnetic susceptibility of the test mass material and to the magnetic
dipole moment from ferromagnetic inclusions. A specific gold-platinum
alloy was chosen because of its vanishing magnetic susceptibility.
A third and, based on LPF experience, likely dominant noise related
to magnetic fields is the interaction and down-conversion of test
mass eddy currents with time dependent magnetic fields, both originating
in the audio band. The current best estimate for LISA is that the
acceleration noise caused by eddy current damping and these magnetic
fields will be at $0.3\,\text{fm}/\text{s}^{2}\sqrt{\text{Hz}}$ above
a few mHz and then increases with $f^{-1}$. Possible mitigation strategies
could include $\mu$-metal and reductions in magnetic fields at audio
frequencies.

All these mitigation steps might allow to reduce the acceleration
noise in the high frequency region by an order of magnitude, as required
by ALIA, and will also increase the frequency band in which the acceleration
noise meets the LISA requirement of $3\,\text{fm}/\text{s}^{2}\sqrt{\text{Hz}}$.
For even lower frequencies, not surprisingly, temperature variations
are among the most crucial disturbances. Temperature fluctuations
and time dependent gradients push and pull on the test mass for example
through differential outgassing and differential radiation pressure
as well as the already discussed gravitational forces. Based on LPF
experience, LISA models assume a steep increase of the temperature
fluctuations towards lower frequencies $(f^{-3.5})$ starting with
$100\,\mu\text{K}/\sqrt{\text{Hz}}$ at $100\,\mu\text{Hz}$. However,
an active temperature sensing and control scheme could significantly
reduce the temperature fluctuations.

To summarize, now is the time to build on the extensive LPF experience
and use the lessons learned to improve acceleration noise beyond its
performance at higher frequencies and expand the performance to even
lower frequencies. The knowledge is fresh, the field is growing, and
starting the technology development now allows to develop a mature
mission concept for the next Decadal survey. But this requires a dedicated
technology development program which builds on the LPF experience
but should also evaluate alternative technologies and methods.

\subsection{Measurement Noise}

A second key technology is required to sense the minute changes in
distance between the two widely separated spacecraft. Laser interferometric
distance measurements are fundamentally limited by the intrinsic phase
noise which is inversely proportional to the amplitude of the metering
laser field. However, reaching this level is often very difficult
especially when the amplitude is large. Ground-based observatories
such as Advanced LIGO use optical cavities to amplify the response
and a Michelson interferometer to suppress the common parts of the
amplitude before they reach and saturate the detector. Both of these
techniques will be very difficult to implement between drag free spacecraft.
LISA is in that sweet spot where the received amplitude is at a level
which allows shot noise limited detection of $\sim10\,\text{MHz}$
beat signals with current technology. LISA still requires clock noise
transfer between the spacecraft and pilot tones to remove phase noise
added by the analog components. A sophisticated timing and ranging
system allows to suppress laser frequency noise using time delay interferometry.
All of these technologies have been demonstrated at the $\sim\text{pm}/\sqrt{\text{Hz}}$-level
in the LISA band sufficient to meet the LISA requirement of $10\,\text{pm}/\sqrt{\text{Hz}}$
equivalent single link displacement noise, twice as high as the allocated
shot noise limit. Figure\,\ref{fig:densing_noise} shows the breakdown
for the sensing noise in the long arm interferometer. It is dominated
by shot noise above 1\,mHz which will also limit the overall LISA
sensitivity at these frequencies. Technical noise sources are allowed
to increase with $f^{-2}$ below 1\,mHz where acceleration noise
is expected to dominate.

\begin{wrapfigure}{r}{0.6\textwidth}%
\noindent \centering{} \includegraphics[width=0.6\textwidth]{Picture2}
\caption{The CBE of the most critical of LISA's long arm interferometer noise
sources \citep{LISAPerfModel}.}
\label{fig:densing_noise} \end{wrapfigure}%

The Folkner mission uses hundred times longer arms which reduces the
received power by four orders of magnitude and increases the shot
noise limit by a factor 100, likely well below the galactic binary
background. If the overall performance of all technical noise sources
in the displacement measurement system continues to raise only with
$f^{-2}$ down to sub-$\mu$Hz, acceleration noise would be a factor
100 above these noise sources. This provides some needed margin to
accommodate expected faster increases in temperature fluctuations
towards the lower end of the $\mu$Hz spectrum which might otherwise
degrade the performance of the phase measurement chain. Note that
increases in the received laser power by either increasing the laser
power itself or the diameter of the telescopes is likely required
to enable for example inter-spacecraft communication or acquiring
lock. There are many other challenges related to this mission including
reaching the orbits, sunshades, and communication, however, the technical
challenges for the principle payload of such a low frequency mission
appear to be mostly related to acceleration noise.

The 200-fold improved shot limit for ALIA present significant challenges
for many subsystems within the interferometric measurement system.
LISA's photo detector noise and phase noise in the analog chain is
already state of the art for the $10-20\,\text{MHz}$ laser beat frequencies.
ALIA's orbit should allow for reduced Doppler shifts resulting in
lower laser beat frequencies if the intrinsic laser noise can be suppressed
below shot noise (1f-RIN-line) at these lower frequencies; the current
laser systems are not shot noise limited below $\sim8\,\text{MHz}$.
However, even if the relative intensity noise can actively be suppressed,
phase noise and timing jitter in the analog parts of the phase measurement
chain still need to be reduced by at least one order of magnitude.
This requires a dedicated effort to develop and study electrical,
electro-optical, and optical components starting from RF cables; the
temperature dependency of the electric susceptibility of the dielectric
inside the cable (often Teflon) is already a concern for phase noise
in LISA, to analog to digital converters and the timing of them.

Beyond the individual components, the phase measurement system in
its entirety needs to be re-evaluated. Are strategies such as comparing
clock noise using the laser links still valid? What alternatives exist?
Are there designs which reduce the laser beat frequencies by a few
orders of magnitude to reduce demands on timing stability. What other
laser sources will be available in time for a new mission? Using a
lower laser wavelength decreases diffraction losses and increases
the phase shift. Furthermore, the amount of scattered light will increase
with the laser power while the sensitivity to scattered light increases
by a factor 200. This requires new strategies to reduce and/or cancel
stray light. Or generally speaking, to take full advantage of the
increased received laser power, LISA's typical $\sim\text{pm}$-requirements
for most technical noise sources in the measurement system turn into
$\sim5\,\text{fm}$-requirements above $\sim20\,\text{mHz}$.

The feasibility of improvements at this magnitude need to be studied
and demonstrated in a dedicated technology development program before
a believable and mature bLISA mission concept can be developed and
presented to the next decadal survey.

\section{Organization, Partnerships, and Current Status}

We imagine that the technology development work proposed here could
be carried out within the structure of the ROSES program, in a manner
analogous to the Beyond Einstein Foundation Science (~2005-2007)
or LISA Preparatory Science (2018) calls. The directed calls could
be formulated by NASA Headquarters APD and PCOS staff with input from
the LISA project team and the concept study. The concept study would
be a panel of experts from the GW and astrophysics community similar
to the many panels created by APD over the last 2 decades. An initial
panel would identify promising concepts and the technology required
early in the decade. A second study would meet later in the decade
to evaluate the science return, the technology development progress
and the outstanding technology risk, in time to make a recommendation
to Astro2030, likely to commence in 2028.

\section{Schedule}

The goal of this proposal is to have a bLISA mission concept ready
for consideration by Astro2030. To accomplish this goal, the technology
challenges for candidate mission concepts have to be understood and
surmountable in order to select the concept. Hence, the technology
development work needs to proceed during the 2020s, accompanied by
mission concept studies to identify the requisite technologies and
support a selection. LISA technology development will continue into
the mid-2020s, when LISA goes into Phase C. The pace of technology
development for bLISA should be set to support a concept down-select
no later than 2028. That down-select will need to appraise the relative
challenges of the candidate missions and the technology development
path to flight for each. That could require only a low TRL, say 3-4,
i.e., proof of concept. After that time, technology development for
bLISA should focus on advancing toward TRLs of 5-6 in anticipation
of a Phase A start early in the 2030s. If LISA launches in 2034, as
currently planned, it would be finishing extended operations in the
mid-2040s. Optimally, a bLISA would launch in the latter half of that
decade.

\section{Cost}

This is a (very) small proposal for space related activities. At this
time, the technologies to be developed are not well-enough understood
to produce a detailed roadmap with schedule and budget profile. At
best, we can only offer ROM costs. Based on our experience with LISA
technology development, we estimate that the technology development
needed in the decade of the 2020s will be about \$20M. The cost of
a mission concept study will be about \$2M. The study might reasonably
span two episodes, an early study to identify the candidate concepts
and their requisite technologies, and then a later study ending by
2028 that evaluates the development progress, the development remaining
and the potential science return of the candidate mission concepts
in order to make the down-select. Then a proposal would be made to
Astro2030.

\pagebreak{}


\section{Introduction}

Many sequential decision problems can be framed as the multi-armed bandit (MAB) problem \citep{Robbins, Lai85}, where a learner sequentially pulls arms and receives random rewards that possibly differ by arms. While the reward compensation mechanism is unknown, the learner can adapt his (her) decision to the past reward feedback so as to maximize the sum of rewards. Since the rewards of the unchosen arms remain unobserved, the learner should carefully balance between ``exploitation", pulling the best arm based on information accumulated so far, and ``exploration", pulling the arm that will assist in future choices, although it does not seem to be the best option at the moment. Application areas include the mobile healthcare system \citep{Tewari and Murphy}, web page ad placement algorithms \citep{Langford}, news article placement algorithms \citep{Li10}, revenue management \citep{Ferreira}, marketing \citep{Schwartz}, and recommendation systems \citep{Kawale}.

 Contextual MAB algorithms make use of side information, called context, given in the form of finite-dimensional covariates. For example, in the news article recommendation example, information on the visiting user as well as the articles are given in the form of a context vector. In 2010, the Yahoo! team \citep{Li10} proposed a contextual MAB algorithm which assumes a linear relationship between the context and reward of each arm. This algorithm achieved a 12.5\% click lift compared to a context-free MAB algorithm. Many other bandit algorithms have been proposed {for linear reward models} \citep{Auer,Dani, Chu11, Abbasi-Yadkori, Agrawal}.

 The aforementioned methods require the dimension of the context not be too large. This is because in these algorithms, the cumulative gap between the rewards of the optimal arm and the chosen arm, namely regret, is shown to be proportional to a polynomial function of the dimension of the context, $d$. In modern applications however, it is often the case that the web or mobile-based contextual variables are high-dimensional. \citet{Li10} applied a dimension reduction method \citep{Chu09} to the context variables before applying their bandit algorithm to the Yahoo! news article recommendation log data.

{In this paper, we introduce a novel approach to the case where the contextual variables are high-dimensional but only a sparse subset of size $s_0$ is correlated with the reward. We specifically consider the stochastic linear contextual bandit (LCB) problem, where the $N$ arms are represented by $N$ different contexts $b_1(t), \cdots, b_N(t)\in\mathbb{R}^d$ at a specific time $t$, and the rewards $r_1(t),\cdots, r_N(t)\in\mathbb{R}$ are controlled by a single linear regression parameter $\beta\in\mathbb{R}^d$, i.e., $r_i(t)$ has mean $b_i(t)^T\beta$ for $i=1,\cdots,N$. The LCB problem is distinguished from the contextual bandit problem with linear rewards (CBL), where the context $b(t)\in\mathbb{R}^d$ does not differ by arms but the the reward of the $i$-th arm is determined by an arm-specific parameter $\beta_i\in\mathbb{R}^d$, i.e., $r_i(t)$ has mean $b(t)^T\beta_i$ for $i=1,\cdots,N$.  }

When the number of arms is large, the CBL approach is not practical due to large number of parameters. Methods for CBL are also not suited to handle the case where the set of arms changes with time. When we recommend a news article or place an advertisement on the web page, the lists of news articles or advertisements change day by day. In this case, it would not be feasible to assign a different parameter for every new incoming item. Therefore, when the number of arms is large and the set of arms changes with time, LCB approaches including the proposed method can be applied while the CBL approaches cannot.

In supervised learning, Lasso \citep{Tibshirani} is a good tool for estimating the linear regression parameter when the covariates are high-dimensional but only a sparse subset is related to the outcome. However, the fast convergence property of Lasso is guaranteed when data are i.i.d. and when the observed covariates are not highly correlated, the latter referred to as the compatibility condition \citep{VanDeGeer}. In the {contextual bandit setting,} these conditions are often violated because the observations are adapted to the past and the context variables for which the rewards are observed converge to a small region of the whole context distribution as the learner updates its arm selection policy.

In the proposed method,  we resolve the non-compatibility issue by coalescing  the methods from missing data literature. We start from the observation that the
bandit problem is a missing data problem since the rewards for the arms that are not chosen are not observed, hence, missing.  The difference is that in missing data settings, missingness is controlled by the environment and given to the learner while in bandit settings, the missingness is controlled by the learner.  Since the learner controls the missingness, the missing mechanism, or arm selection  probability is known in bandit settings.  Given this probability, we can construct unbiased estimates of the rewards for the  arms  not chosen, using the doubly-robust technique \citep{Bang}, which allows  capitalizing on the context information corresponding to the arms that are not chosen.  These data are observed and available, however, were not utilized by {most of the existing contextual bandit algorithms.}

We propose the Doubly-Robust Lasso Bandit algorithm which hinges on the sparse structure as in Lasso \citep{Tibshirani}, but utilizes contexts for unchosen arms by blending the doubly-robust technique \citet{Bang} used in missing data literature. The high-probability upper bound of the regret incurred by the proposed algorithm has order $O(s_0\mathrm{log}(dT)\sqrt{T})$ where $T$ is the total number of time steps. We note that the regret does not depend on the number of arms, $N$, and scales with $\mathrm{log}(d)$ instead of a polynomial function of $d$. Therefore, the regret of the proposed algorithm is sublinear in $T$ even when $N\gg T$ and $d$ scales with $T$.

\citet{Abbasi-Yadkori12}, \citet{Carpentier} and \citet{Gilton} considered the sparse LCB problem as well.
  While the high-probability regret bound of \citet{Abbasi-Yadkori12} does not depend on $N$, it is proportional to $\sqrt{d}$ instead of $\mathrm{log}d$, so is not sublinear in $T$ when $d$ scales with $T$. { \citet{Carpentier} used an explicit exploration phase to identify the support of the regression parameter using techniques from compressed sensing. Their regret bound is tight scaling with $\mathrm{log}d$, but the algorithm is specific to the case where the set of arms is the unit ball for the $||\cdot||_2$ norm and fixed over time. \citet{Gilton} leveraged ideas from linear Thompson Sampling and Relevance Vector Machines \citep{Tipping}. The theoretical results of \citet{Gilton} are weak since they derived the regret bound under the assumption that a sufficiently small superset of the support for the regression parameter is known in advance. }

{\citet{Bastani} addressed the CBL problem with high-dimensional contexts. They proposed the Lasso Bandit algorithm which uses Lasso to estimate the parameter of each arm separately. To solve non-compatibility, \citet{Bastani} used forced-sampling of each arm at predefined time points and maintained two estimators for each arm, one based on the forced-samples and the other based on all samples. They derived a regret bound of order $O\big(Ns_0^2[\mathrm{log}T+\mathrm{log}d]^2\big)$.
{For the same problem}, \citet{Wang} proposed the Minimax Concave Penalized (MCP) Bandit algorithm which uses forced-sampling along with the MCP estimator \citep{Zhang} and improved the bound of \citet{Bastani} to $O\big(Ns_0^2[s_0+\mathrm{log}d]\mathrm{log}T\big)$.

The regrets of  \citet{Bastani} and \citet{Wang} increase linearly with $N$ due to separate estimation for each arm. When the number of arms $N$ is bigger than $T$,  these algorithms should terminate before the forced-sampling of each arm is completed, thus may not be practical. In contrast, the proposed method allows to share information across arms so its performance does not depend on $N$. In cases where all three methods are applicable, the proposed algorithm requires one less tuning parameter than Lasso Bandit and MCP Bandit, which is a significant advantage for an online learning algorithm as it is difficult to  simultaneously tune the hyperparameters and achieve high reward.}

We summarize the main contributions of the paper as follows.
\begin{itemize}
    \item We propose a new linear contextual MAB algorithm for high-dimensional, sparse reward models. The proposed method is simple and requires less tuning parameters than previous works.
    \item We propose a new estimator {for the regression parameter in the reward model} which uses Lasso with the context information of all arms through the doubly-robust technique.
    \item We prove that the high-probability regret upper bound of the proposed algorithm is $O(s_0\mathrm{log}(dT)\sqrt{T})$, which does not depend on $N$ and scales with $\mathrm{log}d$.
    \item {We present experimental results that show the superiority of our method especially when $N$ is big and the contexts of different arms are correlated.}
\end{itemize}

\subsection{Related Work}

{The fast convergence of Lasso {for linear regression parameter} was extenstively studied in \citet{VanDeGeer07}, \citet{Bickel}, and \citet{VanDeGeer}. The doubly-robust technique was originally introduced in \citet{Robins} and extensively analyzed in \citet{Bang} under supervised learning setup. \citet{Dudik}, \citet{Jiang}, and \citet{Farajtabar} are recent works that incorporate the doubly-robust methodology into reinforcement learning but under offline learning settings. \citet{Dimakopoulou} was the first to blend the doubly-robust technique in the online bandit problem. Their work, which focused on low-dimensional settings, proposed to use the doubly-robust technique as a way of balancing the data over the whole context space which makes the online regression less prone to bias when the reward model is misspecified.}

\section{Settings and Assumptions}

In the MAB setting, the learner is repeatedly faced with $N$ alternative arms where at time $t$, the $i$-th arm ($i=1,\cdots,N$) yields a random reward $r_i(t)$ with unknown mean $\theta_i(t)$. We assume that there is a finite-dimensional context vector $b_i(t)\in \mathbb{R}^d$ associated with each arm $i$ at time $t$ and that the mean of $r_i(t)$ depends on $b_i(t)$, i.e., $\theta_i(t)=\theta_t(b_i(t)),$
where $\theta_t(\cdot)$ is an arbitrary function.
Among the $N$ arms, the learner pulls one arm $a(t)$, and observes reward $r_{a(t)}(t)$. The optimal arm at time $t$ is $a^*(t):=\underset{1\leq i \leq N}{\mathrm{argmax}}\{\theta_t(b_i(t))\}$. Let $regret(t)$ be the difference between the expected reward of the optimal arm and the expected reward of the arm chosen by the learner at time $t$, i.e.,
\begin{align*}regret(t)&=\mathbb{E}\big(r_{a^*(t)}(t)-r_{a(t)}(t)~\big|~\{b_i(t)\}_{i=1}^N, a(t)~\big)\\&=\theta_t(b_{a^*(t)}(t))-\theta_t(b_{a(t)}(t)).\end{align*}
Then, the goal of the learner is to minimize the sum of regrets over $T$ steps,
$R(T):=\sum_{t=1}^Tregret(t).$

We specifically consider a sparse {LCB problem,} where $\theta_t(b_i(t))$ is linear in $b_i(t)$,
\begin{align}\theta_t(b_i(t))=b_i(t)^T\beta,~~~i=1,\cdots,N,\label{linear}\end{align}
where $\beta\in\mathbb{R}^d$ is unknown and $||\beta||_0=s_0(\ll d)$, $||\cdot||_p$ denoting the $L_p$-norm. {Hence, only $s_0$ elements in $\beta$ are assumed to be nonzero.} We also make the following assumptions, from \ref{bddness} to \ref{subgauss}.
\begin{assump}\label{bddness}{\bf Bounded norms.}
Without loss of generality, $||b_i(t)||_2\leq 1,~ ||\beta||_2\leq 1.$
\end{assump}

\begin{assump}\label{iid}{\bf IID Context Assumption.}  The distribution of context variables is i.i.d. over time $t$:
\begin{align*}
\{b_1(t), \cdots, b_N(t)\} \overset{i.i.d.}{\sim} \mathcal{P}_b,
\end{align*}
where $\mathcal{P}_b$ is some distribution over $\mathbb{R}^{N\times d}$.
\end{assump}
We note in  \ref{iid} that given time $t$, the contexts from different arms are allowed to be correlated.
\begin{assump}\label{compatibility} {\bf Compatibility Assumption \citep{VanDeGeer}.}  Let $I$ be a set of indices and $\phi$ be a positive constant. Define $$C(I,\phi)=\{M\in\mathbb{R}^{d\times d}: \forall v \in \mathbb{R}^d \text{ such that } ||v_{I^C}||_1\leq 3||v_I||_1, ||v_I||_1^2\leq |I|(v^TMv)/\phi^2\}.$$
Then $\exists \phi_1>0$ such that $$\Sigma:=\mathbb{E}[\bar{b}(t)\bar{b}(t)^T]\in C(supp(\beta), \phi_1),$$
where $\bar{b}(t)=\frac{1}{N}\sum_{i=1}^Nb_i(t),$ and $supp(\beta)$ denotes the support of $\beta$.
\end{assump}
 \ref{compatibility} ensures that the Lasso estimate of $\beta$ converges to the true parameter in a fast rate. This assumption is weaker than the restricted eigenvalue assumption \citep{Bickel}.
\begin{assump}\label{subgauss}{\bf Sub-Gaussian error.}
The error $\eta_i(t):=r_i(t)-b_i(t)^T\beta$ is $R$-sub-Gaussian for some $0<R<O(\sqrt[4]{\mathrm{log}T/T})$, i.e., for every $\lambda\in \mathbb{R},$
\begin{align*}\mathbb{E}[\mathrm{exp}(\lambda\eta_i(t))]\leq \mathrm{exp}({\lambda^2R^2}/{2}).\end{align*}
\end{assump}
Assumption \ref{subgauss} is satisfied whenever $r_i(t)\in [b_i(t)^T\beta-R,b_i(t)^T\beta+R]$.

\section{Doubly-Robust Lasso Bandit}

Lemma 11.2 of \citet{VanDeGeer} shows that when the covariates satisfy the compatibility condition and the noise is i.i.d. Gaussian, the Lasso estimator of the linear regression parameter converges to the true parameter in a fast rate. We restate the lemma.
\begin{lem}\label{lassolemm}{\bf \citep[Lemma 11.2 of ][]{VanDeGeer}} Let $x_{\tau}\in\mathbb{R}^d$ and $y_{\tau}\in\mathbb{R}$ be random variables with $y_{\tau}=x_{\tau}^T\beta+\varepsilon_{\tau},~\tau=1,2,\cdots,t,$ where $\beta\in\mathbb{R}^d$, $||\beta||_0=s_0$, and $\varepsilon_{\tau}$'s are i.i.d. Gaussian with mean zero and variance $R^2$. Let $\hat{\beta}(t)$ be the Lasso estimator of $\beta$ based on the first $t$ observations using $\lambda_t=R\sqrt{\frac{2\mathrm{log}(\mathrm{e}d/\delta)}{t}}$ for the regularization parameter in the $L_1$ penalty, i.e.,
$$\hat{\beta}(t)=\underset{\beta}{\mathrm{argmin}}\big\{\frac{1}{t}\sum_{\tau=1}^t(y_{\tau}-x_{\tau}^T\beta)^2+\lambda_t||\beta||_1\big\}.$$
If $\hat{\Sigma}_t:=\frac{1}{t}\sum_{\tau=1}^t x_{\tau}x_{\tau}^T\in C(supp(\beta),\phi)$ for some $\phi>0$, then for $\forall \delta\in (0,1)$, with probability at least $(1-\delta)$,
\begin{align}||\hat{\beta}(t)-\beta||_1\leq \frac{4\lambda_t s_0}{\phi^2}=\frac{4s_0R}{\phi^2}\sqrt{\frac{2\mathrm{log}(\mathrm{e}d/\delta)}{t}}.\label{iidfast}\end{align}
\end{lem}

Two hurdles that arise when incorporating Lemma \ref{lassolemm} into the contextual MAB setting are that (a) {(non-i.i.d.)} the errors $\varepsilon_{\tau}$'s are not i.i.d., and (b) {(non-compatibility)} the Gram matrix $\hat{\Sigma}_t$ does not satisfy the compatibility condition even under assumption \ref{compatibility} because the context variables of which the rewards are observed do not evenly represent the whole distribution of the context variables. In the {LCB setting,} $x_{\tau}$ corresponds to $b_{a(\tau)}(\tau)$. Since the learner tends to choose the contexts that yield maximum reward as time elapses, the chosen context variables $b_{a(\tau)}(\tau)$'s converge to a small region of the whole context space. We discuss on remedies to (a) and (b) in Section \ref{lassoMAB} and Section \ref{pseudor}, respectively.

\subsection{Lasso Oracle Inequality with non-i.i.d. noise}\label{lassoMAB}

As a remedy to {the non-i.i.d. problem,} \citet{Bastani} proposed a variation of Lemma \ref{lassolemm} which shows that (\ref{iidfast}) holds when $\{\varepsilon_{\tau}\}_{\tau=1}^t$ is a martingale difference sequence. We restate the proposition of \citet{Bastani}.
\begin{lem}\label{lassolemm2}{\bf \citep[Proposition 1 of ][]{Bastani}} Let $\mathcal{F}_{\tau-1}$ denote the filtration up to time $\tau-1$ in the bandit setting, $$\mathcal{F}_{\tau-1}=\{x_1,y_1,\cdots,x_{\tau-1},y_{\tau-1},x_{\tau}\},~~~\tau=1,2,\cdots,t.$$ Suppose the conditions of Lemma \ref{lassolemm} hold except that $\varepsilon_{\tau}$'s are i.i.d. Suppose instead that $\varepsilon_{\tau}|\mathcal{F}_{\tau-1}$ is $R$-sub-Gaussian for $\tau=1,\cdots,t$. Then for $\forall \delta\in (0,1)$, with probability at least $(1-\delta)$,
\begin{align}||\hat{\beta}(t)-\beta||_1\leq \frac{4\lambda_t s_0}{\phi^2}=\frac{4s_0R}{\phi^2}\sqrt{\frac{2\mathrm{log}(\mathrm{e}d/\delta)}{t}}.\label{mtgfast}\end{align}
\end{lem}

\subsection{Doubly-robust pseudo-reward}\label{pseudor}

To overcome the non-compatibility problem, {in the CBL setting,} \citet{Bastani} proposed to impose forced-sampling of each arm at predefined $O(\mathrm{log}T)$ time steps, which produces i.i.d. data. The following lemma shows that when $x_{\tau}$'s are i.i.d. and $\mathbb{E}(x_{\tau}x_{\tau}^T)\in C(supp(\beta),\phi)$, then $\hat{\Sigma}_t\in C(supp(\beta),\phi/\sqrt{2})$ with high probability.
\begin{lem}\label{iidcomp}{\bf \citep[Lemma EC.6. of ][]{Bastani}}
Let $x_1, x_2, \cdots, x_t$ be i.i.d. random vectors in $\mathbb{R}^d$ with $||x_{\tau}||_{\infty}\leq 1$ for all $\tau$. Let $\Sigma=\mathbb{E}[x_{\tau}x_{\tau}^T]$ and $\hat{\Sigma}_t=\frac{1}{t}\sum_{\tau=1}^tx_{\tau}x_{\tau}^T.$ Suppose that $\Sigma\in C(supp(\beta),\phi)$. Then if $t\geq \frac{3}{c^2}\mathrm{log}d$,
\begin{align}\mathbb{P}\Big[\hat{\Sigma}_t\in C\Big(supp(\beta),\frac{\phi}{\sqrt{2}}\Big)\Big]\geq 1-\mathrm{exp}\big(-c^2t\big),\label{iidcompineq}\end{align}
where $c=\mathrm{min}(0.5, \frac{\phi^2}{256s_0})$.
\end{lem}
We derive the following corollary simply by setting the right-hand side of (\ref{iidcompineq}) larger than $1-\frac{\delta'}{t^2}$ .
\begin{cor}\label{cor1} Suppose that the conditions of Lemma \ref{iidcomp} are satisfied. Let $z_T=\mathrm{max}\big(\frac{3}{c^2}\mathrm{log}d, \frac{1}{c^2}\mathrm{log}\frac{T^2}{\delta'}\big)$ for some $\delta'\in (0,1)$. Then with probability at least $1-\delta'$,  \begin{align}
 \hat{\Sigma}_t\in C\Big(supp(\beta),\frac{\phi}{\sqrt{2}}\Big) \text{~~for all~} t\geq z_T. \label{iidcomp2}
\end{align}
\end{cor}
The setting of \citet{Bastani} is different from ours in that they assume the context variable is the same for all arms but the regression parameter differs by arms. \citet{Bastani} maintained two sets of estimators for each arm, the estimator based on the forced-samples and the one based on all samples. Whereas the latter is not based on i.i.d. data, using the forced-sample estimator as a pre-processing step of selecting a subset of arms and then using the all-sample estimator to select the best arm among this subset guarantees that there are $O(T)$ i.i.d. samples for each arm, so the all-sample estimator converges in a fast rate as well.

We propose a different approach to resolve non-compatibility, which is based on the doubly-robust technique in missing data literature. We define the filtration $\mathcal{F}_{t-1}$ as the union of the observations until time $t-1$ and the contexts given at time $t$, i.e.,  $$\mathcal{F}_{t-1}=\{\{b_i(1)\}_{i=1}^N, a(1), r_{a(1)}(1),\cdots,\{b_i(t-1)\}_{i=1}^N, a(t-1), r_{a(t-1)}(t-1),\{b_i(t)\}_{i=1}^N\}.$$ Given $\mathcal{F}_{t-1}$, we pull the arm $a(t)$ randomly according to probability $\pi(t)=[\pi_1(t),\cdots,\pi_N(t)]$, where $\pi_i(t)=\mathbb{P}[a(t)=i|\mathcal{F}_{t-1}]$ is the probability of pulling the $i$-th arm at time $t$ given $\mathcal{F}_{t-1}$. {We specify the values of $\pi$(t) later in Section \ref{sub3}.} Let $\hat{\beta}({t-1})$ be the $\beta$ estimate given $\mathcal{F}_{t-1}$. After the reward $r_{a(t)}(t)$ is observed, we construct a doubly-robust pseudo-reward:
$$\hat{r}(t)=\bar{b}(t)^T\hat{\beta}(t-1)+\frac{1}{N}\frac{r_{a(t)}(t)-b_{a(t)}(t)^T\hat{\beta}(t-1)}{\pi_{a(t)}(t)}.$$
Whether or not $\hat{\beta}(t-1)$ is a valid estimate, this value has conditional expectation $\bar{b}(t)^T\beta$ given that $\pi_i(t)>0$ for all $i$:
\begin{align*}
\mathbb{E}[\hat{r}(t)|\mathcal{F}_{t-1}]&=\mathbb{E}\Big[\frac{1}{N}\sum_{i=1}^N\Big(1-\frac{I(a(t)=i)}{\pi_i(t)}\Big)b_i(t)^T\hat{\beta}(t-1)+\frac{1}{N}\sum_{i=1}^N\frac{I(a(t)=i)}{\pi_i(t)}r_i(t)\Big|\mathcal{F}_{t-1}\Big]\\
&=\mathbb{E}\Big[\frac{1}{N}\sum_{i=1}^Nr_i(t)\Big|\mathcal{F}_{t-1}\Big]=\bar{b}(t)^T\beta.
\end{align*}
Thus by weighting the observed reward $r_{a(t)}(t)$ with the inverse of its observation probability $\pi_{a(t)}(t)$, we obtain an unbiased estimate of the reward corresponding to the average context $\bar{b}(t)$ instead of $b_{a(t)}(t)$. Applying Lasso regression to the pair $(\bar{b}(t),\hat{r}(t))$ instead of $(b_{a(t)}(t), r_{a(t)}(t))$, we can make use of the i.i.d. assumption \ref{iid} and the compatibility assumption  \ref{compatibility} with Corollary \ref{cor1} to achieve a fast convergence rate for $\beta$ as in (\ref{mtgfast}). {We later show that when $x_{\tau}$ in Corollary \ref{cor1} corresponds to $\bar{b}(\tau)$, (\ref{iidcomp2}) holds with $z_T=O(s_0^2\mathrm{log}(dT))$.}

{Since the allocation probability $\pi_{a(t)}(t)$ is known given $\mathcal{F}_{t-1}$, the simple inverse probability weighting (IPW) estimator $\breve{r}(t):=\frac{r_{a(t)}(t)}{N\pi_{a(t)}(t)}$ also gives an unbiased estimate for the reward corresponding to the average context. We however show in the next paragraphs that the doubly-robust estimator has constant variance under minor condition on the weight $\pi_{a(t)}(t)$ while the IPW estimator has variance that increases with $t$ under the same condition. Constant variance is crucial for the performance of the algorithm, since it affects the convergence property of $\hat{\beta}(t)$ in (\ref{mtgfast}) and eventually the regret upper bound in Theorem \ref{thm1}. We can make the IPW estimator have constant variance as well under stronger condition on $\pi_{a(t)}(t)$, but this stronger condition is shown to result in regret linear in $T$.

Let $\tilde{R}^2$ be the variance of $\hat{r}(t)$ given $\mathcal{F}_{t-1}$. We have,
\begin{align}
\tilde{R}^2&:=Var[\hat{r}(t)|\mathcal{F}_{t-1}]=\mathbb{E}\big[\big\{\hat{r}(t)-\bar{b}(t)^T\beta\big\}^2\big|\mathcal{F}_{t-1}\big]\nonumber\\
&=\mathbb{E}\Big[\Big\{\bar{b}(t)^T\big(\hat{\beta}(t-1)-\beta\big)+\frac{\eta_{a(t)}(t)}{N\pi_{a(t)}(t)}+\frac{b_{a(t)}(t)^T(\beta-\hat{\beta}(t-1))}{N\pi_{a(t)}(t)}\Big\}^2\Big|\mathcal{F}_{t-1}\Big].
\label{var}\end{align}
{Suppose $\hat{\beta}(t-1)$ satisfies the Lasso convergence property (\ref{mtgfast}), i.e.,} $||\hat{\beta}(t-1)-\beta||_1\leq O(\sqrt{(\mathrm{log}d+\mathrm{log}
t)/{t}})$ for $t > z_T$.   Then we can bound the first and third terms of (\ref{var}) by a constant under the following restriction,
{\begin{equation}\pi_{a(t)}(t)=\begin{cases}\frac{1}{N} &\text{ for all $t\leq z_T$}\\ O\big(\frac{1}{N}\sqrt{(\mathrm{log}d+\mathrm{log}t)/{t}}\big) &\text{ for all $t>z_T$.}  \end{cases}\label{pi_cond}\end{equation}}
Due to assumption \ref{subgauss} on the sub-Gaussian error $\eta$, the second term of (\ref{var}) is also bounded by a constant under (\ref{pi_cond}). Hence, we have $\tilde{R}^2 \leq O(1)$. Therefore if  (\ref{pi_cond}) holds and if $\hat{\beta}(z_T), \hat{\beta}(z_T+1), \cdots, \hat{\beta}(t-1)$ all satisfy the Lasso convergence property (\ref{mtgfast}), then the pseudo-rewards $\hat{r}(1),\cdots,\hat{r}(t)$ all have constant variance. Consequently, the estimate $\hat{\beta}(t)$ based on $\{\bar{b}(\tau),\hat{r}(\tau)\}_{\tau=1}^t$ satisfies (\ref{mtgfast}). {Meanwhile, the restriction (\ref{pi_cond}) leads to suboptimal choice of arms at each $t$ with probability that decreases with time. In Theorem \ref{thm1}, we prove that the regret due to this suboptimal choice is sublinear in $T$. }

The variance of $\breve{r}(t)$ given filtration $\mathcal{F}_{t-1}$ is $\mathbb{E}\Big[\Big\{\frac{\eta_{a(t)}(t)}{N\pi_{a(t)}(t)}+\frac{b_{a(t)}(t)^T\beta}{N\pi_{a(t)}(t)}-\bar{b}(t)^T\beta\Big\}^2\Big|\mathcal{F}_{t-1}\Big]$.

As opposed to the first and third terms in (\ref{var}), the term ${b_{a(t)}(t)^T\beta}/{N\pi_{a(t)}(t)}$ still increases with $t$ under (\ref{pi_cond}). To achieve a constant variance, we need $\pi_{a(t)}(t)$ be larger than a predetermined constant value, $\frac{1}{N}p_{min}$. {Simply replacing $\pi_{a(t)}(t)$ in $\breve{r}(t)$ with the truncated value $\breve{\pi}_{a(t)}(t)=\mathrm{max}\big\{\frac{1}{N}p_{min}, \pi_{a(t)}(t)\big\}$ will produce a biased estimate of $\bar{b}(t)^T\beta$ when $\pi_{a(t)}(t)$ is actually smaller than $\frac{1}{N}p_{min}$, and the resulting estimate $\hat{\beta}(t)$ will not satisfy convergence property (\ref{mtgfast}).} If we instead directly constrain $\pi_{a(t)}(t)$ to be larger than $\frac{1}{N}p_{min}$, this will lead to suboptimal choice of arms at each $t$ with constant probability so the regret will increase linearly in $T$.}

\subsection{Doubly-Robust Lasso Bandit algorithm}\label{sub3}

To make $\pi_{a(t)}(t)=\frac{1}{N}$ for all $t\leq z_T$, we simply pull arms according to the uniform distribution when $t\leq z_T$. Then to ensure $\pi_{a(t)}(t)\geq O\big(\frac{1}{N}\sqrt{(\mathrm{log}d+\mathrm{log}t)/{t}}\big)$ for all $t>z_T$,  we randomize the arm selection at each step between uniform random selection
and deterministic selection using the $\beta$ estimate of the last step, $\hat{\beta}(t-1).$ This can be implemented  in two stages. First, sample $m_t$ from a Bernouilli distribution with mean $\lambda_{1t}=\lambda_1\sqrt{(\mathrm{log}t+\mathrm{log}d)/t},$ where $\lambda_1>0$ is a tuning parameter. Then if $m_t=1$, pull the arm $a(t)=i$ with probability $1/N$, otherwise, pull the arm $a(t)=\underset{1\leq i\leq N}{\mathrm{argmax}}\{b_i(t)^T\hat{\beta}(t-1)\}$. Under this procedure, we have $\pi_{a(t)}(t)=\lambda_{1t}/N+(1-\lambda_{1t})I\big\{a(t)=\underset{1\leq i\leq N}{\mathrm{argmax}}\{b_i(t)^T\hat{\beta}(t-1)\}\big\},$ which satisfies (\ref{pi_cond}). In practice, we treat $z_T$ as a tuning parameter.

The proposed Doubly-Robust (DR) Lasso Bandit algorithm is presented in Algorithm \ref{prop_alg}. The algorithm selects $a(t)$ via two stages, computes $\pi_{a(t)}(t)$, constructs the pseudo-reward $\hat{r}(t)$ and average context $\bar{b}(t)$, and updates the $\beta$ estimate using Lasso regression on the pairs $(\bar{b}(t),\hat{r}(t))$'s. Recall that the regularization parameter should be set as $\lambda_{2t}=O\big(\sqrt{{\mathrm{log}(d/\delta)}/{t}}\big)$ to guarantee (\ref{mtgfast}) for a fixed $t$. To guarantee (\ref{mtgfast}) for every $t>z_T$, we impute $\delta/t^2$ instead of $\delta$. Hence at time $t$, we update $\lambda_{2t}=\lambda_2\sqrt{(\mathrm{log}t+\mathrm{log}d)/t}$ where $\lambda_2>0$ is a tuning parameter. The algorithm requires three tuning parameters in total, while \citet{Bastani} and \citet{Wang} require four.

\begin{algorithm}
\caption{DR Lasso Bandit}\label{prop_alg}
\begin{algorithmic}
\State Input parameters: $\lambda_1$, $\lambda_2$, $z_T$
\State Set $\hat{\beta}(0)=0_d$, $\mathbb{S}=\{~\}$.
\For{$t=1,2,\cdots, T$}
\State Observe $\{b_1(t), b_2(t), \cdots, b_N(t)\}\sim \mathcal{P}_b$
\If{$t\leq z_T$}
\State Pull arm $a(t)=i$ with probability $\frac{1}{N}$ $(i=1,\cdots,N)$
\State $\pi_{a(t)}(t)\leftarrow 1/{N}$
\Else
\State $\lambda_{1t}\leftarrow \lambda_1\sqrt{(\mathrm{log}t+\mathrm{log}d)/t}$, sample $m_t\sim Ber(\lambda_{1t})$
\If{$m_t=1$}
\State Pull arm $a(t)=i$ with probability $\frac{1}{N}$ $(i=1,\cdots,N)$
\Else
\State Pull arm $a(t)=\underset{1\leq i\leq N}{\mathrm{argmax}}\{b_i(t)^T\hat{\beta}(t-1)\}$
\EndIf
\State $\pi_{a(t)}(t)\leftarrow \lambda_{1t}/{N}+(1-\lambda_{1t})I\Big\{a(t)=\underset{1\leq i\leq N}{\mathrm{argmax}}\{b_i(t)^T\hat{\beta}(t-1)\}\Big\}$
\EndIf
\State Observe $r_{a(t)}(t)$
\State $\bar{b}(t)\leftarrow\frac{1}{N}\sum_{i=1}^Nb_i(t)$,~~ $\hat{r}(t)\leftarrow\bar{b}(t)^T\hat{\beta}(t-1)+\frac{1}{N}\frac{r_{a(t)}(t)-b_{a(t)}(t)^T\hat{\beta}(t-1)}{\pi_{a(t)}(t)}$
\State $\mathbb{S}\leftarrow \mathbb{S} \cup \big\{\big(\bar{b}(t), \hat{r}(t)\big)\big\}$
\State $\lambda_{2t}\leftarrow \lambda_2\sqrt{(\mathrm{log}t+\mathrm{log}d)/t}$
\State $\hat{\beta}(t)\leftarrow\underset{\beta}{\mathrm{argmin}}\Big\{\frac{1}{t}\sum_{(\bar{b},\hat{r})\in \mathbb{S}}(\hat{r}-\bar{b}^T\beta)^2+\lambda_{2t}||\beta||_1\Big\}$
\EndFor
\end{algorithmic}
\end{algorithm}

\section{Regret analysis}

Under (\ref{linear}), the regret at time $t$ is $regret(t)=b_{a^*(t)}(t)^T\beta-b_{a(t)}(t)^T\beta,$
where $a^*(t)=\underset{1\leq i \leq N}{\mathrm{argmax}}\{b_i(t)^T\beta\}.$
We establish the high-probability regret upper bound for Algorithm \ref{prop_alg} in the following theorem.

\begin{thm}\label{thm1} Suppose (\ref{linear}) and assumptions \ref{bddness}, \ref{iid}, \ref{compatibility}, and \ref{subgauss} hold. Then we have for $\forall\delta\in (0,\frac{1}{2})$ and $\forall \delta'\in (0,\delta)$, with probability at least $1-2\delta$, \begin{align*}R(T)&\leq \mathrm{max}\Big(\frac{3}{C(\phi_1)^2}\mathrm{log}d, \frac{1}{C(\phi_1)^2}{\mathrm{log}(\frac{T^2}{\delta'})}\Big)+\frac{\sqrt{128}}{\phi_1^2}s_0\tilde{R}\sqrt{T}\mathrm{log}\Big(\frac{\mathrm{e}dT^2}{\delta-\delta'}\Big)+\lambda_1\sqrt{T}\mathrm{log}(dT)\\
&+\sqrt{\frac{T}{2}\mathrm{log}\big(\frac{1}{\delta}\big)}\\
&=O\big(s_0\sqrt{T}\mathrm{log}(dT)\big), \end{align*}
where $C(\phi_1)=\mathrm{min}(0.5, \frac{\phi_1^2}{256s_0})$ and $\tilde{R}=O(1)$.
\end{thm}

{We provide a sketch of proof in the next section. A complete proof is presented in the Supplementary Material.}

\subsection{Sketch of proof for Theorem \ref{thm1}}
{
In the DR Lasso Bandit algorithm, each of the $T$ decision steps corresponds to one of the following three groups.
\begin{enumerate}[(a)]
    \item  $t\leq z_T~$: the arms are pulled according to the uniform distribution.
    \item  $t>z_T$ and $m_t=1~$: the arms are pulled according to the uniform distribution.
    \item  $t>z_T$ and $m_t=0~$: the arm $a(t)=\underset{1\leq i\leq N}{\mathrm{argmax}}\{b_i(t)^T\hat{\beta}(t-1)\}$ is pulled.
\end{enumerate}
We let $z_T=\mathrm{max}\Big(\frac{3}{C(\phi_1)^2}\mathrm{log}d, \frac{1}{C(\phi_1)^2}{\mathrm{log}(\frac{T^2}{\delta'})}\Big)$. The strategy of the proof is to bound the regrets from each separate group and sum the results.

Due to assumption \ref{bddness} on the norms of $b_i(t)$ and $\beta,$ $regret(t)$ is not larger than $1$ in any case. Therefore, the sum of regrets from group (a) is at most $z_T,$ which corresponds to the first term in Theorem \ref{thm1}. We now denote the sum of regrets from group (b) and group (c) as $R(T,b)$ and $R(T,c),$ respectively. We first bound $R(T,b)$ in Lemma \ref{Rtb}, which follows from the fact that $m_t$ is a Bernouilli variable with mean $\lambda_{1t}$ and  Hoeffding's inequality.

\begin{lem}\label{Rtb} With probability at least $1-\delta,$
$$R(T,b)\leq \sum_{t=z_T}^T\lambda_{1t}+\sqrt{T\mathrm{log}(1/\delta)/2} \leq \lambda_1\sqrt{T}\sqrt{\mathrm{log}(dT)}\sqrt{1+\mathrm{log}T}+ \sqrt{T\mathrm{log}({1}/{\delta})/2}.$$
\end{lem}
To bound $R(T,c)$, we further decompose group (c) into the following two subgroups, where we set $d_t=\frac{\sqrt{128}}{\phi_1^2}s_0\tilde{R}\sqrt{\frac{\mathrm{log}(\mathrm{e}dt^2/(\delta-\delta'))}{t}}$ with $\tilde{R}=O(1)$.

(c-a) $t>z_T$, $m_t=0$~: $a(t)=\underset{1\leq i\leq N}{\mathrm{argmax}}\{b_i(t)^T\hat{\beta}(t-1)\}$ is pulled and $||\hat{\beta}(t-1)-\beta||_1\leq d_{t-1}.$

(c-b) $t>z_T$, $m_t=0$~: $a(t)=\underset{1\leq i\leq N}{\mathrm{argmax}}\{b_i(t)^T\hat{\beta}(t-1)\}$ is pulled and $||\hat{\beta}(t-1)-\beta||_1>d_{t-1}.$

We first show in Lemma \ref{lassolemm3} that the subgroup (c-b) is empty with high probability. Since $\hat{\beta}(t)$ is computed from the pairs $\{\bar{b}(\tau), \hat{r}(\tau)\}_{\tau=1}^t$ where the average contexts $\bar{b}(\tau)$'s satisfy the compatibility condition (\ref{iidcomp2}) and the pseudo-rewards $\hat{r}(\tau)$'s are unbiased with constant variance, Lemma \ref{lassolemm2} can be applied on $\hat{\beta}(t)$.
\begin{lem}\label{lassolemm3}
Suppose (\ref{linear}) and assumptions \ref{bddness}, \ref{iid}, \ref{compatibility}, and \ref{subgauss} hold. Let $d_t=\frac{\sqrt{128}}{\phi_1^2}s_0\tilde{R}\sqrt{\frac{\mathrm{log}(\mathrm{e}dt^2/(\delta-\delta'))}{t}}$. Then we have for $\forall \delta\in (0,1)$, for every $t\geq z_T$,
\begin{align}\mathbb{P}\Big(||\hat{\beta}(t)-\beta||_1 > d_{t} ~\Big|~ ||\hat{\beta}(t-1)-\beta||_1 < d_{t-1}, \cdots, ||\hat{\beta}(z_T)-\beta||_1 < d_{z_T}\Big)\leq \frac{\delta}{t^2}.\label{ineqkey}\end{align}
Hence, with probability at least $1-\delta,$
$$||\hat{\beta}(t)-\beta||_1<d_t~\text{ for every }t\geq z_T.$$
\end{lem}
Finally when $t$ belongs to group (c-a), $regret(t)$ is trivially bounded by $d_t.$ Hence we can bound $R(T,c)$ as in the following lemma.
\begin{lem}With probability at least $1-\delta,$
$$R(T,c)\leq \sum_{t=z_T}^Td_t\leq \frac{\sqrt{128}}{\phi_1^2}s_0\tilde{R}\sqrt{T}\mathrm{log}\Big(\frac{\mathrm{e}dT^2}{\delta-\delta'}\Big).$$
\end{lem}
}

\section{Simulation study}

\begin{figure}
\begin{center}
\includegraphics[width=13.8cm]{plot3.png}
\end{center}
\caption{Median (solid), 1st and 3rd quartiles (dashed) of $R(t)$ over 10 experiments.}\label{fig1}
\end{figure}

We conduct simulation studies to evaluate the proposed DR Lasso Bandit and the Lasso bandit  \citep{Bastani}. We set $N\!\!=\!\!10,$ $20,$ $50,$ or $100$, $d\!=\!100$, and $s_0=5$. For fixed $j=1,\cdots,d,$ we generate $[b_{1j}(t),\cdots,b_{Nj}(t)]^T$ from $\mathcal{N}(0_N,V)$ where $V(i,i)=1$ for every $i$ and $V(i,k)=\rho^2$ for every $i \neq k$. We experiment two cases for $\rho^2$, either $\rho^2=0.3$ (weak correlation) or $\rho^2=0.7$ (strong correlation). We generate  $\eta_i(t)\overset{i.i.d.}{\sim}\mathcal{N}(0,0.05^2)$ and the rewards from (\ref{linear}). We set $||\beta||_0=s_0$ and generate the $s_0$ non-zero elements from a uniform distribution on $[0,1]$.

We conduct 10 replications for each case. The Lasso Bandit algorithm can be applied in our setting by considering a $Nd$-dimensional context vector $b(t)=[b_1(t)^T,\cdots,b_N(t)^T]^T$  and a $Nd$-dimensional regression parameter $\beta_i$ for each arm $i$ where $\beta_i=[I(i=1)\beta^T,\cdots,I(i=N)\beta^T]^T$. For each algorithm, we consider some candidates for the tuning parameters and report the best results. For DR Lasso Bandit, we advise to truncate the value $\hat{r}(t)$  so that it does not explode.

Figure \ref{fig1} shows the cumulative regret $R(t)$ according to time $t$. {When the number of arms $N$ is as small as 10, Lasso Bandit converges faster to the optimal decision rule than DR Lasso Bandit. However, we notice that the regret of Lasso Bandit in the early stage is larger than DR Lasso Bandit, which is due to forced-sampling of arms.} With larger number of arms, DR Lasso Bandit outperforms Lasso Bandit. {Compared to Lasso Bandit,} the regret of DR Lasso Bandit decreases dramatically as the correlation among the contexts of different arms increases. We also verify that the performance of the proposed method is not sensitive to $N$, while the regret of Lasso Bandit increases with $N$.

\section{Concluding remarks}

Sequential decision problems involving web or mobile-based data call for contextual MAB methods that can handle high-dimensional covariates. In this paper, we propose a novel algorithm which enjoys the convergence properties of the Lasso estimator for i.i.d. data via a doubly-robust technique established in the missing data literature. The proposed algorithm attains a tight high-probability regret upper bound which depends on a polylogarithmic function of the dimension and does not depend on the number of arms, overcoming weaknesses of the existing algorithms.