[ beforeskip=.3em plus .6pt,pagenumberformat=]toclinesection

# Time-Uniform Confidence Spheres
for Means of Random Vectors

Ben Chugg  Machine Learning Department  Hongjian Wang  Department of
Statistics and Data Science  Aaditya Ramdas  Machine Learning Department
Department of Statistics and Data Science

(November 2023)

###### Abstract

We derive and study time-uniform confidence spheres—termed _confidence sphere
sequences_ (CSSs)—which contain the mean of random vectors with high
probability simultaneously across all sample sizes. Inspired by the original
work of Catoni and Giulini, we unify and extend their analysis to cover both
the sequential setting and to handle a variety of distributional assumptions.
More concretely, our results include an empirical-Bernstein CSS for bounded
random vectors (resulting in a novel empirical-Bernstein confidence interval),
a CSS for sub-$\psi$ random vectors, and a CSS for heavy-tailed random vectors
based on a sequentially valid Catoni-Giulini estimator. Finally, we provide a
version of our empirical-Bernstein CSS that is robust to contamination by
Huber noise.

###### Contents

  1. 1 Introduction
    1. 1.1 Related work
    2. 1.2 Background and approach
    3. 1.3 Assumptions
  2. 2 Light-tailed random vectors
    1. 2.1 Empirical-Bernstein bound
    2. 2.2 Sub-$\psi$ bound
  3. 3 Heavy-tailed random vectors
  4. 4 Robust confidence spheres
  5. 5 Discussion
    1. 5.1 Adaptivity to anisotropy via ellipsoids
    2. 5.2 Dependence on $d$
  6. 6 Summary
  7. A Stitched bounds with LIL Rate
    1. A.1 Stitched Empirical-Bernstein Bound
    2. A.2 Stitched Sub-Gamma Bound
  8. B A second (robust) time-uniform Catoni-Giulini bound
  9. C A semi-empirical bound for heavy-tails
  10. D Omitted proofs
    1. D.1 Section 2
    2. D.2 Section 3
    3. D.3 Section 4
  11. E Von Mises-Fisher distribution
  12. F Asymptotic width of empirical-Bernstein bounds
    1. F.1 Non-robust case
    2. F.2 Robust case

##  1 Introduction

We consider the classical problem of estimating the mean of random vectors in
nonparametric settings. Perhaps unlike traditional approaches however, we
would like to perform this task sequentially. That is, we are interested
providing estimates of the mean which hold uniformly across time as more data
arrives. To formalize, let $(\boldsymbol{X}_{t})_{t\geqslant 1}$ be a stream
of iid random vectors in $\mathcal{X}\subset\mathbb{R}^{d}$, $d\geqslant 2$,
with mean $\boldsymbol{\mu}$ drawn from some unknown distribution $P$. Our
goal is to construct a non-trivial sequence of sets,
$(\boldsymbol{C}_{t})_{t\geqslant 1}$, $\boldsymbol{C}_{t}\in\mathbb{R}^{d}$,
such that

$P(\forall t\geqslant 1:\boldsymbol{\mu}\in\boldsymbol{C}_{t})\geqslant
1-\alpha,$ (1)

for some predetermined $\alpha\in(0,1)$. In $\mathbb{R}^{1}$, such a sequence
is often called a _confidence sequence_. We will use the more general term
_confidence sphere sequence_ (CSS) to highlight that we are working in higher
dimensions, and because our sets $\boldsymbol{C}_{t}$ are spherical as opposed
to, say, ellipsoidal. Confidence (sphere) sequences are in contradistinction
to confidence _intervals_ —the latter hold only at fixed times $t=n$
determined prior to receiving any samples, while the former hold
_simultaneously_ at all stopping times. Consequently, they are often described
as _time-uniform_ or as _anytime-valid_. Importantly, confidence (sphere)
sequences obviate the need to deploy union bounds to gain uniformity over the
sample size. This is attractive for tasks in which data are continuously
monitored during collection and which may call for data-dependent stopping
times, two properties disallowed by fixed-time constructions such as
confidence intervals (where it is known as p-hacking).

The study of confidence sequences was pioneered by Darling, Robbins, Siegmund
and Lai in the decade from 1967-76 (Darling and Robbins, 1967a; Lai, 1976).
While further inquiry laid somewhat dormant for several decades afterwards,
they have recently regained interest due to novel mathematical techniques and
applicability to sequential decision-making. Indeed, confidence sequences are
used by several companies for A/B testing and adaptive
experimentation.111E.g., see Adobe (2023), Evidently (2023), and Netflix
(Lindon et al., 2022). We refer the reader to Howard et al. (2021); Waudby-
Smith and Ramdas (2023) for more history and modern work on this topic.

Thus far, most work on confidence sequences for the mean of random
observations has been in the scalar setting (with some exceptions discussed
below). However, estimating higher-dimensional means arises naturally in a
wide variety of statistical tasks, such as multi-armed bandits (Dani et al.,
2008; Abbasi-Yadkori et al., 2011), regression (Lai and Wei, 1982), and
stochastic optimization (Feldman et al., 2017) to name a few. Advances in
time-uniform confidence regions for estimators can therefore translate to
sequentially-valid versions of such tasks.

In this work we construct CSSs by means of PAC-Bayesian analysis (introduced
below), inspired by the original work of Catoni and Giulini (2017, 2018) which
is in the fixed-time setting. We combine their insights with recent work which
provides a general framework for understanding how and when PAC-Bayesian
bounds can be made time-uniform (Chugg et al., 2023). The resulting synthesis
yields a general method for constructing CSSs under a variety of
distributional assumptions. Building off of recent univariate work by Wang and
Ramdas (2023b), we then derive two so-called _Huber-robust_ CSSs which
maintain their coverage guarantees in the presence of (a generalization of)
Huber contamination (Huber, 1964).

In slightly more detail, our contributions are as follows. We provide four
main results. Theorem 1 provides an empirical-Bernstein CSS which can be seen
as a multivariate analogue of a recent result for bounded random variables
(Waudby-Smith and Ramdas, 2023, Theorem 2). As a direct consequence, we obtain
a novel empirical-Bernstein confidence interval whose asymptotic width scales
with the (unknown) variance. Theorem 2 then provides a CSS for sub-$\psi$
random vectors (which includes sub-Gaussian, sub-exponential, and sub-Gamma),
which again generalizes known univariate results. Theorem 3 provides a time-
uniform extension of the Catoni-Giulini estimator (Catoni and Giulini, 2018),
which results in a CSS under only second moment assumptions. Finally, Theorem
4 provides a Huber-robust extension of the empirical-Bernstein CSS above. All
of our results are closed-form and easily computable.

The appendix contains several peripheral results, such as a second Catoni-
Giulini-type CSS that can be made robust (but is less explicit than our other
bounds; Appendix B) and a semi-empirical CSS for heavy-tailed distributions
(Appendix C). Of possible independent interest, we prove several results
pertaining to the von-Mises Fisher distribution in Appendix E, such as a bound
on its KL divergence (Lemma 10).

###  1.1 Related work

Given its centrality to many statistical procedures, estimating the mean of
random vectors is a well studied topic. In the fixed-time setting, the sample
mean has nearly optimal behavior for sufficiently light-tailed distributions,
but can perform poorly in heavy-tailed settings (Catoni, 2012). This led to
the development of alternatives such as median-of-means (Lugosi and Mendelson,
2019b) and threshold-based estimators (Catoni and Giulini, 2018). See Lugosi
and Mendelson (2019a) for a survey of such results.

In the sequential setting, the theory of “self-normalized processes” (de la
Peña et al., 2004, 2009) has been leveraged to give anytime-valid confidence
ellipsoids under sub-Gaussian assumptions (Abbasi-Yadkori et al., 2011). Very
recently, Whitehouse et al. (2023) extended this approach to more general
settings. Such results are stated as bounds on
$\|V_{t}^{-1/2}\boldsymbol{S}_{t}\|$, where $(\boldsymbol{S}_{t})_{t\geqslant
1}$, $\boldsymbol{S}_{t}\in\mathbb{R}^{d}$ is a process of interest (e.g.,
$S_{t}=\sum_{i\leqslant t}\boldsymbol{X}_{i}-\boldsymbol{\mu}$) and
$(V_{t})_{t\geqslant 1}$ is the associated “variance process,” where
$V_{t}\in\mathbb{R}^{d\times d}$ is positive semi-definite. Self-normalized
bounds are useful when the variance of the distributions is different in
different directions, in which case one desires ellipsoids instead of spheres.
We briefly discuss how our bounds can sometimes be adapted to such
distributions in Section 5.1, but in general our techniques do not yield such
variance-conscious bounds and are thus not directly comparable to this line of
work. Moreover, such results of Abbasi-Yadkori et al. (2011); Whitehouse et
al. (2023) only apply to sufficiently light-tailed distributions (or to
symmetric heavy-tailed distributions), and not to the more general heavy-
tailed settings we consider in Section 3.

Meanwhile, Manole and Ramdas (2023) provide anytime-bounds on the mean of
random vectors using reverse submartingales. Their bounds achieve the optimal
rate of $\sqrt{\log\log(t)/t}$ (Darling and Robbins, 1967a). Their results are
most comparable to our results on for sub-Gaussian distributions, as they do
not hold under heavy tails. We will discuss the relationship further in
Section 2. As for anytime bounds under heavy-tails, Wang and Ramdas (2023a)
provide confidence sequences which hold under a finite $k$-th moment
assumption for $k\in(1,2]$. Our additional results in Appendix B can be seen
as an extension of their ideas to the multivariate setting, though such an
approach does not lead to closed-form bounds.

Finally, let us briefly discuss robust estimation. In the fixed-time setting
there has been substantial work on the topic (see, e.g., Prasad et al. (2020a,
b); Diakonikolas et al. (2019); Hopkins et al. (2020)). In the sequential
univariate setting, Huber (1965) and Quang (1985) studied robust statistical
tests, which are intimately linked with confidence sequences. They, however,
work in a parametric setting. The nonparametric setting was recently studied
by Wang and Ramdas (2023b) who provide robust confidence sequences for heavy-
tailed distributions. Our work will employ their framework (extended to the
multivariate setting) which generalizes the widely used concept of Huber
contamination.

###  1.2 Background and approach

Let us introduce some necessary technical tools. A (forward) filtration
$\mathcal{F}\equiv(\mathcal{F}_{t})_{t\geqslant 0}$ is a sequence of
$\sigma$-fields such that $\mathcal{F}_{t}\subset\mathcal{F}_{t+1}$,
$t\geqslant 1$. In this work we always take
$\mathcal{F}_{t}=\sigma(\boldsymbol{X}_{1},\dots,\boldsymbol{X}_{t})$, the
$\sigma$-field generated by the first $t$ observations. A stochastic process
$S=(S_{t})_{t\geqslant 1}$ is _adapted_ to $(\mathcal{F}_{t})$ if $S_{t}$ is
$\mathcal{F}_{t}$ measurable for all $t$ (meaning that $S_{t}$ depends on
$\boldsymbol{X}_{1},\dots,\boldsymbol{X}_{t}$), and _predictable_ if $S_{t}$
is $\mathcal{F}_{t-1}$ measurable for all $t$ (meaning that $S_{t}$ depends
only on $\boldsymbol{X}_{1},\dots,\boldsymbol{X}_{t-1}$). A $P$-_martingale_
is an integrable stochastic process $M\equiv(M_{t})_{t\geqslant 1}$ that is
adapted to $\mathcal{F}$, such that
$\mathbb{E}_{P}[M_{t+1}|\mathcal{F}_{t}]=M_{t}$ for all $t$. If ‘=’ is
replaced with ‘$\leqslant$’, then $M$ is a $P$-_supermartingale_. We often
abbreviate nonnegative $P$-supermartingale as $P$-NSM.

Our results rely on PAC-Bayesian theory which has emerged as successful method
to bound the generalization gap of randomized predictors in statistical
learning (Shawe-Taylor and Williamson, 1997; McAllester, 1998). As the name
suggests, PAC-Bayesian bounds employ the “probably approximately correct”
(PAC) learning framework of Valiant (1984) but addressed from a Bayesian
perspective. Roughly speaking, rather than bounding the worst case risk, PAC-
Bayesian bounds being with a prior over the parameter space $\Theta$, and
provide uniform guarantees over all “posterior” distributions. Instead of
containing a measure of the complexity of $\Theta$, the bounds involve a
divergence term between the prior and posterior. We refer the unfamiliar
reader to Guedj (2019) and Alquier (2021) for an introduction to such bounds
and their applications.

For our purposes, the important feature of PAC-Bayesian bounds is the
uniformity over posterior distributions. The insight of Catoni and Giulini
(2018) is to translate this uniformity of distributions into a bound on
$\langle\boldsymbol{\vartheta},\boldsymbol{\epsilon}\rangle$ across all
$\boldsymbol{\vartheta}\in\mathbb{S}^{d-1}$, where $\boldsymbol{\epsilon}$ is
the error of the estimator. Such a bound translates immediately into a
confidence sphere. Recently, uniformity over posteriors was extended to
uniformity over posteriors _and time_ by Chugg et al. (2023).

######  Proposition 1 (Corollary of Theorem 1, Chugg et al. (2023)).

Let $Q(\theta)\equiv(Q_{t}(\theta))_{t\geqslant 1}$ be a nonnegative
supermartingale for each $\theta\in\Theta$. Consider a prior distribution
$\nu$ over $\Theta$ (chosen before seeing the data). Then, with probability at
least $1-\alpha$, we have that simultaneously for all times $t\geqslant 1$ and
posteriors $\rho\in\mathscr{M}(\Theta)$,

$\int_{\Theta}\log Q_{t}(\theta)\rho(\mathrm{d}\theta)\leqslant
D_{\operatorname{KL}}(\rho\|\nu)+\log(1/\alpha).$ (2)

Here, $\mathscr{M}(\Theta)$ is the set of distributions over $\Theta$ and
$D_{\operatorname{KL}}(\cdot\|\cdot)$ is the KL divergence, which we recall is
defined as
$D_{\operatorname{KL}}(\rho\|\nu)=\int_{\Theta}\log(\frac{\mathrm{d}\rho}{\mathrm{d}\nu})\mathrm{d}\rho$
if $\rho\ll\nu$ and $\infty$ otherwise. The uniformity over time and
posteriors refers to the form $\Pr(\forall t,\forall\rho,\eqref{eq:pacb-
master}\text{ holds})\geqslant 1-\alpha$ of the above bound, where the
quantifiers are inside the probability.

At its core, our approach involves a judicious selection of $\Theta$, $Q$, and
family of posteriors $\\{\rho\\}$ and then applying Proposition 1. We will
take $\Theta$ to be either
$\mathbb{S}^{d-1}:=\\{\boldsymbol{x}\in\mathbb{R}^{d}:\norm{\boldsymbol{x}}=1\\}$
(in light-tailed settings) or $\mathbb{R}^{d}$ (in heavy-tailed settings), and
choose a family of posteriors $\\{\rho_{\boldsymbol{\vartheta}}\\}$ indexed by
vectors $\boldsymbol{\vartheta}\in\mathbb{S}^{d-1}$. In light-tailed settings
this will be the von Mises-Fisher (vMF) distribution; in heavy-tailed settings
a multivariate Gaussian. The vMF is defined only for $d\geqslant 2$ which is
why we restrict our attention to non-scalar settings (though it’s perhaps
worth noting that our results for heavy-tailed distributions hold for all
$d\geqslant 1$). We will typically choose supermartingales $Q$ to be functions
of $\sum_{i\leqslant
t}\langle\boldsymbol{\theta},f_{i}(\boldsymbol{X}_{i})-\boldsymbol{\mu}\rangle$
for some functions $f_{i}$, with the aim of obtaining bounds on
$\sup_{\boldsymbol{\vartheta}\in\mathbb{S}^{d-1}}\langle\boldsymbol{\vartheta},\sum_{i}(f_{i}(\boldsymbol{X}_{i})-\boldsymbol{\mu})\rangle)$.
This is in turn a bound on
$\norm{\sum_{i}(f_{i}(\boldsymbol{X}_{i})-\boldsymbol{\mu})}$ which will
furnish a CSS. Throughout, the norm $\norm{\cdot}$ should be taken to be the
$\ell^{2}$-norm $\norm{\cdot}_{2}$.

###  1.3 Assumptions

The most common assumption in the literature is that the data are
independently and identically distributed (iid). This is the case in most
prior work on estimating means (cf. Catoni and Giulini (2017); Lugosi and
Mendelson (2019b); Devroye et al. (2016); Joly et al. (2017)). We label this
Assumption 1:

######  Assumption 1.

$\boldsymbol{X}_{1},\boldsymbol{X}_{2},\dots$ are iid with mean
$\boldsymbol{\mu}$.

In this work, we are able to weaken this condition. We enforce only that the
data stream has a constant _conditional_ mean $\boldsymbol{\mu}$. This allows
for more “adversarial” distributions; see Waudby-Smith and Ramdas (2023,
Section 4.3) for more details.

######  Assumption 2.

$\boldsymbol{X}_{1},\boldsymbol{X}_{2},\dots$ have constant conditional mean
$\mathbb{E}[\boldsymbol{X}_{t}|\mathcal{F}_{t-1}]=\boldsymbol{\mu}$.

Of course, Assumption 1 implies Assumption 2. Assumption 2 is common in work
on confidence sequences, typically because they rely on (super)martingales
which naturally enable such flexibility. While Abbasi-Yadkori et al. (2011)
and Whitehouse et al. (2023) also work under more general conditions than iid
data, we are—as far as we are aware—the first to do so in heavy-tailed
settings.

##  2 Light-tailed random vectors

Throughout this section, we will take our parameter space to be
$\Theta=\mathbb{S}^{d-1}$ (this is not an assumption on the true mean, just a
tool for applying our PAC-Bayes techniques). We will apply Proposition 1 with
the von Mises-Fisher (vMF) family of distributions (Fisher, 1953) (also known
as Langevin distributions (Watamori, 1996)). The vMF distribution has density
$\gamma(\boldsymbol{x};\boldsymbol{\vartheta},\kappa)=C_{d}(\kappa)\exp(\kappa\langle\boldsymbol{\vartheta},\boldsymbol{x}\rangle)$,
where $C_{d}(\kappa)=\kappa^{d/2-1}/[(2\pi)^{d/2}I_{d/2-1}(\kappa)]$ is the
normalization constant. Here, $I_{\ell}$ is the modified Bessel function of
the first kind of order $\ell$, and $\kappa>0$ is a scalar “concentration
parameter”. The expected value of the vMF distribution obeys

$\mathbb{E}_{\boldsymbol{\theta}\sim\gamma_{\boldsymbol{\vartheta}}}\boldsymbol{\theta}=\int_{\mathbb{S}^{d-1}}\boldsymbol{\theta}\gamma(\mathrm{d}\boldsymbol{\theta};\boldsymbol{\vartheta},\kappa)=A_{d}(\kappa)\boldsymbol{\vartheta},$
(3)

where $A_{d}(\kappa)=I_{d/2}(\kappa)/I_{d/2-1}(\kappa)$. Moreover, Lemma 10 in
Appendix E proves that

$D_{\operatorname{KL}}(\gamma(\cdot;\boldsymbol{\vartheta}_{1},\kappa)\|\gamma(\cdot;\boldsymbol{\vartheta}_{2},\kappa))\leqslant
2\kappa A_{d}(\kappa).$ (4)

We refer to Mardia et al. (2000) for more on the vMF distribution, and for a
general introduction to the field of “directional statistics,” which considers
distributions on the sphere. The vMF is only well-defined for $d\geqslant 2$
so we restrict ourselves to this setting.

###  2.1 Empirical-Bernstein bound

We begin by considering bounded random vectors. Suppose there exists some
$B\in\mathbb{R}_{+}$ such that that
$\sup_{\boldsymbol{X}\in\mathcal{X}}\norm{\boldsymbol{X}}\leqslant B$. For any
$t$, let $\bar{\boldsymbol{\mu}}_{t}:=\frac{1}{t}\sum_{i\leqslant
t}\boldsymbol{X}_{i}$ be the empirical mean at time $t$, and define the
function $\psi_{E}(\lambda):=|\log(1-\lambda)+\lambda|$ for $\lambda\in[0,1)$.
The subscript “E” refers to sub-exponential (Howard et al., 2020). The
following supermartingale is a multivariate analogue of that used by Waudby-
Smith and Ramdas (2023, Equation (13)) to construct confidence sequences,
which is based off of work by Howard et al. (2020, 2021), which in turn
generalized a lemma by Fan et al. (2015).

######  Lemma 1.

Let $(\lambda_{t})$ be a predictable sequence in $[0,1)$. For each
$\boldsymbol{\theta}\in\mathbb{S}^{d-1}$, the process
$N(\boldsymbol{\theta})\equiv(N_{t}(\boldsymbol{\theta}))_{t\geqslant 1}$ is a
$P$-NSM, where

$\normalsize N_{t}(\boldsymbol{\theta})=\prod_{i\leqslant
t}\exp\left\\{\frac{\lambda_{i}}{2B}\langle\boldsymbol{\theta},\boldsymbol{X}_{i}-\boldsymbol{\mu}\rangle-\frac{\psi_{E}(\lambda_{i})}{(2B)^{2}}\langle\boldsymbol{\theta},\boldsymbol{X}_{i}-\bar{\boldsymbol{\mu}}_{i-1}\rangle^{2}\right\\}.$

The proof may be found in Appendix D. Let us apply Proposition 1 to Lemma 1.
Taking $\gamma_{\boldsymbol{1}}$ as our prior (in which we specify $\kappa$
later), we employ (4) to obtain that with probability at least $1-\alpha$,
simultaneously for all $t\geqslant 1$ and
$\boldsymbol{\vartheta}\in\mathbb{S}^{d-1}$,

$\displaystyle\int_{\mathbb{S}^{d-1}}\sum_{i\leqslant
t}\frac{\lambda_{i}}{2B}\langle\boldsymbol{\theta},\boldsymbol{X}_{i}-\boldsymbol{\mu}\rangle\gamma_{\boldsymbol{\vartheta}}(\mathrm{d}\boldsymbol{\theta})$
$\displaystyle\leqslant\sum_{i\leqslant
t}\frac{\psi_{E}(\lambda_{i})}{(2B)^{2}}\int_{\mathbb{S}^{d-1}}\langle\boldsymbol{\theta},\boldsymbol{X}_{i}-\bar{\boldsymbol{\mu}}_{i-1}\rangle^{2}\gamma_{\boldsymbol{\vartheta}}(\mathrm{d}\boldsymbol{\theta})+2\kappa
A_{d}(\kappa)+\log\frac{1}{\alpha}.$

Using (3), the left hand side above is equal to
$\frac{A_{d}(\kappa)}{2B}\sum_{i\leqslant
t}\lambda_{i}\langle\boldsymbol{\vartheta},\boldsymbol{X}_{i}-\boldsymbol{\mu}\rangle$.
Since $\norm{\boldsymbol{\theta}}=1$, we have
$\mathbb{E}_{\boldsymbol{\theta}}\langle\boldsymbol{\theta},\boldsymbol{X}_{i}-\bar{\boldsymbol{\mu}}_{i-1}\rangle^{2}\leqslant\norm{\boldsymbol{X}_{i}-\bar{\boldsymbol{\mu}}_{i-1}}^{2}$.
Rearranging gives the following result.

######  Theorem 1.

Suppose there exists $B>0$ such that $\sup\norm{\boldsymbol{X}}\leqslant B$.
Let $(\lambda_{t})_{t\geqslant 1}$ be a predictable sequence of values in
$(0,1)$. Then, under Assumption 2, for all $\kappa\geqslant 0$ and
$\alpha\in(0,1)$, with probability at least $1-\alpha$, simultaneously for all
$t\geqslant 1$,

$\displaystyle\norm{\frac{\sum_{i\leqslant
t}\lambda_{i}\boldsymbol{X}_{i}}{\sum_{i\leqslant
t}\lambda_{i}}-\boldsymbol{\mu}}$
$\displaystyle\leqslant\frac{\frac{1}{2B}\sum_{i\leqslant
t}\psi_{E}(\lambda_{i})\norm{\boldsymbol{X}_{i}-\bar{\boldsymbol{\mu}}_{i-1}}^{2}}{A_{d}(\kappa)\sum_{i\leqslant
t}\lambda_{i}}+\frac{4B\kappa
A_{d}(\kappa)+2B\log\frac{1}{\alpha}}{A_{d}(\kappa)\sum_{i\leqslant
t}\lambda_{i}}.$ (5)

The proof is supplied in Section D for completeness. This is an example of an
empirical Bernstein bound (the constant $B$ is often assumed to be known), as
it adapts to the variance of the distribution. See Figure 1. However, as
presented, the bound is somewhat obscured by the appearance of
$A_{d}(\kappa)$. To gain some intuition, Figure 3 in Appendix E illustrates
the behavior of $A_{d}(\kappa)$ as a function of $\kappa$ for various values
of $d$. Formally, Lemma 9 proves that $\frac{2}{3}<\sqrt{d}A_{d}(\sqrt{d})<1$
for all $d\geqslant 1$ which suggests setting $\kappa=\sqrt{d}$. With this
choice we obtain the following more digestible result.

######  Corollary 1.

Suppose there exists $B>0$ such that $\sup\norm{\boldsymbol{X}}\leqslant B$.
Let $(\lambda_{t})_{t\geqslant 1}$ be a predictable sequence of values in
$(0,1)$. Then, under Assumption 2, with probability at least $1-\alpha$,
simultaneously for all $t\geqslant 1$,

$\normalsize\norm{\frac{\sum_{i}\lambda_{i}\boldsymbol{X}_{i}}{\sum_{i}\lambda_{i}}-\boldsymbol{\mu}}\leqslant\frac{\frac{\sqrt{d}}{2B}\sum_{i}\psi_{E}(\lambda_{i})\overline{Z}_{i}^{2}+4B\sqrt{d}+2B\sqrt{d}\log\frac{1}{\alpha}}{\frac{2}{3}\sum_{i\leqslant
t}\lambda_{i}},$

where
$\overline{Z}_{i}=\norm{\boldsymbol{X}_{i}-\bar{\boldsymbol{\mu}}_{i-1}}$.

Theorem 1 and Corollary 1 immediately yield fixed-time empirical-Bernstein
CIs. In Appendix F.1, we analyze the asymptotic width of the CI under iid data
when taking
$\lambda_{i}=\sqrt{\nicefrac{{c\log(1/\alpha)}}{{\widehat{\sigma}_{i-1}^{2}n}}}$,
where $\widehat{\sigma}_{n}^{2}=\frac{1}{n}\sum_{i\leqslant
n}\norm{\boldsymbol{X}_{i}-\bar{\boldsymbol{\mu}}_{i}}^{2}$ is an empirical
estimate of the variance
$\sigma^{2}=\mathbb{V}(\boldsymbol{X})=\mathbb{E}[\norm{\boldsymbol{X}-\boldsymbol{\mu}}^{2}]$.
For $B=1/2$ and an appropriately chosen $c$, we demonstrate that the width
$W_{n}$ (the RHS of (5)) obeys

$\sqrt{n}W_{n}\xrightarrow{a.s.}W_{\infty}(\alpha,d)\leqslant\zeta_{\alpha}\sigma\sqrt{d\log(1/\alpha)},$

where $\zeta_{\alpha}$ is a constant depending on (and decreasing with)
$\alpha$. For reasonable choices of $\alpha$ ($\alpha\leqslant 0.1$, say)
$\zeta_{\alpha}\leqslant 3.25$. Therefore, the width of the bound scales with
the true unknown variance. This can be seen as a generalization of known
results in $d=1$ to multivariate settings. Indeed, in the scalar setting, the
empirical-Bernstein bound of Waudby-Smith and Ramdas (2023, Theorem 2) gives
an asymptotic width of $\sigma\sqrt{2\log(2/\alpha)}$, and Maurer and Pontil
(2009) provide a bound which scales as $\sigma\sqrt{4\log(2/\alpha)}$.

In the sequential setting, setting
$\lambda_{t}\asymp\sqrt{\nicefrac{{\log(1/\alpha)}}{{t\log t}}}$ results in a
width scaling as

$\widetilde{O}\left(B\sqrt{\frac{d\log(1/\alpha)\log(t)}{t}}\right),$

where $\widetilde{O}$ hides iterated logarithm factors. Combined with the
insights in the fixed-time setting explored above, we suggest taking
$\lambda_{t}\asymp\sqrt{\nicefrac{{\log(1/\alpha)}}{{\widehat{\sigma}_{t-1}^{2}t\log
t}}}$. The extra $\log(t)$ in the denominator is required to ensure the rate
is $\sqrt{\log(t)/t}$. We emphasize that our bounds hold for all predictable
sequences $(\lambda_{t})$ in $(0,1)$, but the precise selection matters for
the asymptotic width and the empirical performance.

The law of the iterated logarithm (LIL) (Darling and Robbins, 1967b) says that
the optimal rate is $O(\sqrt{\log\log(t)/t})$. (See also Wang and Ramdas
(2023a, Proposition 12) for an explicit statement and proof of the LIL lower
bound.) While such a rate is precluded by a single application of our bounds
above, it is possible to have the width shrink at the optimal rate via a
technique known as “stitching” (Howard et al., 2021) which involves deploying
Theorem 1 over geometrically spaced epochs. Details are given in Appendix A.

Let us end this section with an open problem. It was recently shown that
univariate “betting-style” confidence sequences (Waudby-Smith and Ramdas,
2023) for bounded random variables are nearly optimal (Shekhar and Ramdas,
2023). The resulting intervals are, however, not in closed-form. In the scalar
case they can still be solved for efficiently, but the naive extension to the
multivariate case does not have this property. Can a computationally tractable
betting-style CSS for multivariate means be given?

Figure 1: CSS width of the empirical-Bernstein bound in Theorem 1. We use
$d=10$ and we draw independent random variables in each coordinate from a Beta
distribution. Note that Beta(1,1) is uniform whereas Beta(50,10) is
concentrated around $5/6$. Thus, as the variance decreases, we see that the
bound gets tighter.

###  2.2 Sub-$\psi$ bound

Here we relax the boundedness assumption and instead impose a weaker tail
assumption. We assume that there exists a function
$\psi:[0,\lambda_{\max})\to\mathbb{R}$ such that

$\sup_{\boldsymbol{v}\in\mathbb{S}^{d-1}}\mathbb{E}[e^{\lambda\langle\boldsymbol{v},\boldsymbol{X}_{t}-\boldsymbol{\mu}\rangle}|\mathcal{F}_{t-1}]\leqslant
e^{\psi(\lambda)}$ (6)

for all $\lambda\in[0,\lambda_{\max})$. For instance, if
$\psi(\lambda)=\frac{\lambda^{2}}{2(1-b\lambda)}$ then $P$ is $b$-sub-Gamma,
and if $\psi(\lambda)=\lambda^{2}\sigma^{2}/2$ then it is
$(\sigma^{2},\lambda_{\max})$-sub-exponential. Condition (6) was employed by
Manole and Ramdas (2023), and it is reminiscent of the “sub-$\psi$ process” of
Howard et al. (2020). Accordingly, if $P$ obeys such a condition, we will say
it is sub-$\psi$.

As before, let $\gamma_{\boldsymbol{\vartheta}}$ be the density of the vMF
distribution, and let $(\lambda_{t})_{t\geqslant 1}$ be a deterministic
sequence of weights in $[0,\lambda_{\max})$. Due to (6), the process
$(M_{t}(\boldsymbol{\theta}))_{t\geqslant 1}$ defined as
$M_{t}(\boldsymbol{\theta})=\prod_{i=1}^{t}\exp\\{\lambda_{i}\langle\boldsymbol{\theta},\boldsymbol{X}_{i}-\boldsymbol{\mu}\rangle-\psi(\lambda_{i})\\}$,
is a $P$-NSM. Applying Proposition 1 with prior $\nu=\gamma_{\boldsymbol{1}}$,
we obtain that with probability at least $1-\alpha$, for all times $t\geqslant
1$ and all $\boldsymbol{\vartheta}\in\mathbb{S}^{d-1}$,

$\sum_{i\leqslant
t}\int_{\mathbb{S}^{d-1}}\lambda_{i}\langle\boldsymbol{\theta},\boldsymbol{X}_{i}-\boldsymbol{\mu}\rangle\gamma_{\boldsymbol{\vartheta}}(\mathrm{d}\boldsymbol{\theta})\leqslant\sum_{i\leqslant
t}\psi(\lambda_{i})+2\kappa A_{d}(\kappa)+\log(\frac{1}{\alpha}).$

As before, we have
$\mathbb{E}_{\boldsymbol{\theta}\sim\gamma_{\boldsymbol{\vartheta}}}\langle\boldsymbol{\theta},\boldsymbol{X}_{i}-\boldsymbol{\mu}\rangle=A_{d}(\kappa)\langle\boldsymbol{\vartheta},\boldsymbol{X}_{i}-\boldsymbol{\mu}\rangle$.
This yields a result similar to Theorem 1 involving $A_{k}(d)$. Instead of
stating it, let us proceed directly its more concrete cousin. Taking
$\kappa=\sqrt{d}$ as was done in Section 2.1 and rearranging the above display
gives the following.

######  Theorem 2.

Let $\boldsymbol{X}_{1},\boldsymbol{X}_{2},\dots$ be sub-$\psi$. Let
$(\lambda_{t})_{t\geqslant 1}$ be a deterministic sequence in
$[0,\lambda_{\max})$. Then, under Assumption 2, for all $\alpha\in(0,1)$, with
probability at least $1-\alpha$, for all $t\geqslant 1$,

$\norm{\frac{\sum_{i\leqslant
t}\lambda_{i}\boldsymbol{X}_{i}}{\sum_{i\leqslant
t}\lambda_{i}}-\boldsymbol{\mu}}\leqslant\frac{\sqrt{d}\sum_{i\leqslant
t}\psi(\lambda_{i})+2\sqrt{d}+\sqrt{d}\log\frac{1}{\alpha}}{\frac{2}{3}\sum_{i\leqslant
t}\lambda_{i}}.$

A particularly relevant case is when $P$ is conditionally $\sigma$-sub-
Gaussian, in which case $\lambda_{\max}=\infty$ and
$\psi(\lambda)=\lambda^{2}\sigma^{2}/2$. In this case, if we take
$\lambda_{t}\asymp\sqrt{\nicefrac{{\log(1/\alpha)}}{{\sigma^{2}t\log t}}}$,
then the RHS above scales as
$\widetilde{O}\big{(}\sigma\sqrt{\nicefrac{{d\log(1/\alpha)\log(t)}}{{t}}}\big{)}$.
Again, similarly to Section 2.1, the optimal rate of $\sqrt{\log\log(t)/t}$
can be obtained with stitching in which case we recover a result with the same
rate (and dependence on $d$ and $\sigma$) as Manole and Ramdas (2023,
Corollary 23), but ours holds under slightly more general distributional
assumptions.

##  3 Heavy-tailed random vectors

Let us now assume that the random vectors $\boldsymbol{X}$ are less well-
behaved. We assume only that there exists some $p\geqslant 2$ such that

$\mathbb{E}_{P}[\norm{\boldsymbol{X}}^{p}|\mathcal{F}_{t-1}]\leqslant
v<\infty,\quad\forall t\geqslant 1.$ (7)

Note that this condition implies that the norm of the mean $\boldsymbol{\mu}$
is bounded by $v^{1/p}$:
$\mathbb{E}[\norm{\boldsymbol{X}}|\mathcal{F}_{t-1}]\leqslant\mathbb{E}[\norm{\boldsymbol{X}}^{p}|\mathcal{F}_{t-1}]^{1/p}\leqslant
v^{1/p}$ and, similarly, the second moment by $v^{2/p}$.

This section provides a time-uniform generalization of the bound given by
Catoni and Giulini (2018). We note that a time-uniform version of a second
Catoni-Giulini bound (Catoni and Giulini, 2017) is possible, but the resulting
CSS is not closed-form. This makes it less practical, so we move the
discussion of this result to Appendix B, where it lies in wait for the
mathematically curious. Throughout this section, we will take our parameter
space to be $\Theta=\mathbb{R}^{d}$ and place a family of Gaussian
distributions
$\\{\phi_{\boldsymbol{\vartheta}}:\boldsymbol{\vartheta}\in\mathbb{S}^{d-1}\\}$
over $\Theta$, where $\phi_{\boldsymbol{\vartheta}}$ is the density of the
normal distribution $N(\boldsymbol{\vartheta},\beta^{-1}\boldsymbol{I}_{d})$
and $\beta>0$. When applying Proposition 1, we will set
$\nu=\phi_{\boldsymbol{0}}$, in which case
$D_{\operatorname{KL}}(\rho_{\boldsymbol{\vartheta}}\|\nu)=\beta/2$. (Note
this formula follows from the standard KL divergence between multivariate
Gaussians.)

To motivate the bound, observe that the theorems given thus far roughly have
the form of a “point estimate + uncertainty.” In particular, they rely on the
estimators
$\widehat{\boldsymbol{\mu}}_{t}=\sum_{i}\lambda_{i}\boldsymbol{X}_{i}/\sum_{i}\lambda_{i}$.
However, it is known that such simple linear combinations of samples like the
(possibly weighted) empirical average have sub-optimal finite-sample
performance under heavy tails (like distributions with only two moments)
(Catoni, 2012). Other estimators are therefore required.

To begin, for a sequence of positive scalars $(\lambda_{t})$, define the
threshold function

$\mathsf{th}_{t}(\boldsymbol{x})=\frac{\lambda_{t}\norm{\boldsymbol{x}}\wedge
1}{\lambda_{t}\norm{\boldsymbol{x}}}\boldsymbol{x}.$ (8)

Our estimator will be $\widehat{\boldsymbol{\mu}}_{t}=\sum_{i\leqslant
t}\lambda_{i}\mathsf{th}_{i}(\boldsymbol{X}_{i})/\sum_{i\leqslant
t}\lambda_{i}$. The analysis proceeds by separating the quantity
$\langle\boldsymbol{\vartheta},\mathsf{th}_{t}(\boldsymbol{X}_{t})-\boldsymbol{\mu}\rangle$
into two terms:

$\normalsize\langle\boldsymbol{\vartheta},\mathsf{th}_{t}(\boldsymbol{X}_{t})-\boldsymbol{\mu}\rangle=\underbrace{\langle\boldsymbol{\vartheta},\mathsf{th}_{t}(\boldsymbol{X}_{t})-\boldsymbol{\mu}_{t}\rangle}_{(i)}+\underbrace{\langle\boldsymbol{\vartheta},\boldsymbol{\mu}_{t}-\boldsymbol{\mu}\rangle}_{(ii)},$
(9)

where
$\boldsymbol{\mu}_{t}=\mathbb{E}[\mathsf{th}_{t}(\boldsymbol{X})|\mathcal{F}_{t-1}]$.
The second term is dealt with by the following technical lemma, which is
proved in Appendix D.

######  Lemma 2.

Term (ii) of (9) obeys:

$\sup_{\boldsymbol{\vartheta}\in\mathbb{S}^{d-1}}\langle\boldsymbol{\vartheta},\boldsymbol{\mu}_{t}-\boldsymbol{\mu}\rangle=\norm{\boldsymbol{\mu}_{t}-\boldsymbol{\mu}}\leqslant\lambda_{t}v^{2/p}.$
(10)

To handle term (i) of (9), we appeal to a PAC-Bayesian theorem. The following,
due to Chugg et al. (2023, Corollary 15), is a time-uniform extension of the
bound by Catoni (2004, Equation (5.2.1)). It is based on applying Proposition
1 to the supermartingale defined by
$U_{t}(\boldsymbol{\theta})=\prod_{i\leqslant
t}\exp\\{\lambda_{i}f_{i}(\boldsymbol{\theta})-\log\mathbb{E}\exp(\lambda_{i}f_{i}(\boldsymbol{\theta}))\\}$.

######  Lemma 3.

Let $\\{f_{t}:\mathcal{X}\times\Theta\to\mathbb{R}\\}$ be a sequence of
measurable functions. Fix a prior $\nu$ over $\Theta$. Then, with probability
$1-\alpha$ over $P$, for all $t\geqslant 1$ and all distributions $\rho$ over
$\Theta$,

$\displaystyle\sum_{i\leqslant
t}\int_{\mathbb{R}^{d}}f_{i}(\boldsymbol{X}_{i},\boldsymbol{\theta})\rho(\mathrm{d}\boldsymbol{\theta})\leqslant\sum_{i\leqslant
t}\int_{\mathbb{R}^{d}}\log\mathbb{E}_{i-1}e^{f_{i}(\boldsymbol{X},\boldsymbol{\theta})}\rho(\mathrm{d}\boldsymbol{\theta})+D_{\operatorname{KL}}(\rho\|\nu)+\log\frac{1}{\alpha}.$

We apply Lemma 3 with the functions
$f_{i}(\boldsymbol{X}_{i},\boldsymbol{\theta})=\langle\boldsymbol{\theta},\mathsf{th}_{i}(\boldsymbol{X}_{i})-\boldsymbol{\mu}_{t}\rangle$.
Keeping in mind that
$\mathbb{E}_{\boldsymbol{\theta}\sim\phi_{\boldsymbol{\vartheta}}}\langle\boldsymbol{\theta},\mathsf{th}_{i}(\boldsymbol{X}_{i})-\boldsymbol{\mu}_{i}\rangle=\langle\boldsymbol{\vartheta},\mathsf{th}_{i}(\boldsymbol{X}_{i})-\boldsymbol{\mu}_{i}\rangle$,
we obtain that with probability $1-\alpha$, simultaneously for all $t\geqslant
1$,

$\displaystyle\sup_{\boldsymbol{\vartheta}\in\mathbb{S}^{d-1}}\sum_{i\leqslant
t}\lambda_{i}\langle\boldsymbol{\vartheta},\mathsf{th}_{i}(\boldsymbol{X}_{i})-\boldsymbol{\mu}_{i}\rangle\leqslant\sum_{i\leqslant
t}\operatorname*{\mathbb{E}}_{\boldsymbol{\theta}\sim\phi_{\boldsymbol{\vartheta}}}\log\mathbb{E}\left\\{\mathrm{e}^{\lambda_{i}\left\langle\boldsymbol{\theta},\mathsf{th}_{i}(\boldsymbol{X}_{i})-\boldsymbol{\mu}_{i}\right\rangle}|\mathcal{F}_{i-1}\right\\}+\frac{\beta}{2}+\log\frac{1}{\alpha}.$
(11)

Our second technical lemma, following Catoni and Giulini (2018), helps bound
the right hand side of the above. The proof is given in Appendix D.

######  Lemma 4.

For all $t\geqslant 1$,

$\displaystyle\operatorname*{\mathbb{E}}_{\boldsymbol{\theta}\sim\phi_{\boldsymbol{\vartheta}}}\log\operatorname*{\mathbb{E}}_{\boldsymbol{X}\sim
P}\left\\{\mathrm{e}^{\left\langle\boldsymbol{\theta},\mathsf{th}_{t}(\boldsymbol{X})-\boldsymbol{\mu}_{t}\right\rangle}|\mathcal{F}_{t-1}\right\\}\leqslant\frac{1}{4}v^{2/p}\lambda_{t}^{2}e^{2/\beta+2}.$

Applying Lemma 4 to (11) furnishes the main result. The arithmetic is provided
in Appendix D.2.

######  Theorem 3.

Let $\boldsymbol{X}_{1},\boldsymbol{X}_{2},\dots\sim P$ with mean
$\boldsymbol{\mu}$ obey $\mathbb{E}[\norm{\boldsymbol{X}}^{p}]\leqslant
v<\infty$ for some $p\geqslant 2$. Fix any $\beta>0$. Then, under Assumption
2, with probability at least $1-\alpha$, for all $t\geqslant 1$,

$\displaystyle\left\|\frac{\sum_{i\leqslant
t}\lambda_{i}\mathsf{th}_{i}(\boldsymbol{X}_{i})}{\sum_{i\leqslant
t}^{t}\lambda_{i}}-\boldsymbol{\mu}\right\|\leqslant\frac{v^{2/p}\left(2e^{\frac{2}{\beta}+2}+1\right)\sum_{i\leqslant
t}\lambda_{i}^{2}+\frac{\beta}{2}+\log(\frac{1}{\alpha})}{\sum_{i\leqslant
t}\lambda_{i}}.$

Following a similar strategy to Section C, consider setting
$\lambda_{t}\asymp\sqrt{\frac{\log(1/\alpha)}{v^{2/p}t\log t}}$, and
$\beta=1$.222For such a choice of $\lambda$, the optimal value of $\beta$ can
be unearthed with a bit of calculus, but the final expression includes the
difficult-to-interpret Lambert $W$ function. The improvement is minor, so we
do not pursue this further. In this case, the RHS of the above display is then
bounded as

$\widetilde{O}\left(\sqrt{\frac{v^{2/p}\log(1/\alpha)\log t}{t}}\right).$ (12)

Figure 2: Comparison of the width of the CSS in Theorem 3 and the median-of-
means estimator (made anytime-valid with a union bound). We use $d=20$, mean
0, and the identity covariance matrix.

In the fixed time setting, ideal estimators are sub-Gaussian (Lugosi and
Mendelson, 2019a, Equation (3.1)), in the sense that they have rate

$O\left(\sqrt{\frac{\lambda_{\max}\log(1/\alpha)}{n}}+\sqrt{\frac{\Tr(\Sigma)}{n}}\right),$
(13)

where $\lambda_{\max}=\lambda_{\max}(\Sigma)$ is the maximum eigenvalue of the
covariance matrix $\Sigma$. For us, setting
$\lambda=\lambda_{i}=\sqrt{\frac{\log(1/\alpha)}{v^{2/p}n}}$ for all $i$,
yields a slightly larger rate $O\left(\sqrt{v^{2/p}\log(1/\alpha)/n}\right)$.
Indeed,
$\lambda_{\max}\leqslant\mathbb{V}(\boldsymbol{X}|\mathcal{F}_{t-1})=\Tr(\Sigma)\leqslant
v^{2/p}$. This is no surprise: it is well-known that the Catoni-Giulini
estimator is not sub-Gaussian, hence the development of multivariate median-
of-means estimators (Lugosi and Mendelson, 2019b). An interesting open problem
is to develop a CSS for heavy-tailed distributions that precisely attains the
ideal sub-Gaussian performance (13). We compare the radius empirically to the
fixed-time bound (turned into a CSS using a union bound with
$\alpha/(t+t^{2})$ at time $t$)333See e.g. Wang and Ramdas (2023a, Remark 1).
of the median-of-means estimator (Lugosi and Mendelson, 2019a, Proposition 1),
which bears the smallest explicit constant we find in the literature, in
Figure 2.

##  4 Robust confidence spheres

Here we discuss how we can modify our empirical Bernstein CSS of Section 2.1
to ensure it is robust to Huber contamination (Tukey, 1960; Huber, 1973). To
formalize robustness, recall that the _total variation (TV) distance_ is
defined as
$D_{\operatorname{TV}}(P,Q)=\sup_{A\in\mathbb{B}(\mathbb{R}^{d})}|P(A)-Q(A)|$,
where $\mathbb{B}(\mathbb{R}^{d})$ is the set of Borel sets in
$\mathbb{R}^{d}$. The TV distance takes values in $[0,1]$ and equals 0 iff
$P=Q$ almost everywhere. A particularly useful equivalent definition of the TV
distance is: For any $a<b$,

$\normalsize
D_{\operatorname{TV}}(P,Q)=\frac{1}{b-a}\sup_{f:\mathbb{R}^{d}\to[a,b]}\left|\mathbb{E}_{P}f(\boldsymbol{X})-\mathbb{E}_{Q}f(\boldsymbol{X})\right|.$
(14)

Given a distribution $P$, following Wang and Ramdas (2023b) we define the TV
$\epsilon$-ball around $P$ as

$\mathscr{B}_{\operatorname{TV}}(P,\epsilon):=\left\\{Q\in\mathscr{M}(\mathbb{R}^{d}):D_{\operatorname{TV}}(P,Q)\leqslant\epsilon\right\\}.$
(15)

Let $\boldsymbol{\mu}_{P}$ denote the mean of distribution $P$. For
simplicity, throughout this section we will work under Assumption 1, i.e., we
assume iid observations. We would like to develop CSSs for
$\boldsymbol{\mu}_{P}$ that are robust to data generated from any
$Q\in\mathscr{B}_{\operatorname{TV}}(P,\epsilon)$, for some predetermined
$\epsilon>0$. Formally, let us say that a sequence
$(\boldsymbol{C}_{t})_{t\geqslant 1}$ is a $(1-\alpha)$-CSS for
$\boldsymbol{\mu}_{P}$ if

$\forall Q\in\mathscr{B}_{\operatorname{TV}}(P,\epsilon),\;Q(\forall
t\geqslant 1:\boldsymbol{\mu}_{P}\in\boldsymbol{C}_{t})\geqslant 1-\alpha.$
(16)

We note that the contamination model used by Huber (cf. Huber (1964, 1973))
was the “$\epsilon$-contamination model,” in which data was generated from $P$
with probability $(1-\epsilon)P$ and from some other arbitrary distribution
$Q$ with probability $\epsilon$. As noted by Wang and Ramdas (2023b), such
mixture distributions can be shown to be a subset of
$\mathscr{B}_{\operatorname{TV}}(P,\epsilon)$.

Fix a distribution $P$ with observations $\boldsymbol{X}$ such that
$\norm{\boldsymbol{X}}\leqslant B$. We begin by modifying the process
$(N_{t}(\boldsymbol{\theta}))$ of Lemma 1 to be a $Q$-NSM for all
$Q\in\mathscr{B}_{\operatorname{TV}}(P,\epsilon)$. Let
$Z_{t}(\boldsymbol{\theta})=\langle\boldsymbol{\theta},\boldsymbol{X}_{t}-\boldsymbol{\mu}\rangle$.
For $\lambda_{t}\in[0,0.8]$, we have

$\normalsize
F_{t}(\boldsymbol{\theta}):=\frac{\lambda_{t}}{2B}Z_{t}(\boldsymbol{\theta})-\frac{\psi_{E}(\lambda_{t})}{(2B)^{2}}\langle\boldsymbol{\theta},\boldsymbol{X}_{t}-\bar{\boldsymbol{\mu}}_{t-1}\rangle^{2}\in[-1,2],$

since the absolute value of both terms is at most 1 and the second term is
nonnegative. Hence $\exp F_{t}(\boldsymbol{\theta})\in[e^{-1},e^{2}]$ and we
can apply definition (14) with the function $f(\boldsymbol{X}_{t})=\exp
F_{t}(\boldsymbol{\theta})$ to obtain the following result. Details are given
in Appendix D.3.

######  Lemma 5.

Let $(\lambda_{t})$ be a sequence of scalars in $(0,0.8]$ and let
$\epsilon>0$. For all $\boldsymbol{\theta}\in\mathbb{S}^{d-1}$, the process
$N^{R}(\boldsymbol{\theta})\equiv(N_{t}^{R}(\boldsymbol{\theta}))_{t\geqslant
1}$ given by

$N_{t}^{R}(\boldsymbol{\theta})=\prod_{i\leqslant
t}\frac{\exp(F_{i}(\boldsymbol{\theta}))}{1+e^{2}\epsilon},$

is a $Q$-NSM for all $Q\in\mathscr{B}_{\operatorname{TV}}(P,\epsilon)$.

Here we use the superscript $R$ to signal robustness. Since we are dealing
with bounded observations, we return to the setting and tools of Section 2 and
employ the vMF distribution alongside Proposition 1. Taking as our prior
$\gamma_{\boldsymbol{1}}$ gives the resulting robust CSS.

######  Theorem 4.

Let $P$ be a distribution over
$\\{\boldsymbol{x}\in\mathbb{R}^{d}:\norm{\boldsymbol{x}}\leqslant B\\}$ with
mean $\boldsymbol{\mu}_{P}$. Let $(\lambda_{t})$ be a sequence in $(0,0.8]$.
For $\epsilon>0$, let $\boldsymbol{X}_{1},\boldsymbol{X}_{2},\dots\sim Q$ for
$Q\in\mathscr{B}_{\operatorname{TV}}(P,\epsilon)$. Then for any
$\alpha\in(0,1)$, with probability $1-\alpha$, for all $t\geqslant 1$,

$\displaystyle\norm{\frac{\sum_{i}\lambda_{i}\boldsymbol{X}_{i}}{\sum_{i}\lambda_{i}}-\boldsymbol{\mu}}$
$\displaystyle\leqslant\frac{\frac{\sqrt{d}}{2B}\sum_{i}\psi_{E}(\lambda_{i})\overline{Z}_{i}^{2}+6B\sqrt{d}\log\frac{1}{\alpha}}{\frac{2}{3}\sum_{i\leqslant
t}\lambda_{i}}+\frac{2B\sqrt{d}t\log(1+e^{2}\epsilon)}{\frac{2}{3}\sum_{i\leqslant
t}\lambda_{i}},$

where we recall that $\psi_{E}(\lambda_{i})=|\log(1-\lambda_{i})+\lambda_{i}|$
and we’ve set
$\overline{Z}_{i}=\norm{\boldsymbol{X}_{i}-\bar{\boldsymbol{\mu}}_{i-1}}$.

The proof is similar to that of Theorem 1 but is provided in Appendix D.3 for
completeness. The factor of $t$ in the numerator of the bound may seem
alarming, but it is a necessary price to pay for robustness. Indeed, it is
well-known that robust confidence sequences and intervals, and thus CSSs,
cannot have widths that shrink to zero (Wang and Ramdas, 2023b, Lemma 1). In
fact, error scaling as a function of $\epsilon$ is impossible to avoid (Chen
et al., 2018; Lugosi and Mendelson, 2021). Consequently, we search for
alternative strategies of how to set $\lambda_{t}$ than those of previous
sections.

Taking $\lambda$ to be constant results in a CSS with asymptotic width
$O(\sqrt{d}\sigma^{2}/B+\sqrt{B}\epsilon)$ where
$\sigma^{2}=\mathbb{V}(\boldsymbol{X})=\mathbb{E}{\norm{\boldsymbol{X}-\boldsymbol{\mu}}}^{2}$.
In the fixed-time setting, we can optimize $\lambda$ for time $t=n$. This
results in a bound with width $O(\sigma\sqrt{d\epsilon})$ as $n\to\infty$. In
the sequential setting, we might set $\lambda_{t}\asymp
1/\widehat{\sigma}_{t-1}$ where
$\widehat{\sigma}_{t}^{2}=\frac{1}{t}\sum_{i}\overline{Z}_{i}^{2}$ is the
empirical variance proxy. If $\lambda_{t}\leqslant 0.68$, then the width of
the CSS converges almost surely to $O(\sigma\sqrt{B}(\sqrt{d}+\epsilon))$ as
$t\to\infty$, which matches the lower bound (in both $d$ and $\epsilon$) due
to Chen et al. (2018, Theorem 2.2). A longer discussion can be found in
Appendix F.2.

##  5 Discussion

###  5.1 Adaptivity to anisotropy via ellipsoids

We have thus far made isotropic assumptions on the scale of the random vectors
(2-norm bounds, $\sigma^{2}$-sub-Gaussianity, and second moments of the
2-norm). However, it is not hard to see that by applying an appropriate
“whitening transformation” to the random vectors, we can obtain confidence
_ellipsoid_ sequences that adapt to the possible anisotropy of distributions.
For example, if $\boldsymbol{X}_{t}$’s are conditionally $\Sigma$-sub-Gaussian
(i.e.
$\mathbb{E}[e^{\langle\boldsymbol{u},\boldsymbol{X}_{t}-\boldsymbol{\mu}\rangle}|\mathcal{F}_{t-1}]\leqslant
e^{\langle\boldsymbol{u},\Sigma\boldsymbol{u}\rangle/2}$ for all
$\boldsymbol{u}\in\mathbb{R}^{d}$), applying Theorem 2 to the 1-sub-Gaussian
vectors $\Sigma^{-1/2}\boldsymbol{X}$ gives rise to the following
$(1-\alpha)$-confidence ellipsoid sequence for $\boldsymbol{\mu}$:

$\norm{\frac{\sum_{i\leqslant
t}\lambda_{i}\boldsymbol{X}_{i}}{\sum_{i\leqslant
t}\lambda_{i}}-\boldsymbol{\mu}}_{\Sigma}\leqslant\frac{\frac{\sqrt{d}}{2}\sum_{i\leqslant
t}\lambda_{i}^{2}+2+\sqrt{d}\log(1/\alpha)}{0.8\sum_{i\leqslant
t}\lambda_{i}},$

where $\norm{\boldsymbol{x}}_{\Sigma}=\norm{\Sigma^{-1/2}\boldsymbol{x}}$
denotes the Mahalanobis norm. Similarly, Theorem 3 implies the following
$(1-\alpha)$-confidence ellipsoid sequence

$\displaystyle\left\|\frac{\sum_{i\leqslant
t}\lambda_{i}\mathsf{th}^{\mathsf{MH}}_{i}(\boldsymbol{X}_{i})}{\sum_{i\leqslant
t}^{t}\lambda_{i}}-\boldsymbol{\mu}\right\|_{\Sigma}$
$\displaystyle\leqslant\frac{\sqrt{d}\left(2e^{\frac{2}{\beta}+2}+1\right)\sum_{i\leqslant
t}\lambda_{i}^{2}+\frac{\beta}{2}+\log(\frac{1}{\alpha})}{\sum_{i\leqslant
t}\lambda_{i}},$

where
$\mathsf{th}^{\mathsf{MH}}_{t}(\boldsymbol{x})=\frac{\lambda_{t}\norm{\boldsymbol{x}}_{\Sigma}\wedge
1}{\lambda_{t}\norm{\boldsymbol{x}}_{\Sigma}}\boldsymbol{x}$, under the
assumption that
$\mathbb{E}[\boldsymbol{X}\boldsymbol{X}^{\mathsf{T}}|\mathcal{F}_{t-1}]\preceq\Sigma$,
where $\preceq$ represents the Loewner ordering ($A\preceq B$ means that $B-A$
is a positive semidefinite matrix). The $\sqrt{d}$ term above, which matches
the $\sqrt{d}$ in the light-tailed bounds and _seems to be worse_ than the
“dimension-free” Theorem 3, shall be discussed immediately next.

###  5.2 Dependence on $d$

We now turn to the question of why some of our CSSs have explicit $\sqrt{d}$
terms (Theorems 1 and 2) whereas others appear to be “dimension-free”
(Theorems 3 and 5). To compare them, we identify two different types of
factors on the scales (or scatteredness) of distributions: (1) Trace-type
factors which sum up the scatteredness along all $d$ dimensions, including the
$\mathbb{E}\norm{\boldsymbol{X}}^{2}\leqslant v$ assumption for Theorem 3; (2)
$\lambda_{\max}$-type factors which measure the distribution along its _most_
scattered direction, including the $\norm{\boldsymbol{X}}\leqslant B$
assumption for Theorem 1 and the $\sigma^{2}$-sub-Gaussian assumption for
Theorem 2. (Indeed, for any positive definite matrix $\Sigma$, a distribution
is $\Sigma$-sub-Gaussian if and only if it is $\lambda_{\max}(\Sigma)$-sub-
Gaussian.)

Therefore, when the distribution is roughly isotropic, one should see the
trace-type $v$ as bearing an implicit $d$ within, leading to all our spheres’
radii scaling as $\sqrt{d}$. Recall the “ideal” sub-Gaussian radius of
dependence $\sqrt{\Tr(\Sigma)}$ we cite in (13). When the distribution is
moderately anisotropic, trace-type factors will be tighter than
$\lambda_{\max}$-type factors multiplied by $\sqrt{d}$, but the gap again
disappears if one follows the recipe we recommend in Section 5.1 to obtain
ellipsoids that all adaptively scale as
$\norm{\widehat{\boldsymbol{\mu}}_{t}-\boldsymbol{\mu}}_{\Sigma}\lesssim\sqrt{d}$.
For high-dimensional cases with sparsity where the distribution almost
entirely lies within an $o(d)$-dimensional affine subspace of
$\mathbb{R}^{d}$, Theorem 3 is preferred as its $d$-free trace-type assumption
adapts to the “effective” dimensionality.

##  6 Summary

We have studied how time-uniform confidence sphere sequences can be obtained
via PAC-Bayesian inequalities. We extended the initial work of Catoni and
Giulini (2017, 2018) to (i) be time-uniform, and (ii) to handle a variety of
distributional assumptions. Our results include an empirical-Bernstein CSS for
bounded random vectors, a corollary of which is a novel empirical-Bernstein
confidence interval for multivariate means, and a method for making this
result robust to Huber noise. We also provide CSSs under more general
sub-$\psi$ tail assumptions (which includes sub-Gaussianity) and, finally,
under a heavy-tailed (second moment only) assumption.

This work suggests several intriguing research directions. First, as discussed
at the end of Section 2.1, given their near-optimality in the univariate case
(Shekhar and Ramdas, 2023), can one develop a computationally tractable
betting-style CSS? Second, while we provide a semi-empirical CSS for heavy-
tailed distributions in Appendix C, can a fully empirical CSS be given? Third,
can a CSS with sub-Gaussian rates be given in the heavy-tailed setting?
Finally, Appendix B gives a robust CSS under heavy-tails but the bound is not
closed-form. Can an explicit robust CSS be constructed in this setting?

#### Acknowledgements.

BC and AR acknowledge support from NSF grants IIS-2229881 and DMS-2310718. BC
was supported in part by the NSERC PGS D program, grant no. 567944.

## References

  * Abbasi-Yadkori et al. (2011) Yasin Abbasi-Yadkori, Dávid Pál, and Csaba Szepesvári.  Improved algorithms for linear stochastic bandits.  _Advances in Neural Information Processing Systems_ , 24, 2011.
  * Adobe (2023) Adobe.  Adobe Experimentation Panel.  https://experienceleague.adobe.com/docs/analytics-platform/using/cja-workspace/panels/experimentation.html?lang=en, 2023\.  [Online; accessed 13-October-2023].
  * Alquier (2021) Pierre Alquier.  User-friendly introduction to PAC-Bayes bounds.  _arXiv preprint arXiv:2110.11216_ , 2021.
  * Catoni (2004) Olivier Catoni.  _Statistical learning theory and stochastic optimization: Ecole d’Eté de Probabilités de Saint-Flour, XXXI-2001_ , volume 1851.  Springer Science & Business Media, 2004.
  * Catoni (2012) Olivier Catoni.  Challenging the empirical mean and empirical variance: a deviation study.  In _Annales de l’IHP Probabilités et statistiques_ , volume 48, pages 1148–1185, 2012.
  * Catoni and Giulini (2017) Olivier Catoni and Ilaria Giulini.  Dimension-free PAC-Bayesian bounds for matrices, vectors, and linear least squares regression.  _arXiv preprint arXiv:1712.02747_ , 2017.
  * Catoni and Giulini (2018) Olivier Catoni and Ilaria Giulini.  Dimension-free PAC-Bayesian bounds for the estimation of the mean of a random vector.  _(Almost) 50 Shades of Bayesian Learning: PAC-Bayesian Trends and Insights. NeurIPS Workshop._ , 2018.
  * Chen et al. (2018) Mengjie Chen, Chao Gao, and Zhao Ren.  Robust covariance and scatter matrix estimation under Huber’s contamination model.  _The Annals of Statistics_ , 46(5):1932–1960, 2018.
  * Chugg et al. (2023) Ben Chugg, Hongjian Wang, and Aaditya Ramdas.  A unified recipe for deriving (time-uniform) PAC-Bayes bounds.  _arXiv preprint arXiv:2202.01250_ , 2023.
  * Dani et al. (2008) Varsha Dani, Thomas P Hayes, and Sham M Kakade.  Stochastic linear optimization under bandit feedback.  _21st Annual Conference on Learning Theory_ , 2008.
  * Darling and Robbins (1967a) Donald A Darling and Herbert Robbins.  Confidence sequences for mean, variance, and median.  _Proceedings of the National Academy of Sciences_ , 58(1):66–68, 1967a.
  * Darling and Robbins (1967b) Donald A Darling and Herbert Robbins.  Iterated logarithm inequalities.  _Proceedings of the National Academy of Sciences_ , 57(5):1188–1192, 1967b.
  * de la Peña et al. (2004) Victor H de la Peña, Michael J Klass, and Tze Leung Lai.  Self-normalized processes: exponential inequalities, moment bounds and iterated logarithm laws.  _Annals of Probability_ , 2004.
  * de la Peña et al. (2009) Victor H de la Peña, Tze Leung Lai, and Qi-Man Shao.  _Self-normalized processes: Limit theory and Statistical Applications_.  Springer, 2009.
  * Delyon (2009) Bernard Delyon.  Exponential inequalities for sums of weakly dependent variables.  _Electronic Journal of Probability_ , 2009.
  * Devroye et al. (2016) Luc Devroye, Matthieu Lerasle, Gabor Lugosi, and Roberto I Oliveira.  Sub-Gaussian mean estimators.  _The Annals of Statistics_ , 2016.
  * Diakonikolas et al. (2019) Ilias Diakonikolas, Gautam Kamath, Daniel Kane, Jerry Li, Ankur Moitra, and Alistair Stewart.  Robust estimators in high-dimensions without the computational intractability.  _SIAM Journal on Computing_ , 48(2):742–864, 2019\.
  * Evidently (2023) Evidently.  How Evidently calculates results.  https://docs.aws.amazon.com/AmazonCloudWatch/latest/monitoring/CloudWatch-Evidently-calculate-results.html, 2023\.  [Online; accessed 13-October-2023].
  * Fan et al. (2015) Xiequan Fan, Ion Grama, and Quansheng Liu.  Exponential inequalities for martingales with applications.  _Electronic Journal of Probability_ , 20:1–22, 2015.
  * Feldman et al. (2017) Vitaly Feldman, Cristobal Guzman, and Santosh Vempala.  Statistical query algorithms for mean vector estimation and stochastic convex optimization.  In _Proceedings of the Twenty-Eighth Annual ACM-SIAM Symposium on Discrete Algorithms_ , pages 1265–1277. SIAM, 2017.
  * Fisher (1953) Ronald Aylmer Fisher.  Dispersion on a sphere.  _Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences_ , 217(1130):295–305, 1953\.
  * Guedj (2019) Benjamin Guedj.  A primer on PAC-Bayesian learning.  In _Proceedings of the second congress of the French Mathematical Society_ , volume 33, 2019.
  * Haddouche and Guedj (2023) Maxime Haddouche and Benjamin Guedj.  PAC-Bayes generalisation bounds for heavy-tailed losses through supermartingales.  _Transactions on Machine Learning Research_ , 2023.
  * Hopkins et al. (2020) Sam Hopkins, Jerry Li, and Fred Zhang.  Robust and heavy-tailed mean estimation made simple, via regret minimization.  _Advances in Neural Information Processing Systems_ , 33:11902–11912, 2020.
  * Howard et al. (2020) Steven R Howard, Aaditya Ramdas, Jon McAuliffe, and Jasjeet Sekhon.  Time-uniform Chernoff bounds via nonnegative supermartingales.  _Probability Surveys_ , 17:257–317, 2020.
  * Howard et al. (2021) Steven R Howard, Aaditya Ramdas, Jon McAuliffe, and Jasjeet Sekhon.  Time-uniform, nonparametric, nonasymptotic confidence sequences.  _The Annals of Statistics_ , 49(2), 2021.
  * Huber (1964) Peter J Huber.  Robust estimation of a location parameter.  _Annals of Mathematical Statistics_ , 35:73–101, 1964.
  * Huber (1965) Peter J Huber.  A robust version of the probability ratio test.  _The Annals of Mathematical Statistics_ , pages 1753–1758, 1965.
  * Huber (1973) Peter J Huber.  Robust regression: asymptotics, conjectures and Monte Carlo.  _The Annals of Statistics_ , pages 799–821, 1973.
  * Joly et al. (2017) Emilien Joly, Gábor Lugosi, and Roberto Imbuzeiro Oliveira.  On the estimation of the mean of a random vector.  _Electronic Journal of Statistics_ , 11(1):440, 2017.
  * Lai (1976) Tze Leung Lai.  On confidence sequences.  _The Annals of Statistics_ , pages 265–280, 1976.
  * Lai and Wei (1982) Tze Leung Lai and Ching Zong Wei.  Least squares estimates in stochastic regression models with applications to identification and control of dynamic systems.  _The Annals of Statistics_ , 10(1):154–166, 1982\.
  * Lindon et al. (2022) Michael Lindon, Chris Sanden, and Vaché Shirikian.  Rapid regression detection in software deployments through sequential testing.  In _Proceedings of the 28th ACM SIGKDD Conference on Knowledge Discovery and Data Mining_ , pages 3336–3346, 2022.
  * Lugosi and Mendelson (2019a) Gábor Lugosi and Shahar Mendelson.  Mean estimation and regression under heavy-tailed distributions: A survey.  _Foundations of Computational Mathematics_ , 19(5):1145–1190, 2019a.
  * Lugosi and Mendelson (2019b) Gábor Lugosi and Shahar Mendelson.  Sub-Gaussian estimators of the mean of a random vector.  _The Annals of Statistics_ , 47(2):783–794, 2019b.
  * Lugosi and Mendelson (2021) Gabor Lugosi and Shahar Mendelson.  Robust multivariate mean estimation: the optimality of trimmed mean.  _The Annals of Statistics_ , 2021.
  * Manole and Ramdas (2023) Tudor Manole and Aaditya Ramdas.  Martingale methods for sequential estimation of convex functionals and divergences.  _IEEE Transactions on Information Theory_ , 2023.
  * Mardia et al. (2000) Kanti V Mardia, Peter E Jupp, and KV Mardia.  _Directional statistics_ , volume 2.  Wiley Online Library, 2000.
  * Maurer and Pontil (2009) Andreas Maurer and Massimiliano Pontil.  Empirical bernstein bounds and sample variance penalization.  _arXiv preprint arXiv:0907.3740_ , 2009.
  * McAllester (1998) David A McAllester.  Some PAC-Bayesian theorems.  In _Proceedings of the eleventh annual conference on Computational learning theory_ , pages 230–234, 1998.
  * Prasad et al. (2020a) Adarsh Prasad, Sivaraman Balakrishnan, and Pradeep Ravikumar.  A robust univariate mean estimator is all you need.  In _International Conference on Artificial Intelligence and Statistics_ , pages 4034–4044. PMLR, 2020a.
  * Prasad et al. (2020b) Adarsh Prasad, Arun Sai Suggala, Sivaraman Balakrishnan, and Pradeep Ravikumar.  Robust estimation via robust gradient estimation.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 82(3):601–627, 2020b.
  * Quang (1985) Pham Xuan Quang.  Robust sequential testing.  _The Annals of Statistics_ , pages 638–649, 1985.
  * Shawe-Taylor and Williamson (1997) John Shawe-Taylor and Robert C Williamson.  A PAC analysis of a Bayesian estimator.  In _Proceedings of the Tenth Annual Conference on Computational Learning Theory_ , pages 2–9, 1997.
  * Shekhar and Ramdas (2023) Shubhanshu Shekhar and Aaditya Ramdas.  On the near-optimality of betting confidence sets for bounded means.  _arXiv preprint arXiv:2310.01547_ , 2023.
  * Tukey (1960) John Wilder Tukey.  A survey of sampling from contaminated distributions.  _Contributions to probability and statistics_ , pages 448–485, 1960\.
  * Valiant (1984) Leslie G Valiant.  A theory of the learnable.  _Communications of the ACM_ , 27(11):1134–1142, 1984.
  * Wang and Ramdas (2023a) Hongjian Wang and Aaditya Ramdas.  Catoni-style confidence sequences for heavy-tailed mean estimation.  _Stochastic Processes and their Applications_ , 2023a.
  * Wang and Ramdas (2023b) Hongjian Wang and Aaditya Ramdas.  Huber-robust confidence sequences.  In _International Conference on Artificial Intelligence and Statistics_ , pages 9662–9679. PMLR, 2023b.
  * Watamori (1996) Yoko Watamori.  Statistical inference of Langevin distribution for directional data.  _Hiroshima Mathematical Journal_ , 26(1):25–74, 1996.
  * Waudby-Smith and Ramdas (2023) Ian Waudby-Smith and Aaditya Ramdas.  Estimating means of bounded random variables by betting.  _Journal of the Royal Statistical Society, Series B (with discussion)_ , 2023.
  * Whitehouse et al. (2023) Justin Whitehouse, Steven Wu, and Aaditya Ramdas.  Time-uniform self-concentration for vector-valued processes.  _arXiv preprint arXiv:2310.09100_ , 2023.

##  Appendix A Stitched bounds with LIL Rate

Here we demonstrate how to modify our bounds in Section 2 to have the rate
designated by the law of the iterated logarithm (LIL), which is
$\sqrt{\log\log(t)/t}$. This is done via “stitching,” which is a common
technique in the literature on confidence sequences and originates in Howard
et al. (2021).

Stitching involves applying distinct bounds over geometrically spaced epochs.
Bounds in each epoch can be engineered to be tighter than a bound uniform over
all time. The bounds are then carefully combined together using a union bound.

To proceed, we define our “stitching function” as

$\ell(m)=(m+1)^{2}\zeta(2),\text{~{} where
~{}}\zeta(2)=\sum_{k=1}^{\infty}k^{-2}\approx 1.645.$ (17)

Note that $\sum_{m\geqslant 0}\frac{1}{\ell(m)}=1$. In what follows, any
function $\ell$ could be used such that $\sum_{m\geqslant
0}\frac{1}{\ell(m)}\leqslant 1$ but here we fix a particular choice for
convenience. As we did throughout Section 2, let us take
$\Theta=\mathbb{S}^{d-1}$ and we will place a vMF distribution over $\Theta$.
Let us consider a general $P$-supermartingale of the form

$S_{t}^{m}(\boldsymbol{\theta})=\prod_{i\leqslant
t}\exp\\{\lambda_{m}\langle\boldsymbol{\theta},\boldsymbol{X}_{i}-\boldsymbol{\mu}\rangle-\psi(\lambda_{m})v_{i}\\},$

where $v_{i}$ is some implicit variance process. For instance, in our
empirical-Bernstein bounds (Section 2.1) we have
$v_{i}=\frac{1}{2B}\norm{\boldsymbol{X}_{i}-\bar{\boldsymbol{\mu}}_{i-1}}^{2}$,
and in our sub-$\psi$ bounds (Section 2.2) we have $v_{i}=1$. Applying
Proposition 1 to the above process gives that for all $\alpha\in(0,1)$ and
$\lambda>0$, with probability $1-\alpha$, simultaneously for all $t\geqslant
1$ and $\boldsymbol{\vartheta}\in\mathbb{S}^{d-1}$,

$\displaystyle\lambda A_{d}(\kappa)\sum_{i\leqslant
t}\langle\boldsymbol{\vartheta},\boldsymbol{X}_{i}-\boldsymbol{\mu}\rangle$
$\displaystyle=\lambda\int_{\mathbb{S}^{d-1}}\sum_{i\leqslant
t}\langle\boldsymbol{\theta},\boldsymbol{X}_{i}-\boldsymbol{\mu}\rangle\gamma_{\boldsymbol{\vartheta}}(\mathrm{d}\boldsymbol{\theta})$
$\displaystyle\leqslant\psi(\lambda)\sum_{i\leqslant t}v_{i}+2\kappa
A_{d}(\kappa)+\log(1/\alpha).$

We will apply the above bound over the epochs $[2^{m},2^{m+1})$ for
$m=0,1,2,\dots$, using a distinct $\lambda=\lambda_{m}$ and
$\alpha=\alpha_{m}$ in each epoch. Consider the epoch $[2^{m},2^{m+1})$ for
some $m$. Set $V_{t}=\sum_{i\leqslant t}v_{i}$ and $r_{m}=\log(1/\alpha_{m})$.
Dividing the above bound by $\lambda A_{d}(\kappa)t$ on both sides gives that
with probability $1-\alpha_{m}$, simultaneously for all $t\in[2^{m},2^{m+1})$:

$\norm{\frac{1}{t}\sum_{i\leqslant
t}\boldsymbol{X}_{i}-\boldsymbol{\mu}}\leqslant\frac{\psi(\lambda_{m})V_{t}+2\kappa
A_{d}(\kappa)+r_{k}}{\lambda_{m}A_{d}(\kappa)t}=:g_{\lambda_{m},r_{m}}(V_{t}),\quad\forall
t\in[2^{m},2^{m+1}).$

Taking a union bound gives that with probability
$1-\sum_{m=0}^{\infty}\alpha_{m}$, simulaneously for all $t\geqslant 1$:

$\norm{\frac{1}{t}\sum_{i\leqslant
t}\boldsymbol{X}_{i}-\boldsymbol{\mu}}\leqslant
g_{\lambda_{m},r_{m}}(V_{t}),\text{~{} where
~{}}m\leqslant\log_{2}(t)\leqslant m+1.$

We take $\alpha_{m}=\alpha/\ell(m)$ so that
$\sum_{m=0}^{\infty}\alpha_{m}=\alpha$. The remaining work resides in ensuring
that $g_{\lambda_{m},r_{m}}(V_{t})$ shrinks at the desired iterated logarithm
rate. We handle this differently for the empirical-Bernstein bound and sub-
Gamma bounds.

###  A.1 Stitched Empirical-Bernstein Bound

Here we use the upper bound $\psi_{E}(\lambda)\leqslant\lambda^{2}$ for
$\lambda\leqslant 0.68$ and $V_{t}=\frac{1}{2B}\sum_{i\leqslant
t}\norm{\boldsymbol{X}_{i}-\bar{\boldsymbol{\mu}}_{i-1}}^{2}\leqslant\frac{1}{2B}\sum_{i\leqslant
t}(2B)^{2}=2tB$. Consider taking $\lambda_{m}=\sqrt{r_{m}/2^{m}}$. Then, for
large enough $m$, we have

$\displaystyle g_{\lambda_{m},r_{m}}(V_{t})$
$\displaystyle\leqslant\frac{\lambda_{m}^{2}2tB+2\kappa
A_{d}(\kappa)+r_{m}}{\lambda_{m}A_{d}(\kappa)t}\lesssim\frac{1}{A_{d}(\kappa)}\sqrt{\frac{r_{m}}{2^{m}}}+\frac{\kappa
2^{m/2}}{t\sqrt{r_{m}}}+\frac{\sqrt{r_{m}2^{m}}}{tA_{d}(\kappa)}.$

As we did throughout Section 2, take $\kappa=\sqrt{d}$ so that
$A_{d}(\kappa)\gtrsim 1/\sqrt{d}$. Take $r_{m}=\log(\ell(m)/\alpha)$. Notice
that since $t\in[2^{m},2^{m+1})$, we have $m\leqslant\log_{2}(t)$ so
$r_{m}\leqslant\log(\ell(\log_{2}(t))/\alpha)$. Also, $2^{m}\geqslant t/2$,
therefore

$\sqrt{\frac{r_{m}}{2^{m}}}\lesssim\sqrt{\frac{\log(\ell(\log_{2}t))+\log(1/\alpha)}{t}},\text{~{}
and ~{}}\sqrt{r_{m}2^{m}}\lesssim\sqrt{t\log(\ell(\log_{2}(t))/\alpha)}.$

It follows that

$\displaystyle g_{\lambda_{m},r_{m}}(V_{t})$
$\displaystyle\lesssim\sqrt{\frac{d(\log(\ell(\log_{2}(t)))+\log(1/\alpha))}{t}},$

which is the claimed iterated logarithm rate. Note that this analysis holds as
soon as $\lambda_{m}=\sqrt{r_{m}/2^{m}}\leqslant 0.68$, i.e., as soon as
$\log(\ell(m))+\log(1/\alpha)\leqslant(0.68)^{2}2^{m}$. In the case of
$\alpha=0.05$, this is satisfied for all $m\geqslant 4$.

###  A.2 Stitched Sub-Gamma Bound

Here we consider the sub-$\psi$ bounds of Section 2.2 with
$\psi(\lambda)=\psi_{G}(\lambda)=\frac{\lambda^{2}}{2(1-c\lambda)}$. This
generalizes both the sub-Gaussian and sub-exponential cases.

Note that in the sub-$\psi$ case we have $V_{t}=t$. We take
$\lambda_{m}=\psi^{-1}(r_{m}/2^{m})$, where

$\psi^{-1}(u)=\frac{2}{c+\sqrt{c^{2}+2/u}},$

so, noting as above that $r_{m}/2^{m}\asymp l/t$ where
$l=\log(\ell(\log_{2}(t))/\alpha)$, we have

$\psi^{-1}(r_{m}/2^{m})\asymp\frac{1}{c+\sqrt{c^{2}+t/l}}.$

Set $w_{m}=2\kappa A_{d}(\kappa)+r_{m}$ and rewrite the boundary as

$\displaystyle g_{\lambda_{m},r_{m}}(t)$
$\displaystyle\asymp\frac{\sqrt{d}r_{m}}{2^{m+1/2}\psi^{-1}(r_{m}/2^{m+1/2})}+\frac{\sqrt{d}+r_{m}}{\psi^{-1}(r_{m}/2^{m+1/2})t}$
$\displaystyle\lesssim\frac{\sqrt{d}l}{t\psi^{-1}(r_{m}/2^{m+1/2})}\lesssim\frac{l\sqrt{d}(c+\sqrt{c^{2}+t/l})}{t}$
$\displaystyle=\frac{\sqrt{d}(lc+\sqrt{l^{2}c^{2}+tl})}{t}\lesssim\frac{\sqrt{dtl}}{t}$
$\displaystyle=\sqrt{\frac{d(\log(\ell(\log_{2}(t)))+\log(1/\alpha))}{t}},$

which is the desired rate.

##  Appendix B A second (robust) time-uniform Catoni-Giulini bound

Now we give a time-uniform extension of a second Catoni-Giulini estimator for
heavy-tailed random vectors (Catoni and Giulini, 2017). The estimator is less
explicit than that presented in Section 3, and instead estimates the mean in
each direction $\boldsymbol{\vartheta}\in\mathbb{S}^{d-1}$. However, it lends
itself to a robust CSS, so we feel the result is worth presenting.

Let $\varphi$ be any symmetric influence function such that

$-\log(1-x+x^{2}/2)\leqslant\varphi(x)\leqslant\log(1+x+x^{2}/2),$

and $|\varphi(x)|\leqslant C$ for some constant $C>0$. For instance, Catoni
and Giulini (2017) consider $\varphi(x)=x-x^{3}/6$ for $-\sqrt{2}\leqslant
x\leqslant\sqrt{2}$, and $\varphi(x)=\text{sgn}(x)\frac{2\sqrt{2}}{3}$
otherwise. Consider the function

$\xi_{t}(\boldsymbol{\vartheta}):=\sum_{i\leqslant
t}\int\varphi(\lambda_{i}\langle\boldsymbol{\theta},\boldsymbol{X}_{i}\rangle)\phi_{\boldsymbol{\vartheta}}(\mathrm{d}\boldsymbol{\theta}).$
(18)

The function $\varphi$, and hence the estimator
$\xi_{t}(\boldsymbol{\theta})$, can be explicitly computed. We refer to Catoni
and Giulini (2017, Lemma 3.2) for details. The definition of $\varphi$ and
basic properties of the logarithm give us the following result.

######  Lemma 6.

The process
$M(\boldsymbol{\theta})\equiv(M_{t}(\boldsymbol{\theta}))_{t\geqslant 1}$
defined as

$\displaystyle M_{t}(\boldsymbol{\theta})=\prod_{i\leqslant
t}\frac{\exp(\varphi(\lambda_{i}\langle\boldsymbol{\theta},\boldsymbol{X}_{i}\rangle))}{1+\lambda_{i}\langle\boldsymbol{\theta},\boldsymbol{\mu}\rangle+\frac{\lambda_{i}^{2}}{2}\mathbb{E}_{i-1}[\langle\boldsymbol{\theta},\boldsymbol{X}_{i}\rangle^{2}]},$

where $\mathbb{E}_{i-1}[\cdot]=\mathbb{E}[\cdot|\mathcal{F}_{i-1}]$, is a
$P$-NSM.

###### Proof.

For all $\boldsymbol{\theta}\in\mathbb{R}^{d}$, observe that

$\displaystyle\mathbb{E}\left[\exp\\{\varphi(\lambda_{t}\langle\boldsymbol{\theta},\boldsymbol{X}_{i}\rangle)\\}|\mathcal{F}_{t-1}\right]$
$\displaystyle\leqslant\mathbb{E}\left[1+\lambda_{t}\langle\boldsymbol{\theta},\boldsymbol{X}_{t}\rangle+\frac{\lambda_{t}^{2}}{2}\langle\boldsymbol{\theta},\boldsymbol{X}_{t}\rangle^{2}\bigg{|}\mathcal{F}_{t-1}\right]$
$\displaystyle\leqslant
1+\lambda_{t}\langle\boldsymbol{\theta},\boldsymbol{\mu}\rangle+\frac{\lambda_{t}^{2}}{2}\mathbb{E}[\langle\boldsymbol{\theta},\boldsymbol{X}_{t}\rangle^{2}|\mathcal{F}_{t-1}].$

Dividing both sides by the quantity on the right gives the result. ∎

Taking $\phi_{\boldsymbol{0}}$ as our prior and applying Proposition 1 to the
process defined in Lemma 6 gives that with probability $1-\alpha$, for all
$t\geqslant 1$ and $\boldsymbol{\vartheta}\in\mathbb{S}^{d-1}$,

$\displaystyle\xi_{t}(\boldsymbol{\vartheta})$ $\displaystyle=\sum_{i\leqslant
t}\int\varphi(\lambda_{i}\langle\boldsymbol{\theta},\boldsymbol{X}_{i}\rangle)\phi_{\boldsymbol{\vartheta}}(\mathrm{d}\boldsymbol{\theta})$
$\displaystyle\leqslant\sum_{i\leqslant
t}\int\log(1+\lambda_{i}\langle\boldsymbol{\theta},\boldsymbol{\mu}\rangle+\frac{\lambda_{i}^{2}}{2}\mathbb{E}_{i-1}\langle\boldsymbol{\theta},\boldsymbol{X}_{i}\rangle^{2})\phi_{\boldsymbol{\vartheta}}(\mathrm{d}\boldsymbol{\theta})+\frac{\beta}{2}+\log(1/\alpha)$
$\displaystyle\leqslant\sum_{i\leqslant
t}\int\left\\{\lambda_{i}\langle\boldsymbol{\theta},\boldsymbol{\mu}\rangle+\frac{\lambda_{i}^{2}}{2}\mathbb{E}_{i-1}\langle\boldsymbol{\theta},\boldsymbol{X}_{i}\rangle^{2}\right\\}\phi_{\boldsymbol{\vartheta}}(\mathrm{d}\boldsymbol{\theta})+\frac{\beta}{2}+\log(1/\alpha)$
$\displaystyle=\sum_{i\leqslant
t}\lambda_{i}\langle\boldsymbol{\vartheta},\boldsymbol{\mu}\rangle+\sum_{i\leqslant
t}\frac{\lambda_{i}^{2}}{2}\bigg{(}\mathbb{E}_{i-1}\langle\boldsymbol{\vartheta},\boldsymbol{X}\rangle^{2}+\frac{\mathbb{E}_{i-1}\norm{\boldsymbol{X}}^{2}}{\beta}\bigg{)}+\frac{\beta}{2}+\log(1/\alpha)$
$\displaystyle\leqslant\sum_{i\leqslant
t}\lambda_{i}\langle\boldsymbol{\vartheta},\boldsymbol{\mu}\rangle+\sum_{i\leqslant
t}\frac{\lambda_{i}^{2}}{2}\bigg{(}v^{2/p}+\frac{v^{2/p}}{\beta}\bigg{)}+\frac{\beta}{2}+\log(1/\alpha).$

Therefore, with probability $1-\alpha$, for all $t\geqslant 1$ and all
$\boldsymbol{\vartheta}\in\mathbb{S}^{d-1}$,

$\xi_{t}(\boldsymbol{\vartheta})-\sum_{i\leqslant
t}\lambda_{i}\langle\boldsymbol{\vartheta},\boldsymbol{\mu}\rangle\leqslant\sum_{i\leqslant
t}\frac{\lambda_{i}^{2}}{2}\bigg{(}v^{2/p}+\frac{v^{2/p}}{\beta}\bigg{)}+\frac{\beta}{2}+\log\frac{1}{\alpha}.$

Taking $\beta=1$ and noting that $\xi$ is symmetric due to the (assumed)
symmetry of $\varphi$, the above display rearranges to

$\displaystyle\sup_{\boldsymbol{\vartheta}\in\mathbb{S}^{d-1}}\bigg{|}\xi_{t}(\boldsymbol{\vartheta})-\sum_{i\leqslant
t}\lambda_{i}\langle\boldsymbol{\vartheta},\boldsymbol{\mu}\rangle\bigg{|}\leqslant
v^{2/p}\sum_{i\leqslant t}\lambda_{i}^{2}+\frac{1}{2}+\log\frac{1}{\alpha},$

which in turn is

$\displaystyle\sup_{\boldsymbol{\vartheta}\in\mathbb{S}^{d-1}}\bigg{|}\frac{\xi_{t}(\boldsymbol{\vartheta})}{\sum_{i\leqslant
t}\lambda_{i}}-\langle\boldsymbol{\vartheta},\boldsymbol{\mu}\rangle\bigg{|}\leqslant\frac{v^{2/p}\sum_{i\leqslant
t}\lambda_{i}^{2}+\frac{1}{2}+\log\frac{1}{\alpha}}{\sum_{i\leqslant
t}\lambda_{i}}.$

We have proved the following theorem.

######  Theorem 5.

Let $\boldsymbol{X}_{1},\boldsymbol{X}_{2},\dots\sim P$ with mean
$\boldsymbol{\mu}$ obey $\mathbb{E}[\norm{\boldsymbol{X}}^{p}]\leqslant
v<\infty$ for some $p\geqslant 2$. Fix any $\beta>0$. Then, under Assumption
2, with probability at least $1-\alpha$, for all $t\geqslant 1$,

$\displaystyle\sup_{\boldsymbol{\vartheta}\in\mathbb{S}^{d-1}}\bigg{|}\frac{\xi_{t}(\boldsymbol{\vartheta})}{\sum_{i\leqslant
t}\lambda_{i}}-\langle\boldsymbol{\vartheta},\boldsymbol{\mu}\rangle\bigg{|}\leqslant\frac{v^{2/p}\sum_{i\leqslant
t}\lambda_{i}^{2}+\frac{1}{2}+\log\frac{1}{\alpha}}{\sum_{i\leqslant
t}\lambda_{i}}.$

We observe that, modulo constants, the bound has a similar form to that in
Theorem 3. The same discussion on optimizing $(\lambda_{t})$ thus applies
here.

We now demonstrate how to make this CSS robust. As in Section 4, this will
hinge on altering the underlying supermartingale, given by Lemma 6 in this
case.

Recall that by construction, the influence function $\varphi$ is bounded by
some constant $C$. By (14), this implies that

$\mathbb{E}_{Q}\exp\\{\varphi(\lambda_{t}\langle\boldsymbol{\theta},\boldsymbol{X}\rangle)\\}\leqslant\mathbb{E}_{P}\exp\\{\varphi(\lambda_{t}\langle\boldsymbol{\theta},\boldsymbol{X}\rangle)\\}+2C.$

This inequality enables the following result, which is in the same vein is
Lemma 5.

######  Lemma 7.

Let $(\lambda_{t})$ be a sequence of positive scalars and let $\epsilon>0$.
Suppose $|\varphi(x)|\leqslant C$ for all $x\in\mathbb{R}$. Then, for all
$\boldsymbol{\theta}\in\mathbb{R}^{d}$, the process
$M^{R}(\boldsymbol{\theta})\equiv(M_{t}^{R}(\boldsymbol{\theta}))_{t\geqslant
1}$ given by

$M_{t}^{R}(\boldsymbol{\theta})=\prod_{i\leqslant
t}\frac{\exp\\{\varphi(\lambda_{i}(\boldsymbol{\theta},\boldsymbol{X}_{i}\rangle))}{1+\lambda_{i}\langle\boldsymbol{\theta},\boldsymbol{\mu}\rangle+\frac{\lambda_{i}^{2}}{2}\mathbb{E}_{P}\langle\boldsymbol{\theta},\boldsymbol{X}\rangle^{2}+2C},$

is a $Q$-NSM for all $Q\in\mathscr{B}_{\operatorname{TV}}(P,\epsilon)$.

Note that the expectation on the denominator of
$M_{t}^{R}(\boldsymbol{\theta})$ is with respect to $P$. Following how Lemma 6
was leveraged to generate Theorem 5, we may leverage Lemma 7 to obtain the
following result.

######  Theorem 6.

Let $P$ be a distribution over
$\\{\boldsymbol{x}\in\mathbb{R}^{d}:\norm{\boldsymbol{x}}\leqslant B\\}$ with
mean $\boldsymbol{\mu}_{P}$. For $\epsilon>0$, let
$\boldsymbol{X}_{1},\boldsymbol{X}_{2},\dots\sim Q$ for
$Q\in\mathscr{B}_{\operatorname{TV}}(P,\epsilon)$. Then for any
$\alpha\in(0,1)$, with probability $1-\alpha$, for all $t\geqslant 1$,

$\displaystyle\small\sup_{\boldsymbol{\vartheta}\in\mathbb{S}^{d-1}}\bigg{|}\frac{\xi_{t}(\boldsymbol{\vartheta})}{\sum_{i\leqslant
t}\lambda_{i}}-\langle\boldsymbol{\vartheta},\boldsymbol{\mu}\rangle\bigg{|}\leqslant\frac{v^{2/p}\sum_{i\leqslant
t}\lambda_{i}^{2}+\frac{1}{2}+2C+\log\frac{1}{\alpha}}{\sum_{i\leqslant
t}\lambda_{i}}.$

Once again, the rate discussion is similar to Section 3 and is omitted here.
The proof of Theorem 6 is substantively similar to Theorem 5 and is omitted.

##  Appendix C A semi-empirical bound for heavy-tails

Here we present a semi-empirical bound for random vectors with heavy tails.
The rate is strictly larger than the Catoni-estimators, so we omitted the
result from the main paper. However, perhaps it demonstrates the possibility
of obtaining fully empirical bounds in the heavy-tailed setting, so we include
it here.

The bound is based on a supermartingale that can be derived using an
inequality furnished by Delyon (2009, Proposition 12); see also Howard et al.
(2020, Lemma 3). Variants of the resulting process have been studied by Wang
and Ramdas (2023a, Lemma 5) in the context of confidence sequences, and both
Haddouche and Guedj (2023, Lemma 1.3) and Chugg et al. (2023, Corollary 17) in
the context of PAC-Bayesian bounds. Throughout, let
$\Sigma=\text{Cov}(\boldsymbol{X})$.

######  Lemma 8.

Let
$Z_{i}(\boldsymbol{\theta})=\langle\boldsymbol{\theta},\boldsymbol{X}_{i}-\boldsymbol{\mu}\rangle$.
For each fixed $\boldsymbol{\theta}\in\mathbb{R}^{d}$, the process
$S(\boldsymbol{\theta})=(S_{t}(\boldsymbol{\theta}))_{t\geqslant 1}$ is a
nonnegative supermartingale, where

$S_{t}(\boldsymbol{\theta})=\prod_{i\leqslant
t}\exp\\{\lambda_{i}Z_{i}(\boldsymbol{\theta})-\frac{\lambda_{i}^{2}}{6}\left[Z_{i}^{2}(\boldsymbol{\theta})+2\norm{\boldsymbol{\theta}}^{2}_{2}\Tr(\Sigma)\right]\bigg{\\}}$
(19)

###### Proof.

Delyon (2009, Equation (52)) demonstrates that for all $x\in\mathbb{R}$,
$\exp(x-x^{2}/2)\leqslant 1+x+x^{2}/3$. Set
$x=\lambda_{i}Z_{i}(\boldsymbol{\theta})$. Applying expectations to this
inequality, we have

$\displaystyle\mathbb{E}\left[\exp\left\\{\lambda_{i}Z_{i}(\boldsymbol{\theta})-\frac{\lambda_{i}^{2}}{2}Z_{i}^{2}(\boldsymbol{\theta})\right\\}\right]$
$\displaystyle\leqslant
1+\lambda_{i}\mathbb{E}[Z_{i}(\boldsymbol{\theta})]+\frac{\lambda_{i}^{2}}{3}\mathbb{E}[Z_{i}^{2}(\boldsymbol{\theta})]$
$\displaystyle=1+\frac{\lambda_{i}^{2}}{3}\mathbb{E}[Z_{i}^{2}(\boldsymbol{\theta})]\leqslant\exp\left\\{\frac{\lambda_{i}^{2}}{3}\mathbb{E}[Z_{i}^{2}(\boldsymbol{\theta})]\right\\}.$

Cauchy-Schwarz gives
$\mathbb{E}Z_{i}^{2}(\boldsymbol{\theta})\leqslant\norm{\boldsymbol{\theta}}_{2}^{2}\mathbb{E}\norm{\boldsymbol{X}_{i}-\boldsymbol{\mu}}^{2}=\norm{\boldsymbol{\theta}}^{2}\Tr(\Sigma)$.
Applying this to the display above and rearranging yields that

$\exp\\{\lambda_{i}Z_{i}(\boldsymbol{\theta})-\frac{\lambda_{i}^{2}}{2}Z_{i}^{2}(\boldsymbol{\theta})-\frac{\lambda_{i}^{2}}{3}\norm{\boldsymbol{\theta}}^{2}\Tr(\Sigma)\\}\leqslant
1,$

which in turn implies that $S(\boldsymbol{\theta})$ is a supermartingale. ∎

Applying Proposition 1 to Lemma 8 with the prior $\phi_{\boldsymbol{0}}$, we
obtain that with probability at least $1-\alpha$, for all $t\geqslant 1$ and
all $\boldsymbol{\vartheta}\in\mathbb{S}^{d-1}$,

$\sum_{i\leqslant
t}\lambda_{i}\mathbb{E}_{\boldsymbol{\theta}}Z_{i}(\boldsymbol{\theta})\leqslant\frac{1}{6}\sum_{i\leqslant
t}\lambda_{i}^{2}\mathbb{E}_{\boldsymbol{\theta}}Z_{i}^{2}(\boldsymbol{\theta})+\frac{1}{3}\mathbb{E}_{\boldsymbol{\theta}}\norm{\boldsymbol{\theta}}^{2}\Tr(\Sigma)\sum_{i\leqslant
t}\lambda_{i}^{2}+\frac{\beta}{2}+\log(1/\alpha).$ (20)

We upper bound the right hand side by writing
$\mathbb{E}_{\rho}Z_{i}^{2}(\boldsymbol{\theta})\leqslant\norm{\boldsymbol{\theta}}^{2}\norm{\boldsymbol{X}_{i}-\boldsymbol{\mu}}^{2}\leqslant\norm{\boldsymbol{\theta}}^{2}(\norm{\boldsymbol{X}_{i}}+v^{2/p})^{2}$.
Then, note that
$\mathbb{E}_{\boldsymbol{\theta}\sim\phi_{\boldsymbol{\vartheta}}}\norm{\boldsymbol{\theta}}^{2}=d/\beta^{2}$.
Therefore, (20) becomes

$\sup_{\boldsymbol{\vartheta}\in\mathbb{S}^{d-1}}\sum_{i\leqslant
t}\lambda_{i}Z_{i}(\boldsymbol{\vartheta})\leqslant\frac{d}{6\beta^{2}}\sum_{i\leqslant
t}\lambda_{i}^{2}\norm{\boldsymbol{X}_{i}-\boldsymbol{\mu}}^{2}+\frac{d\Tr(\Sigma)}{3\beta^{2}}\sum_{i\leqslant
t}\lambda_{i}^{2}+\beta/2+\log(1/\alpha).$ (21)

Consider taking $\beta=\sqrt{d}$ so that dependence on $d$ disappears from the
first two terms on the RHS. Rearranging yields the following theorem which,
owing to the reliance on the norm of the observations, can be considered
another empirical Bernstein bound.

######  Theorem 7.

Let $\boldsymbol{X}_{1},\boldsymbol{X}_{2},\dots\overset{\mathrm{iid}}{\sim}P$
with mean $\boldsymbol{\mu}$ obey (7). Let $(\lambda_{t})$ be a sequence of
positive scalars. Then, with probability $1-\alpha$, for all $t\geqslant 1$,

$\norm{\frac{\sum_{i\leqslant
t}\lambda_{i}\boldsymbol{X}_{i}}{\sum_{i\leqslant
t}\lambda_{i}}-\boldsymbol{\mu}}\leqslant\frac{\sum_{i\leqslant
t}\lambda_{i}^{2}(\norm{\boldsymbol{X}_{i}}_{2}^{2}+v^{2/p})}{6\sum_{i\leqslant
t}\lambda_{i}}+\frac{\Tr(\Sigma)\sum_{i\leqslant
t}\lambda_{i}^{2}}{3\sum_{i\leqslant
t}\lambda_{i}}+\frac{\sqrt{d}/2+\log(1/\alpha)}{\sum_{i\leqslant
t}\lambda_{i}}.$

Consider setting $\lambda_{t}\asymp\sqrt{\frac{\log(1/\alpha)}{v^{2/p}t\log
t}}$. Then, noting that $\sum_{i\leqslant t}\frac{1}{i\log
i}\asymp\log\log(t)$ and $\sum_{i\leqslant t}\frac{1}{\sqrt{i\log
i}}\asymp\sqrt{t/\log(t)}$, the _expected_ width of the RHS above can be
bounded as

$\displaystyle\widetilde{O}\left(\sqrt{\frac{v^{2/p}\log(t)\log(1/\alpha)}{t}}\max\left\\{1,\frac{\sqrt{d}}{\log(1/\alpha)}\right\\}\right),$

where $\widetilde{O}$ hides iterated logarithm factors. The expected width
refers to the fact that we’ve applied expectations over $P$ and bounded
$\mathbb{E}\norm{\boldsymbol{X}_{i}}^{2}$ as $v^{2/p}$.

##  Appendix D Omitted proofs

###  D.1 Section 2

#### Proof of Lemma 1

Fan et al. (2015, Equation (4.11)) demonstrates that for all
$\upsilon\geqslant-1$ and $\lambda\in[0,1)$,

$\exp\\{\lambda\upsilon-\psi_{E}(\lambda)\upsilon^{2}\\}\leqslant
1+\lambda\upsilon,$ (22)

where we recall that
$\psi_{E}(\lambda)=|\log(1-\lambda)+\lambda|=-\lambda-\log(1-\lambda)$. In
order to demonstrate that $(N_{t}(\boldsymbol{\theta}))$ is a $P$-NSM, it
suffices to demonstrate that

$\small\mathbb{E}_{P}\left[\exp\left\\{\frac{\lambda_{t}}{2B}\langle\boldsymbol{\theta},\boldsymbol{X}_{t}-\boldsymbol{\mu}\rangle-\frac{\psi_{E}(\lambda_{t})}{(2B)^{2}}\langle\boldsymbol{\theta},\boldsymbol{X}_{t}-\bar{\boldsymbol{\mu}}_{t-1}\rangle^{2}\right\\}\bigg{|}\mathcal{F}_{t-1}\right]\leqslant
1,$

for all $t\geqslant 1$. Inspired by a trick of Howard et al. (2021, Section
A.8), set

$\boldsymbol{Y}_{t}:=\frac{\boldsymbol{X}_{t}-\boldsymbol{\mu}}{2B},\quad\boldsymbol{\delta}_{t}=\frac{\bar{\boldsymbol{\mu}}_{t}-\boldsymbol{\mu}}{2B},$

and observe that
$\boldsymbol{Y}_{t}-\boldsymbol{\delta}_{t-1}=\frac{1}{2B}(\boldsymbol{X}_{t}-\bar{\boldsymbol{\mu}}_{t-1})$.
Then,

$\displaystyle\quad\exp\left\\{\frac{\lambda_{t}}{2B}\langle\boldsymbol{\theta},\boldsymbol{X}_{t}-\boldsymbol{\mu}\rangle-\frac{\psi_{E}(\lambda_{t})}{(2B)^{2}}\langle\boldsymbol{\theta},\boldsymbol{X}_{t}-\bar{\boldsymbol{\mu}}_{t-1}\rangle^{2}\right\\}$
$\displaystyle=\exp\left\\{\lambda_{t}\langle\boldsymbol{\theta},\boldsymbol{Y}_{t}\rangle-\psi_{E}(\lambda_{t})\langle\boldsymbol{\theta},\boldsymbol{Y}_{t}-\boldsymbol{\delta}_{t-1}\rangle^{2}\right\\}$
$\displaystyle=\exp\left\\{\lambda_{t}\langle\boldsymbol{\theta},\boldsymbol{Y}_{t}-\boldsymbol{\delta}_{t-1}\rangle-\psi_{E}(\lambda_{t})\langle\boldsymbol{\theta},\boldsymbol{Y}_{t}-\boldsymbol{\delta}_{t-1}\rangle^{2}\right\\}\exp(\lambda_{t}\langle\boldsymbol{\theta},\boldsymbol{\delta}_{t-1}\rangle).$

Since $\norm{\boldsymbol{\theta}}=1$ we have, by construction,

$\displaystyle|\langle\boldsymbol{\theta},\boldsymbol{Y}_{t}-\boldsymbol{\delta}_{t-1}\rangle|$
$\displaystyle\leqslant\frac{1}{2B}(\norm{\boldsymbol{X}_{t}}+\norm{\bar{\boldsymbol{\mu}}_{t-1}})\leqslant
1,$

since $\bar{\boldsymbol{\mu}}_{t-1}$ is an average of vectors each with norm
at most $B$. Therefore, we may apply (22) with
$\upsilon=\langle\boldsymbol{\theta},\boldsymbol{Y}_{t}-\boldsymbol{\delta}_{t-1}\rangle$,
which, paired with the display above yields

$\displaystyle\exp\left\\{\frac{\lambda_{t}}{2B}\langle\boldsymbol{\theta},\boldsymbol{X}_{t}-\boldsymbol{\mu}\rangle-\frac{\psi_{E}(\lambda_{t})}{(2B)^{2}}\langle\boldsymbol{\theta},\boldsymbol{X}_{t}-\bar{\boldsymbol{\mu}}_{t-1}\rangle^{2}\right\\}\leqslant(1+\lambda_{t}\langle\boldsymbol{\theta},\boldsymbol{Y}_{t}-\boldsymbol{\delta}_{t-1}\rangle)\exp(\lambda_{t}\langle\boldsymbol{\theta},\boldsymbol{\delta}_{t-1}\rangle).$

Taking expectations and noticing that
$\mathbb{E}_{P}[\boldsymbol{Y}_{t}|\mathcal{F}_{t-1}]=0$ gives

$\displaystyle\mathbb{E}_{P}\left[\exp\left\\{\frac{\lambda_{t}}{2B}\langle\boldsymbol{\theta},\boldsymbol{X}_{t}-\boldsymbol{\mu}\rangle-\frac{\psi_{E}(\lambda_{t})}{(2B)^{2}}\langle\boldsymbol{\theta},\boldsymbol{X}_{t}-\bar{\boldsymbol{\mu}}_{t-1}\rangle^{2}\right\\}\bigg{|}\mathcal{F}_{t-1}\right]$
$\displaystyle\leqslant(1-\lambda_{t}\langle\boldsymbol{\theta},\boldsymbol{\delta}_{t-1}\rangle)\exp(\lambda_{t}\langle\boldsymbol{\theta},\boldsymbol{\delta}_{t-1}\rangle)\leqslant
1,$

where the final inequality uses the fact that $1+x\leqslant e^{x}$ for all
$x\in\mathbb{R}$. This completes the proof.

#### Proof of Theorem 1

Recall that Lemma 1 proves the process defined by the following quantity is a
$P$-NSM:

$N_{t}(\boldsymbol{\theta})=\prod_{i\leqslant
t}\exp\left\\{\frac{\lambda_{i}}{2B}\langle\boldsymbol{\theta},\boldsymbol{X}_{i}-\boldsymbol{\mu}\rangle-\frac{\psi_{E}(\lambda_{i})}{(2B)^{2}}\langle\boldsymbol{\theta},\boldsymbol{X}_{i}-\bar{\boldsymbol{\mu}}_{i-1}\rangle^{2}\right\\}.$

Applying Proposition 1 to $N_{t}(\boldsymbol{\theta})$ gives that with
probability $1-\alpha$, for all $t\geqslant 1$ and all
$\gamma_{\boldsymbol{\vartheta}}$:

$\displaystyle\int_{\mathbb{S}^{d-1}}\sum_{i\leqslant
t}\frac{\lambda_{i}}{2B}\langle\boldsymbol{\theta},\boldsymbol{X}_{i}-\boldsymbol{\mu}\rangle\gamma_{\boldsymbol{\vartheta}}(\mathrm{d}\boldsymbol{\theta})$
$\displaystyle\leqslant\int_{\mathbb{S}^{d-1}}\sum_{i\leqslant
t}\frac{\psi_{E}(\lambda_{i})}{(2B)^{2}}\langle\boldsymbol{\theta},\boldsymbol{X}_{i}-\bar{\boldsymbol{\mu}}_{i-1}\rangle^{2}\gamma_{\boldsymbol{\vartheta}}(\mathrm{d}\boldsymbol{\theta})+D_{\operatorname{KL}}(\gamma_{\boldsymbol{\vartheta}}\|\gamma_{\boldsymbol{1}})+\log(1/\alpha)$
$\displaystyle\leqslant\int_{\mathbb{S}^{d-1}}\sum_{i\leqslant
t}\frac{\psi_{E}(\lambda_{i})}{(2B)^{2}}\norm{\boldsymbol{X}_{i}-\bar{\boldsymbol{\mu}}_{i-1}}^{2}\gamma_{\boldsymbol{\vartheta}}(\mathrm{d}\boldsymbol{\theta})+2\kappa
A_{d}(\kappa)+\log(1/\alpha)$ $\displaystyle=\sum_{i\leqslant
t}\frac{\psi_{E}(\lambda_{i})}{(2B)^{2}}\norm{\boldsymbol{X}_{i}-\bar{\boldsymbol{\mu}}_{i-1}}^{2}+2\kappa
A_{d}(\kappa)+\log(1/\alpha).$

Integrating the left hand side above gives

$\int_{\mathbb{S}^{d-1}}\sum_{i\leqslant
t}\frac{\lambda_{i}}{2B}\langle\boldsymbol{\theta},\boldsymbol{X}_{i}-\boldsymbol{\mu}\rangle\gamma_{\boldsymbol{\vartheta}}(\mathrm{d}\boldsymbol{\theta})=A_{d}(\kappa)\sum_{i\leqslant
t}\frac{\lambda_{i}}{2B}\langle\boldsymbol{\vartheta},\boldsymbol{X}_{i}-\boldsymbol{\mu}\rangle.$

Multiplying both sides by $2B$ and dividing by $A_{d}(\kappa)$, we have shown
that with probability $1-\alpha$, for all $t\geqslant 1$,

$\sup_{\boldsymbol{\vartheta}\in\mathbb{S}^{d-1}}\bigg{\langle}\boldsymbol{\vartheta},\sum_{i\leqslant
t}(\boldsymbol{X}_{i}-\boldsymbol{\mu})\bigg{\rangle}\leqslant\frac{\frac{1}{2B}\sum_{i\leqslant
t}\psi_{E}(\lambda_{i})\norm{\boldsymbol{X}_{i}-\bar{\boldsymbol{\mu}}_{i-1}}^{2}+4B\kappa
A_{d}(\kappa)+2B\log(1/\alpha)}{A_{d}(\kappa)}.$

The left hand side is equal to $\norm{\sum_{i\leqslant
t}\lambda_{i}(\boldsymbol{X}_{i}-\boldsymbol{\mu})}$. Dividing by
$\sum_{i}\lambda_{i}$ results in the bound:

$\norm{\frac{\sum_{i\leqslant
t}\lambda_{i}\boldsymbol{X}_{i}}{\sum_{i\leqslant
t}\lambda_{i}}-\boldsymbol{\mu}}\leqslant\frac{\frac{1}{2B}\sum_{i\leqslant
t}\psi_{E}(\lambda_{i})\norm{\boldsymbol{X}_{i}-\bar{\boldsymbol{\mu}}_{i-1}}^{2}+4B\kappa
A_{d}(\kappa)+2B\log(1/\alpha)}{A_{d}(\kappa)\sum_{i\leqslant t}\lambda_{i}},$

which is the desired result.

###  D.2 Section 3

#### Proof of Lemma 2

First, let us observe the relationship

$0\leqslant 1-\frac{a\wedge 1}{a}\leqslant a,\quad\forall a>0.$

This is easily seen by case analysis. Indeed, for $a\geqslant 1$, we have
$(a\wedge 1)/a=1/a$ and $1-1/a\leqslant 1\leqslant a$. For $a<1$, we have
$1-(a\wedge 1)/a=1-1=0\leqslant a$.

Now, let

$\alpha(\boldsymbol{X})=\frac{\lambda\norm{\boldsymbol{X}}\wedge
1}{\lambda\norm{\boldsymbol{X}}},$

and note that the above analysis demonstrates that

$|\alpha(\boldsymbol{X})-1|=1-\alpha(\boldsymbol{X})\leqslant\lambda\norm{\boldsymbol{X}}.$
(23)

Therefore,

$\displaystyle\quad\langle\boldsymbol{\vartheta},\boldsymbol{\mu}_{t}-\boldsymbol{\mu}\rangle$
$\displaystyle=\langle\boldsymbol{\vartheta},\mathbb{E}[\alpha(\boldsymbol{X})\boldsymbol{X}]-\mathbb{E}[\boldsymbol{X}]\rangle$
$\displaystyle=\mathbb{E}(\alpha(\boldsymbol{X})-1)\langle\boldsymbol{\vartheta},\boldsymbol{X}\rangle$
$\displaystyle\leqslant\mathbb{E}|\alpha(\boldsymbol{X})-1||\langle\boldsymbol{\vartheta},\boldsymbol{X}\rangle|$
$\displaystyle\leqslant\mathbb{E}\lambda\norm{\boldsymbol{X}}\norm{\boldsymbol{\vartheta}}\norm{\boldsymbol{X}}$
by (23) and CS $\displaystyle=\lambda\mathbb{E}\norm{\boldsymbol{X}}^{2}$
$\displaystyle\norm{\boldsymbol{\vartheta}}=1$ $\displaystyle\leqslant\lambda
v^{2/p}$ $\displaystyle\text{by Jensen}.$

This proves the claim.

#### Proof of Lemma 4

Recall that we are trying to prove that

$\displaystyle\operatorname*{\mathbb{E}}_{\boldsymbol{\theta}\sim\phi_{\boldsymbol{\vartheta}}}\log\operatorname*{\mathbb{E}}_{\boldsymbol{X}\sim
P}\left\\{\mathrm{e}^{\left\langle\boldsymbol{\theta},\mathsf{th}_{t}(\boldsymbol{X})-\boldsymbol{\mu}_{t}\right\rangle}\right\\}\leqslant\frac{1}{2}v^{2/p}\lambda_{i}^{2}e^{2/\beta+2}$

To begin, notice that Jensen’s inequality gives

$\displaystyle\quad\int\log\operatorname*{\mathbb{E}}_{\boldsymbol{X}\sim
P}\left\\{\mathrm{e}^{\lambda_{i}\left\langle\boldsymbol{\theta},\mathsf{th}_{i}(\boldsymbol{X})-\boldsymbol{\mu}_{i}\right\rangle}\right\\}\rho_{\boldsymbol{\vartheta}}(\mathrm{d}\boldsymbol{\theta})$
$\displaystyle\leqslant\log\int\operatorname*{\mathbb{E}}_{\boldsymbol{X}\sim
P}\left\\{\mathrm{e}^{\lambda_{i}\left\langle\boldsymbol{\theta},\mathsf{th}_{i}(\boldsymbol{X})-\boldsymbol{\mu}_{i}\right\rangle}\right\\}\rho_{\boldsymbol{\vartheta}}(\mathrm{d}\boldsymbol{\theta})$
$\displaystyle=\log\operatorname*{\mathbb{E}}_{\boldsymbol{X}\sim
P}\left\\{\int\mathrm{e}^{\lambda_{i}\left\langle\boldsymbol{\theta},\mathsf{th}_{i}(\boldsymbol{X})-\boldsymbol{\mu}_{i}\right\rangle}\rho_{\boldsymbol{\vartheta}}(\mathrm{d}\boldsymbol{\theta})\right\\}$
$\displaystyle=\log\mathbb{E}\exp\left(\lambda_{i}\left\langle\boldsymbol{\vartheta},\mathsf{th}_{i}(\boldsymbol{X})-\boldsymbol{\mu}_{i}\right\rangle+\frac{\lambda_{i}^{2}}{2\beta}\|\mathsf{th}_{i}(\boldsymbol{X})-\boldsymbol{\mu}_{i}\|^{2}\right),$

where the final line uses the usual closed-form expression of the multivariate
Gaussian MGF.

Define the functions on $\mathbb{R}$

$g_{1}(x):=\frac{1}{x}(e^{x}-1),$

and

$g_{2}(x):=\frac{2}{x^{2}}(e^{x}-x-1),$

(with $g_{1}(0)=g_{2}(0)=1$ by continuous extension). Both $g_{1}$ and $g_{2}$
are increasing. Notice that

$e^{x+y}=1+x+\frac{x^{2}}{2}g_{2}(x)+g_{1}(y)ye^{x}.$ (24)

Consider setting $x$ and $y$ to be the two terms in the CGF above, i.e.,

$x=\lambda_{i}\langle\boldsymbol{\vartheta},\mathsf{th}_{i}(\boldsymbol{X})-\boldsymbol{\mu}_{i}\rangle,$
$y=\frac{\lambda_{i}^{2}}{2\beta}\|\mathsf{th}_{i}(\boldsymbol{X})-\boldsymbol{\mu}_{i}\|^{2},$

where we recall that $\boldsymbol{\vartheta}\in\mathbb{S}^{d-1}$. By before
applying (24) we would like to develop upper bounds on $x$ and $y$. Observe
that

$\norm{\mathsf{th}_{i}(\boldsymbol{X})}=\frac{\lambda_{i}\norm{\boldsymbol{X}}\wedge
1}{\lambda_{i}}\leqslant\frac{1}{\lambda_{i}},$

and consequently,

$\norm{\boldsymbol{\mu}_{i}}=\norm{\mathbb{E}\mathsf{th}_{i}(\boldsymbol{X})}\leqslant\mathbb{E}\norm{\mathsf{th}_{i}(\boldsymbol{X})}\leqslant\frac{1}{\lambda_{i}},$

by Jensen’s inequality. Therefore, by Cauchy-Schwarz and the triangle
inequality,

$x\leqslant\lambda_{i}\norm{\boldsymbol{\vartheta}}\norm{\mathsf{th}_{i}(\boldsymbol{X})-\boldsymbol{\mu}_{i}}\leqslant\lambda_{i}(\norm{\mathsf{th}_{i}(\boldsymbol{X})}+\norm{\boldsymbol{\mu}_{i}})\leqslant
2.$

Via similar reasoning, we can bound $y$ as

$\displaystyle
y\leqslant\frac{\lambda_{i}^{2}}{2\beta}(\norm{\mathsf{th}_{i}(\boldsymbol{X})}+\norm{\boldsymbol{\mu}_{i}})^{2}\leqslant\frac{2}{\beta}.$

Finally, substituting in these values of $x$ and $y$ to (24), taking
expectations, and using the fact that $g_{1}$ and $g_{2}$ are increasing gives

$\displaystyle\quad\mathbb{E}\left[\exp\left(\left\langle\boldsymbol{\vartheta},\mathsf{th}_{i}(\boldsymbol{X})-\boldsymbol{\mu}_{i}\right\rangle+\frac{1}{2\beta}\|\mathsf{th}_{i}(\boldsymbol{X})-\boldsymbol{\mu}_{i}\|^{2}\right)\bigg{|}\mathcal{F}_{i-1}\right]$
$\displaystyle\leqslant
1+\lambda_{i}\mathbb{E}[\langle\boldsymbol{\vartheta},\mathsf{th}_{i}(\boldsymbol{X})-\boldsymbol{\mu}_{i}\rangle|\mathcal{F}_{i-1}]+g_{2}\left(2\right)\frac{\lambda_{i}^{2}}{2}\mathbb{E}[\langle\boldsymbol{\vartheta},\mathsf{th}_{i}(\boldsymbol{X})-\boldsymbol{\mu}_{i}\rangle^{2}|\mathcal{F}_{i-1}]$
$\displaystyle\qquad+g_{1}\left(\frac{2}{\beta}\right)\frac{\lambda_{i}^{2}e^{2}}{2\beta}\mathbb{E}[\norm{\mathsf{th}_{i}(\boldsymbol{X})-\boldsymbol{\mu}_{i}}^{2}|\mathcal{F}_{i-1}]$
$\displaystyle\leqslant
1+g_{2}\left(2\right)\frac{\lambda_{i}^{2}}{2}\mathbb{E}[\norm{\mathsf{th}_{i}(\boldsymbol{X})-\boldsymbol{\mu}_{i}}^{2}|\mathcal{F}_{i-1}]+g_{1}\left(\frac{2}{\beta}\right)\frac{\lambda_{i}^{2}e^{2}}{2\beta}\mathbb{E}[\norm{\mathsf{th}_{i}(\boldsymbol{X})-\boldsymbol{\mu}_{i}}^{2}|\mathcal{F}_{i-1}],$

where we’ve used that
$\mathbb{E}[\mathsf{th}_{i}(\boldsymbol{X})-\boldsymbol{\mu}_{i}|\mathcal{F}_{i-1}]=\boldsymbol{0}$.
Denote by $\boldsymbol{X}^{\prime}$ an iid copy of $\boldsymbol{X}$. Using the
notation $\mathbb{E}_{i-1}[\cdot]=\mathbb{E}[\cdot|\mathcal{F}_{i-1}]$, we can
bound the norm as follows:

$\displaystyle\mathbb{E}_{i-1}[\norm{\mathsf{th}_{i}(\boldsymbol{X})-\boldsymbol{\mu}_{i}}^{2}]$
$\displaystyle=\mathbb{E}_{i-1}[\norm{\mathsf{th}_{i}(\boldsymbol{X})}^{2}]-\norm{\boldsymbol{\mu}_{i}}^{2}$
$\displaystyle=\frac{1}{2}\mathbb{E}_{i-1}\left[\norm{\mathsf{th}_{i}(\boldsymbol{X})}^{2}-2\langle\mathsf{th}_{i}(\boldsymbol{X}),\boldsymbol{\mu}_{i}\rangle+\mathbb{E}_{i-1}[\norm{\mathsf{th}_{i}(\boldsymbol{X})}^{2}]\right]$
$\displaystyle=\frac{1}{2}\operatorname*{\mathbb{E}}_{\boldsymbol{X}}\bigg{[}\norm{\mathsf{th}_{i}(\boldsymbol{X})}^{2}-2\left\langle\mathsf{th}_{i}(\boldsymbol{X}),\operatorname*{\mathbb{E}}_{\boldsymbol{X}^{\prime}}[\mathsf{th}_{i}(\boldsymbol{X}^{\prime})|\mathcal{F}_{i-1}]\right\rangle$
$\displaystyle+\qquad\operatorname*{\mathbb{E}}_{\boldsymbol{X}^{\prime}}\left[\norm{\mathsf{th}_{i}(\boldsymbol{X}^{\prime})}^{2}|\mathcal{F}_{i-1}\right]\bigg{|}\mathcal{F}_{i-1}\bigg{]}$
$\displaystyle=\frac{1}{2}\operatorname*{\mathbb{E}}_{\boldsymbol{X},\boldsymbol{X}^{\prime}}\left[\norm{\mathsf{th}_{i}(\boldsymbol{X})-\mathsf{th}_{i}(\boldsymbol{X}^{\prime})}^{2}|\mathcal{F}_{i-1}\right]$
$\displaystyle\leqslant\frac{1}{2}\operatorname*{\mathbb{E}}_{\boldsymbol{X},\boldsymbol{X}^{\prime}}\left[\norm{\boldsymbol{X}-\boldsymbol{X}^{\prime}}^{2}|\mathcal{F}_{i-1}\right]$
$\displaystyle=\mathbb{E}_{i-1}\norm{\boldsymbol{X}-\mathbb{E}_{i-1}\boldsymbol{X}}^{2}\leqslant\mathbb{E}_{i-1}\norm{\boldsymbol{X}}^{2}\leqslant
v^{2/p},$

where the first inequality uses the basic fact from convex analysis that, for
a closed a convex set $D\subset\mathbb{R}^{n}$ and any
$\mathbf{x},\mathbf{y}\in\mathbb{R}^{n}$,

$\norm{\Pi_{D}(\mathbf{x})-\Pi_{D}(\mathbf{y})}\leqslant\norm{\mathbf{x}-\mathbf{y}},$

where $\Pi_{D}$ is the projection onto $D$. Putting everything together thus
far, we have shown that

$\displaystyle\int\log\operatorname*{\mathbb{E}}_{\boldsymbol{X}\sim
P}\left\\{\mathrm{e}^{\lambda_{i}\left\langle\boldsymbol{\theta},\mathsf{th}_{i}(\boldsymbol{X})-\boldsymbol{\mu}_{i}\right\rangle}\bigg{|}\mathcal{F}_{i-1}\right\\}\rho_{\boldsymbol{\vartheta}}(\mathrm{d}\boldsymbol{\theta})$
$\displaystyle\leqslant\log\left\\{1+g_{2}\left(2\right)\frac{\lambda_{i}^{2}v^{2/p}}{2}+g_{1}\left(\frac{2}{\beta}\right)\frac{\lambda_{i}^{2}e^{2}}{2\beta}v^{2/p}\right\\}$
$\displaystyle\leqslant
g_{2}\left(2\right)\frac{\lambda_{i}^{2}v^{2/p}}{2}+g_{1}\left(\frac{2}{\beta}\right)\frac{\lambda_{i}^{2}e^{2}}{2\beta}v^{2/p}$
$\displaystyle=v^{2/p}\frac{\lambda_{i}^{2}}{4}\left\\{e^{2/\beta+2}-3\right\\}$
$\displaystyle\leqslant\frac{1}{4}v^{2/p}\lambda_{i}^{2}e^{2/\beta+2},$

which is the desired inequality.

#### Proof of Theorem 3

We simply apply Lemmas 4 and 2. For all
$\boldsymbol{\vartheta}\in\mathbb{S}^{d-1}$,

$\displaystyle\sum_{i\leqslant
t}\lambda_{i}\langle\boldsymbol{\vartheta},\mathsf{th}_{i}(\boldsymbol{X}_{i})-\boldsymbol{\mu}\rangle$
$\displaystyle=\sum_{i\leqslant
t}\lambda_{i}(\langle\boldsymbol{\vartheta},\mathsf{th}_{i}(\boldsymbol{X}_{i})-\boldsymbol{\mu}_{i}\rangle+\langle\boldsymbol{\vartheta},\boldsymbol{\mu}_{i}-\boldsymbol{\mu}\rangle)$
$\displaystyle\leqslant\sum_{i\leqslant
t}\lambda_{i}\langle\boldsymbol{\vartheta},\mathsf{th}_{i}(\boldsymbol{X}_{i})-\boldsymbol{\mu}_{i}\rangle+v^{2/p}\sum_{i\leqslant
t}\lambda_{i}^{2}$
$\displaystyle\leqslant\frac{v^{2/p}e^{\frac{2}{\beta}+2}}{4}\sum_{i\leqslant
t}\lambda_{i}^{2}+v^{2/p}\sum_{i\leqslant
t}\lambda_{i}^{2}+\frac{\beta}{2}+\log(\frac{1}{\alpha})$
$\displaystyle\leqslant
v^{2/p}\left(2e^{\frac{2}{\beta}+2}+1\right)\sum_{i\leqslant
t}\lambda_{i}^{2}+\frac{\beta}{2}+\log(\frac{1}{\alpha}),$

Noting that

$\sup_{\boldsymbol{\vartheta}\in\mathbb{S}^{d-1}}\sum_{i\leqslant
t}\lambda_{i}\langle\boldsymbol{\vartheta},\mathsf{th}_{i}(\boldsymbol{X}_{i})-\boldsymbol{\mu}\rangle=\norm{\sum_{i\leqslant
t}\lambda_{i}(\mathsf{th}_{i}(\boldsymbol{X}_{i})-\boldsymbol{\mu})},$

and dividing through by $\sum_{i\leqslant t}\lambda_{i}$ gives the desired
result.

###  D.3 Section 4

#### Proof of Lemma 5

Recall that from the proof of Lemma 1 that $\mathbb{E}_{P}[\exp
F_{t}(\boldsymbol{\theta})|\mathcal{F}_{t-1}]\leqslant 1$. Further, as
discussed in the main paper, we have

$F_{t}(\boldsymbol{\theta}):=\frac{\lambda_{t}}{2B}Z_{t}(\boldsymbol{\theta})-\frac{\psi_{E}(\lambda_{t})}{(2B)^{2}}\langle\boldsymbol{\theta},\boldsymbol{X}_{i}-\bar{\boldsymbol{\mu}}_{i-1}\rangle^{2}\in[-1,2],$

and consequently $\exp F_{t}(\boldsymbol{\theta})\in[1/e,e^{2}]$. Therefore,
applying (14) with the function

$f_{t}(\boldsymbol{x})=\exp\left\\{\frac{\lambda_{t}}{2B}\langle\boldsymbol{\theta},\boldsymbol{x}-\boldsymbol{\mu}\rangle-\frac{\psi_{E}(\lambda_{t})}{(2B)^{2}}\langle\boldsymbol{\theta},\boldsymbol{x}-\bar{\boldsymbol{\mu}}_{t-1}\rangle^{2}\right\\},$

(note the dependence on $t$) which is bounded in $[1/e,e^{2}]$, implies that

$\displaystyle\mathbb{E}_{Q}\exp F_{t}(\boldsymbol{\theta})$
$\displaystyle=\mathbb{E}_{\boldsymbol{X}\sim
Q}f_{t}(\boldsymbol{X})\leqslant\mathbb{E}_{\boldsymbol{X}\sim
P}f_{t}(\boldsymbol{X})+\epsilon(e^{2}-1/e)=\mathbb{E}_{P}\exp
F_{t}(\boldsymbol{\theta})+\epsilon(e^{2}-1/e).$

Take any $Q\in\mathscr{B}_{\operatorname{TV}}(P,\epsilon)$. Then,

$\displaystyle\mathbb{E}_{Q}\left[\frac{F_{t}(\boldsymbol{\theta})}{1+e^{2}\epsilon}\bigg{|}\mathcal{F}_{t-1}\right]$
$\displaystyle=\frac{\mathbb{E}_{Q}[\exp
f(F_{t}(\boldsymbol{\theta}))|\mathcal{F}_{t-1}]}{1+e^{2}\epsilon}$
$\displaystyle\leqslant\frac{\mathbb{E}_{P}[f(\exp
F_{t}(\boldsymbol{\theta})|\mathcal{F}_{t-1}]+\epsilon(e^{2}-1/e)}{1+3\epsilon}$
$\displaystyle\leqslant\frac{1+e^{2}\epsilon}{1+e^{2}\epsilon}=1,$

which suffices to demonstrate that
$N_{t}^{R}(\boldsymbol{\theta})=\prod_{i\leqslant
t}\exp\\{F_{i}(\boldsymbol{\theta})\\}(1+e^{2}\epsilon)^{-1}$ is a $Q$-NSM.

#### Proof of Theorem 4

Apply Proposition 1 to the $Q$-NSM

$N_{t}^{R}(\boldsymbol{\theta})=\prod_{i\leqslant
t}\frac{\exp\\{F_{i}(\boldsymbol{\theta})\\}}{1+e^{2}\epsilon},$

we obtain that with probability $1-\alpha$, simultaneously for all $t\geqslant
1$,

$\int_{\mathbb{S}^{d-1}}\sum_{i\leqslant
t}F_{i}(\boldsymbol{\theta})\gamma_{\boldsymbol{\vartheta}}(\mathrm{d}\boldsymbol{\theta})\leqslant
t\log(1+e^{2}\epsilon)+2\kappa A_{d}(\kappa)+\log\frac{1}{\alpha}.$

From here, after expanding $F_{i}(\boldsymbol{\theta})$, the arithmetic is
nearly identical to the proof of Theorem 1, but there is an extra factor of
$t\log(1+e^{2}\epsilon)$ on the RHS. Therefore, following the steps of that
proof we arrange the above the display to obtain

$\norm{\frac{\sum_{i\leqslant
t}\lambda_{i}\boldsymbol{X}_{i}}{\sum_{i\leqslant
t}\lambda_{i}}-\boldsymbol{\mu}}\leqslant\frac{\frac{1}{2B}\sum_{i\leqslant
t}\psi_{E}(\lambda_{i})\overline{Z}_{i}^{2}+2Bt\log(1+e^{2}\epsilon)+4B\kappa
A_{d}(\kappa)+2B\log(1/\alpha)}{A_{d}(\kappa)\sum_{i\leqslant t}\lambda_{i}}.$

From here we take $\kappa=\sqrt{d}$ and apply Lemma 9 to obtain

$\norm{\frac{\sum_{i\leqslant
t}\lambda_{i}\boldsymbol{X}_{i}}{\sum_{i\leqslant
t}\lambda_{i}}-\boldsymbol{\mu}}\leqslant\frac{\frac{1}{2B}\sum_{i\leqslant
t}\psi_{E}(\lambda_{i})\overline{Z}_{i}^{2}+2Bt\log(1+e^{2}\epsilon)+4B+2B\log(1/\alpha)}{\frac{2}{3\sqrt{d}}\sum_{i\leqslant
t}\lambda_{i}}.$ (25)

The final bound then follows from using that $4B\leqslant 4B\log(1/\alpha)$.

##  Appendix E Von Mises-Fisher distribution

Figure 3: Behavior of the constant in the mean of the vMF distribution,
$A_{d}(\kappa)$.

With mean direction $\boldsymbol{\vartheta}\in\mathbb{S}^{d-1}$ and
concentration parameter $\kappa>0$, the density of the $(d-1)$-dimensional von
Mises-Fisher (vMF) distribution $\gamma(\boldsymbol{\vartheta},\kappa)$ (with
respect to the uniform distribution) on $\mathbb{S}^{d-1}$ is

$\gamma(\boldsymbol{x};\boldsymbol{\vartheta},\kappa)=\frac{\kappa^{d/2-1}}{(2\pi)^{d/2}I_{d/2-1}(\kappa)}\exp(\kappa\langle\boldsymbol{\vartheta},\boldsymbol{x}\rangle).$

Its mean is $A_{d}(\kappa)\boldsymbol{\vartheta}\in\mathbb{R}^{d}$, where

$A_{d}(\kappa)=\frac{I_{d/2}(\kappa)}{I_{d/2-1}(\kappa)},$ (26)

and $I_{\alpha}$ is the modified Bessel function of the first kind and order
$\alpha$,

$I_{\alpha}(x)=\sum_{m=0}^{\infty}\frac{(x/2)^{2m+\alpha}}{m!\Gamma(m+\alpha+1)}.$
(27)

######  Lemma 9.

For any $d\geqslant 1$, $\frac{2}{3}<\sqrt{d}A_{d}(\sqrt{d})<1$.

###### Proof.

Note that,

$\displaystyle\sqrt{d}A_{d}(\sqrt{d})=\frac{\sum_{m=0}^{\infty}\sqrt{d}\frac{(\sqrt{d}/2)^{2m+d/2}}{m!\Gamma(m+d/2+1)}}{\sum_{m=0}^{\infty}\frac{(\sqrt{d}/2)^{2m+d/2-1}}{m!\Gamma(m+d/2)}}.$

Denote the $m$th summands of the numerator and denominator by

$W(m,d)=\sum_{m=0}^{\infty}\sqrt{d}\frac{(\sqrt{d}/2)^{2m+d/2}}{m!\Gamma(m+d/2+1)},\quad
V(m,d)=\frac{(\sqrt{d}/2)^{2m+d/2-1}}{m!\Gamma(m+d/2)}.$

Then,

$\frac{W(m,d)}{V(m,d)}=\frac{d/2}{m+d/2}\leqslant 1,$

equality only when $m=0$. Therefore,

$\sqrt{d}A_{d}(\sqrt{d})=\frac{\sum_{m=0}^{\infty}V(m,d)\frac{d/2}{m+d/2}}{\sum_{m=0}^{\infty}V(m,d)}<1.$

Further, observe that

$\frac{W(m,d)}{V(m+1,d)}={2(m+1)}.$

Therefore

$\sqrt{d}A_{d}(\sqrt{d})=\frac{\sum_{m=1}^{\infty}2mV(m,d)}{V(0,d)+\sum_{m=1}^{\infty}V(m,d)},$

and

$\displaystyle\frac{1}{\sqrt{d}A_{d}(\sqrt{d})}$
$\displaystyle=\frac{V(0,d)+\sum_{m=1}^{\infty}V(m,d)}{\sum_{m=1}^{\infty}2mV(m,d)}<\frac{V(0,d)+\sum_{m=1}^{\infty}V(m,d)}{\sum_{m=1}^{\infty}2V(m,d)}$
$\displaystyle<\frac{V(0,d)}{2V(1,d)}+\frac{1}{2}=1+\frac{1}{2}=\frac{3}{2}.$

Thus we conclude that $\frac{2}{3}<\sqrt{d}A_{d}(\sqrt{d})<1$ for all
$d\geqslant 1$. ∎

######  Lemma 10.

The Kullback-Leibler divergence from
$\gamma(\boldsymbol{\vartheta}_{0},\kappa)$ to
$\gamma(\boldsymbol{\vartheta}_{1},\kappa)$ satisfies

$D_{\operatorname{KL}}(\gamma(\boldsymbol{\vartheta}_{1},\kappa)\|\gamma(\boldsymbol{\vartheta}_{0},\kappa))\leqslant
2\kappa A_{d}(\kappa).$

###### Proof.

Let $\boldsymbol{X}\sim\gamma(\boldsymbol{\vartheta}_{1},\kappa)$. By a direct
calculation,

$\displaystyle
D_{\operatorname{KL}}(\gamma(\boldsymbol{\vartheta}_{1},\kappa)\|\gamma(\boldsymbol{\vartheta}_{0},\kappa))=\mathbb{E}\left[\log\frac{\gamma(\boldsymbol{X};\boldsymbol{\vartheta}_{1},\kappa)}{\gamma(\boldsymbol{X};\boldsymbol{\vartheta}_{0},\kappa)}\right]=\mathbb{E}\left[\kappa\langle\boldsymbol{\vartheta}_{1}-\boldsymbol{\vartheta}_{0},\boldsymbol{X}\rangle\right]=\kappa\langle\boldsymbol{\vartheta}_{1}-\boldsymbol{\vartheta}_{0},A_{d}(\kappa)\boldsymbol{\vartheta}_{1}\rangle.$

Since both $\boldsymbol{\vartheta}_{0},\boldsymbol{\vartheta}_{1}$ are on the
unit sphere $\mathbb{S}^{d-1}$, the inner product
$\langle\boldsymbol{\vartheta}_{1}-\boldsymbol{\vartheta}_{0},\boldsymbol{\vartheta}_{1}\rangle$
is upper bounded by 2, which concludes the proof. ∎

##  Appendix F Asymptotic width of empirical-Bernstein bounds

###  F.1 Non-robust case

Here we assume the data are iid and we set $B=1/2$ for theoretical (and
notational) convenience. Let
$\sigma^{2}=\mathbb{V}(\boldsymbol{X})=\mathbb{E}[\norm{\boldsymbol{X}-\boldsymbol{\mu}}^{2}]$.
Fix a sample size $n$, let
$\widehat{\sigma}_{n}^{2}=\frac{1}{n}\sum_{i=1}^{n}\norm{\boldsymbol{X}_{i}-\bar{\boldsymbol{\mu}}_{i-1}}^{2}$,
and consider

$\lambda_{t}=\sqrt{\frac{c\log(1/\alpha)}{\widehat{\sigma}_{t-1}^{2}n}},\quad\forall
t\geqslant 1,$

for some constant $c$. Here we define $\widehat{\sigma}_{0}^{2}=1$. Our goal
is to show that there exists a $c$ such that $\sqrt{n}W_{n}$, where $W_{n}$ is
the width of the empirical-Bernstein interval in Theorem 1, converges almost
surely to a quantity upper bounded by
$\zeta_{\alpha}\sigma\sqrt{d\log(1/\alpha)}$, where $\zeta_{\alpha}$ is a
constant (dependent on $\alpha$) to be determined later. Our proof follows
similar steps to Waudby-Smith and Ramdas (2023, Appendix E.2) who prove the
result in the scalar setting. Most of the relevant mechanics still go through,
however.

We build up to the result via a sequence of lemmas.

######  Lemma 11.

$\widehat{\sigma}_{n}^{2}$ converges to $\sigma^{2}$ almost surely.

###### Proof.

Decompose $\widehat{\sigma}_{n}^{2}$ as follows:

$\displaystyle\widehat{\sigma}_{n}^{2}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\norm{\boldsymbol{X}_{i}-\bar{\boldsymbol{\mu}}_{i-1}}^{2}=\frac{1}{n}\sum_{i=1}^{n}\norm{\boldsymbol{X}_{i}-\boldsymbol{\mu}+\boldsymbol{\mu}-\bar{\boldsymbol{\mu}}_{i-1}}^{2}$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\norm{\boldsymbol{X}_{i}-\boldsymbol{\mu}}^{2}+\frac{2}{n}\sum_{i=1}^{n}\langle\boldsymbol{X}_{i}-\boldsymbol{\mu},\boldsymbol{\mu}-\bar{\boldsymbol{\mu}}_{i-1}\rangle+\frac{1}{n}\sum_{i=1}^{n}\norm{\boldsymbol{\mu}-\bar{\boldsymbol{\mu}}_{i-1}}^{2}.$
(28)

Now, by the SLLN, the first sum converges to $\sigma^{2}$ almost surely. As
for the second sum, write

$\displaystyle\left|\frac{2}{n}\sum_{i=1}^{n}\langle\boldsymbol{X}_{i}-\boldsymbol{\mu},\boldsymbol{\mu}-\bar{\boldsymbol{\mu}}_{i-1}\rangle\right|\leqslant\frac{2}{n}\sum_{i=1}^{n}\norm{\boldsymbol{X}_{i}-\boldsymbol{\mu}}\norm{\boldsymbol{\mu}-\bar{\boldsymbol{\mu}}_{i-1}}\leqslant\frac{2}{n}\sum_{i=1}^{n}\norm{\boldsymbol{\mu}-\bar{\boldsymbol{\mu}}_{i-1}}\xrightarrow{a.s.}0,$

where almost sure convergence follows from combining the following three
observations: (i) If the absolute value of a sequence converges to zero, then
the sequence converges to zero; (ii) If a sequence converges to a value, then
its partial sums converge to that value; (iii) $\bar{\boldsymbol{\mu}}_{i-1}$
converges to $\boldsymbol{\mu}$ a.s., so by continuous mapping theorem
$\norm{\boldsymbol{\mu}-\bar{\boldsymbol{\mu}}_{i-1}}\xrightarrow{a.s}0$. The
third sum in (28) converges almost surely to 0 for the same reason. This
completes the proof. ∎

######  Lemma 12.

$\sum_{i\leqslant
n}\psi_{E}(\lambda_{i})\norm{\boldsymbol{X}_{i}-\bar{\boldsymbol{\mu}}_{i-1}}^{2}$
converges to $\frac{c}{2}\log(1/\alpha)$ almost surely.

###### Proof.

Define the function $\psi_{H}(\lambda_{i})=\lambda_{i}^{2}/8$. The “H” stands
for Hoeffding, and the notation is borrowed from Howard et al. (2020). Now,
Waudby-Smith and Ramdas (2023) demonstrate that
$\frac{\psi_{E}(\lambda)}{4\psi_{H}(\lambda_{i})}\to 1$ as $\lambda\to 0$.
Therefore, we can write $\psi_{E}(\lambda_{i})/\psi_{H}(\lambda_{i})=4+4v_{i}$
where $v_{i}\to 0$ as $i\to\infty$, since $\lambda_{i}\xrightarrow{a.s.}0$.
Furthermore, by Lemma 11, we can write
$\sigma^{2}/\widehat{\sigma}_{n}^{2}=1+u_{n}$ for some
$u_{n}\xrightarrow{a.s.}0$. Therefore,

$\displaystyle\sum_{i=1}^{n}\psi_{E}(\lambda_{i})\norm{\boldsymbol{X}_{i}-\bar{\boldsymbol{\mu}}_{i-1}}^{2}$
$\displaystyle=\sum_{i=1}^{n}\frac{\psi_{E}(\lambda_{i})}{\psi_{H}(\lambda_{i})}\psi_{H}(\lambda_{i})\norm{\boldsymbol{X}_{i}-\bar{\boldsymbol{\mu}}_{i-1}}^{2}$
$\displaystyle=\frac{1}{2}\sum_{i=1}^{n}\lambda_{i}^{2}\norm{\boldsymbol{X}_{i}-\bar{\boldsymbol{\mu}}_{i-1}}^{2}(1+v_{i})$
$\displaystyle=\frac{c\log(1/\alpha)}{2n}\sum_{i=1}^{n}\frac{\norm{\boldsymbol{X}_{i}-\bar{\boldsymbol{\mu}}_{i-1}}^{2}}{\widehat{\sigma}_{i}^{2}}(1+v_{i})$
$\displaystyle=\underbrace{\frac{c\log(1/\alpha)}{2n}\sum_{i=1}^{n}\frac{\norm{\boldsymbol{X}_{i}-\bar{\boldsymbol{\mu}}_{i-1}}^{2}}{\sigma^{2}}}_{:=K_{i}}(1+u_{i})(1+v_{i}).$

Write $(1+u_{i})(1+v_{i})=1+w_{i}$ where $w_{i}\xrightarrow{a.s.}0$. Then, the
above sum becomes

$\displaystyle\sum_{i=1}^{n}K_{i}(1+w_{i})$
$\displaystyle=\frac{c\log(1/\alpha)}{2\sigma^{2}}\left\\{\frac{1}{n}\sum_{i=1}^{n}\norm{\boldsymbol{X}_{i}-\bar{\boldsymbol{\mu}}_{i-1}}^{2}(1+w_{i})\right\\}$
$\displaystyle=\frac{c\log(1/\alpha)}{2\sigma^{2}}\left\\{\widehat{\sigma}_{n}^{2}+\widehat{\sigma}_{n}^{2}w_{i}\right\\}\xrightarrow{a.s.}\frac{c\log(1/\alpha)}{2},$

since $\widehat{\sigma}_{n}^{2}\xrightarrow{a.s.}\sigma^{2}$ by Lemma 11 and
$w_{i}\xrightarrow{a.s.}0$.

∎

######  Lemma 13.

$\frac{1}{\sqrt{n}}\sum_{i\leqslant n}\lambda_{i}$ converges to
$\sqrt{\frac{c\log(1/\alpha)}{\sigma^{2}}}$ almost surely.

###### Proof.

Via similar reasoning as above, we can write $\sigma/\widehat{\sigma}_{i}\to
1+u_{i}$ for some $u_{i}\xrightarrow{a.s.}0$. (Here we’ve employed the
continuous mapping theorem on Lemma 11.) Therefore,

$\displaystyle\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\lambda_{i}$
$\displaystyle=\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\sqrt{\frac{c\log(1/\alpha)}{\widehat{\sigma}_{i-1}^{2}n}}=\frac{\sqrt{c\log(1/\alpha)}}{n\sigma}\sum_{i=1}^{n}\sqrt{\frac{\sigma^{2}}{\widehat{\sigma}_{i-1}^{2}}}\xrightarrow{a.s.}\sqrt{\frac{c\log(1/\alpha)}{\sigma^{2}}},$

where we used the fact that if a sequence converges to a value then so too do
its partial sums. ∎

Now, let

$W_{n}:=\frac{\sum_{i\leqslant
t}\psi_{E}(\lambda_{i})\norm{\boldsymbol{X}_{i}-\bar{\boldsymbol{\mu}}_{i-1}}^{2}}{A_{d}(\kappa)\sum_{i\leqslant
t}\lambda_{i}}+\frac{2\kappa
A_{d}(\kappa)+\log\frac{1}{\alpha}}{A_{d}(\kappa)\sum_{i\leqslant
t}\lambda_{i}},$ (29)

be the width of our CI in Theorem 1. We have

$\displaystyle\sqrt{n}W_{n}$ $\displaystyle=\frac{\sum_{i\leqslant
t}\psi_{E}(\lambda_{i})\norm{\boldsymbol{X}_{i}-\bar{\boldsymbol{\mu}}_{i-1}}^{2}+2\kappa
A_{d}(\kappa)+\log\frac{1}{\alpha}}{\frac{1}{\sqrt{n}}A_{d}(\kappa)\sum_{i\leqslant
t}\lambda_{i}}$
$\displaystyle\xrightarrow{a.s.}\frac{\frac{c}{2}\log(1/\alpha)+2\kappa
A_{d}(\kappa)+\log(1/\alpha)}{A_{d}(\kappa)\sqrt{c\log(1/\alpha)/\sigma^{2}}}$
$\displaystyle=\sigma\left(\frac{\sqrt{c\log(1/\alpha)}}{2A_{d}(\kappa)}+\frac{2\kappa}{\sqrt{c\log(1/\alpha)}}+\frac{\sqrt{\log(1/\alpha)}}{A_{d}(\kappa)\sqrt{c}}\right)=:G(d,\alpha,c).$

As suggested by Lemma 9, consider setting $\kappa=\sqrt{d}$. Then
$A_{d}(\kappa)\geqslant\frac{2}{3\sqrt{d}}$ so after some rearranging we
obtain that

$G(d,\alpha,c)\leqslant\sigma\sqrt{d\log(1/\alpha)}\left(\frac{3\sqrt{c}}{4}+\frac{2}{\sqrt{c}\log(1/\alpha)}+\frac{3}{2\sqrt{c}}\right).$
(30)

Some elementary calculus shows that the value of $c$ which minimizes the above
right hand side is

$c^{*}=2+\frac{8}{3\log(1/\alpha)}.$

However, even with such a choice, the multiplier on
$\sigma\sqrt{d\log(1/\alpha)}$ in (30) is unbounded as $\alpha\to 1$. However,
for reasonable choices of $\alpha$, such as $\alpha\leqslant 0.1$, we have

$\zeta_{\alpha}:=\frac{3\sqrt{c^{*}}}{4}+\frac{2}{\sqrt{c^{*}}\log(1/\alpha)}+\frac{3}{2\sqrt{c^{*}}}=\sqrt{\frac{9}{2}+\frac{6}{\log(1/\alpha)}}\leqslant\begin{cases}2.67,&\alpha\leqslant
0.1\\\ 2.55,&\alpha\leqslant 0.05\\\ 2.41,&\alpha\leqslant 0.01.\end{cases}$

Thus, we obtain an asymptotic width bounded by
$\sigma\sqrt{\zeta_{\alpha}d\log(1/\alpha)}$ for all $\alpha\in(0,1)$, and for
$\alpha\leqslant 0.1$ we can simplify the bound to
$3.25\sigma\sqrt{d\log(1/\alpha)}$. We emphasize that $d$ must be considered
fixed in the above analysis. That is, $d$ cannot grow with $n$. Obtaining such
a result is an interesting open problem.

###  F.2 Robust case

Now we turn our attention to the robust empirical-Bernstein bound given in
Section 4. Consider $\lambda_{i}=\lambda$. Assuming that
$\psi_{E}(\lambda)\leqslant\lambda^{2}$, the width of the CSS in Theorem 4 is
bounded by, up to constants,

$g(\lambda):=\frac{\sqrt{d}\lambda\sum_{i}\overline{Z}_{i}^{2}}{Bt}+\frac{B\sqrt{d}\log(1/\alpha)}{t\lambda}+\frac{B\sqrt{d}\epsilon}{\lambda}.$

Therefore, a naive constant $\lambda$ independent of any parameters
$(\lambda=0.5)$ say, gives a width of

$W_{t}\lesssim\frac{\sqrt{d}}{B}\frac{1}{t}\sum_{i}\overline{Z}_{i}^{2}+\frac{B\sqrt{d}\log(1/\alpha)}{t}+B\sqrt{d}\epsilon\xrightarrow{a.s.}\sigma^{2}\frac{\sqrt{d}}{B}+B\sqrt{d}\epsilon.$

However, we can optimize $\lambda$ by differentiating $g$. Optimizing for a
particular time $t=n$,

$\displaystyle
g^{\prime}(\lambda)=\frac{\sqrt{d}\sum_{i}\overline{Z}_{i}^{2}}{Bn}-\frac{B\sqrt{d}\log(1/\alpha)}{n\lambda^{2}}-\frac{B\sqrt{d}\epsilon}{\lambda^{2}}.$

Setting $g^{\prime}(\lambda)=0$ and solving for $\lambda$ gives

$\displaystyle\lambda^{2}$
$\displaystyle=\frac{B^{2}\log(1/\alpha)}{\sum_{i}\overline{Z}_{i}^{2}}+\frac{B^{2}n\epsilon}{\sum_{i}\overline{Z}_{i}^{2}}$
$\displaystyle=\frac{B^{2}\log(1/\alpha)}{n\frac{1}{n}\sum_{i}\overline{Z}_{i}^{2}}+\frac{B^{2}\epsilon}{\sum_{i}\frac{1}{t}\overline{Z}_{i}^{2}}$
$\displaystyle\approx\frac{B^{2}\log(1/\alpha)}{n\sigma^{2}}+\frac{B^{2}\epsilon}{\sigma^{2}},$

so that the optimal $\lambda$ is approximately

$\lambda^{*}:=\frac{B}{\sigma}\sqrt{\frac{\log(1/\alpha)}{n}+\epsilon}.$

Plugging this back into $g(\lambda)$ gives an asymptotic width of
$g(\lambda^{*})\lesssim\sigma\sqrt{d\epsilon}$. In the fixed-time setting,
this suggests taking

$\lambda_{t}=\frac{B}{\widehat{\sigma}_{t-1}}\sqrt{\frac{\log(1/\alpha)}{n}+\epsilon}.$

Now let us consider the sequential setting. We are considering

$\lambda_{t}=\frac{1}{b\widehat{\sigma}_{t-1}},$

for some $b>0$. We will assume that $\lambda_{t}\leqslant 0.68$ so that
$\psi_{E}(\lambda_{i})\leqslant\lambda_{i}^{2}$. Plugging this into the bound
of Theorem 4 gives

$\displaystyle\norm{\frac{\sum_{i\leqslant
t}\lambda_{i}\boldsymbol{X}_{i}}{\sum_{i\leqslant
t}\lambda_{i}}-\boldsymbol{\mu}}$
$\displaystyle\leqslant\frac{\frac{\sqrt{d}}{2B}\sum_{i\leqslant
t}\psi_{E}(\lambda_{i})\overline{Z}_{i}^{2}+t\log(1+e^{2}\epsilon)+6B\sqrt{d}\log(1/\alpha)}{\frac{2}{3}\sum_{i\leqslant
t}\lambda_{i}}$
$\displaystyle\leqslant\underbrace{\frac{\frac{\sqrt{d}}{2B}\sum_{i\leqslant
t}\lambda_{i}^{2}\overline{Z}_{i}^{2}}{\frac{2}{3}\sum_{i\leqslant
t}\lambda_{i}}}_{(i)}+\underbrace{\frac{6B\sqrt{d}\log(1/\alpha)}{\frac{2}{3}\sum_{i\leqslant
t}\lambda_{i}}}_{(ii)}+\underbrace{\frac{3te^{2}\epsilon}{2\sum_{i\leqslant
t}\lambda_{i}}}_{(iii)}.$

Before analyzing these terms, first note that

$\displaystyle\frac{1}{t}\sum_{i\leqslant
t}\lambda_{i}=\frac{1}{bt}\sum_{i\leqslant
t}\frac{1}{\widehat{\sigma}_{i-1}}=\frac{1}{b\sigma t}\sum_{i\leqslant
t}\frac{\sigma}{\widehat{\sigma}_{i-1}}\xrightarrow{a.s.}\frac{1}{b\sigma},$

using the arguments of Section F.1. The convergence of terms (ii) and (iii)
follow as:

$\displaystyle\frac{6B\sqrt{d}\log(1/\alpha)}{\frac{2}{3}\sum_{i\leqslant
t}\lambda_{i}}=\frac{9B\sqrt{d}\log(1/\alpha)}{t}\frac{1}{\frac{1}{t}\sum_{i\leqslant
t}\lambda_{i}}\xrightarrow{a.s.}0,$ (31)

due to the continuous mapping theorem. For term (iii), we have

$\displaystyle\frac{3te^{2}\epsilon}{2\sum_{i\leqslant
t}\lambda_{i}}=\frac{3e^{2}\epsilon}{\frac{2}{t}\sum_{i\leqslant
t}\lambda_{i}}\xrightarrow{a.s.}\frac{3b\sigma e^{2}\epsilon}{2}.$ (32)

As for term (i), notice that

$\displaystyle\frac{1}{t}\sum_{i\leqslant
t}\lambda_{i}^{2}\overline{Z}_{i}^{2}$
$\displaystyle=\frac{1}{b^{2}t}\sum_{i\leqslant
t}\frac{\norm{\boldsymbol{X}_{i}-\bar{\boldsymbol{\mu}}_{i-1}}^{2}}{\widehat{\sigma}_{t-1}^{2}}$
$\displaystyle=\frac{1}{\sigma^{2}b^{2}t}\sum_{i\leqslant
t}\norm{\boldsymbol{X}_{i}-\bar{\boldsymbol{\mu}}_{i-1}}^{2}\frac{\sigma^{2}}{\widehat{\sigma}_{t-1}^{2}}$
$\displaystyle=\frac{1}{\sigma^{2}b^{2}t}\sum_{i\leqslant
t}\norm{\boldsymbol{X}_{i}-\bar{\boldsymbol{\mu}}_{i-1}}^{2}(1+u_{i}),$

where $u_{i}\xrightarrow{a.s.}0$ via the arguments in Section F.1. Moreover,
by Lemma 11,
$\frac{1}{t}\sum_{i}\norm{\boldsymbol{X}_{i}-\bar{\boldsymbol{\mu}}_{i-1}}^{2}\xrightarrow{a.s.}\sigma^{2}$.
Therefore, rewriting term (i),

$\displaystyle\frac{\frac{\sqrt{d}}{2B}\sum_{i\leqslant
t}\lambda_{i}^{2}\overline{Z}_{i}^{2}}{\frac{2}{3}\sum_{i\leqslant
t}\lambda_{i}}$
$\displaystyle=\frac{3\sqrt{d}}{B}\frac{\frac{1}{t}\sum_{i\leqslant
t}\lambda_{i}^{2}\overline{Z}_{i}^{2}}{\frac{1}{t}\sum_{i\leqslant
t}\lambda_{i}}\xrightarrow{a.s.}\frac{3\sqrt{d}}{B}\frac{1/b^{2}}{(b\sigma)^{-1}}=\frac{3\sigma\sqrt{d}B}{b}.$
(33)

Putting together (31), (32) and (33) gives that the width converges almost
surely to a quantity upper bounded by

$\frac{3\sigma\sqrt{d}B}{b}+\frac{3b\sigma^{2}e^{2}\epsilon}{2}\lesssim\sigma\sqrt{Bd}+\sqrt{B}\sigma\epsilon,$

if $b\propto\sqrt{B}$.