# Risk Bounds for Learning via Hilbert Coresets

Spencer Douglas1, Piyush Kumar1, R.K. Prasanth1
1Systems & Technology Research
600 West Cummings Park, Suite 6500, Woburn, MA 01801, USA

###### Abstract

We develop a formalism for constructing stochastic upper bounds on the
expected full sample risk for supervised classification tasks via the Hilbert
coresets approach within a transductive framework. We explicitly compute tight
and meaningful bounds for complex datasets and complex hypothesis classes such
as state-of-the-art deep neural network architectures. The bounds we develop
exhibit nice properties: i) the bounds are non-uniform in the hypothesis space
$\mathcal{H}$, ii) in many practical examples, the bounds become effectively
deterministic by appropriate choice of prior and training data-dependent
posterior distributions on the hypothesis space, and iii) the bounds become
significantly better with increase in the size of the training set. We also
lay out some ideas to explore for future research.

##  I Introduction

Generalization bounds for learning provide a theoretical guarantee on the
performance of a learning algorithm on unseen data. The goal of such bounds is
to provide control of the error on unseen data with pre-specified confidence.
In certain situations, such bounds may also help in designing new learning
algorithms. With the great success of deep neural networks (DNNs) on a variety
of machine learning tasks, it is natural to try to construct theoretical
bounds for such models. Recent years have seen several notable efforts to
construct these bounds.

One classical approach, that of Probably-Approximately-Correct (PAC) learning,
aims to develop generalization bounds that hold (with high probability)
uniformly over all hypotheses $h$ (parametrized with parameters $\theta$) in
some hypothesis class $\mathcal{H}$ and over all data distributions $D$ in
some class of distributions $\mathcal{D}$. The well-known VC-dimension bounds
[1] fall in this category. Since these are uniform bounds, they are _worst
case_ and very general. However, the flip side is that they tend to be very
loose and not useful in practical situations. For example, the VC dimension of
neural networks is typically bounded by the number of parameters of the
network. Since modern DNNs have millions of parameters, these uniform VC
dimension bounds are typically not useful at all.

Better prospects for developing useful and non-trivial bounds arise when one
develops _data-dependent_ and _non-uniform_ bounds. For example, there has
been lot of recent work on bounds based on _Rademacher complexity_ [2, 3].
These bounds are dependent on the particular data distribution $D$
($\in\mathcal{D}$) but still depend uniformly on the entire hypothesis class
$\mathcal{H}$. Recent work on computing these bounds for modern DNNs still
gives rise to bounds that may not be useful for at least some classes of
modern DNNs [4].

Another data-dependent approach to developing useful generalization bounds is
that of the _PAC-Bayesian_ framework [5, 6, 7, 8]. In this case, one starts
with some ’prior’ distribution $\pi_{0}$ (that is independent of the data)
over the hypothesis space $\mathcal{H}$, and then constructs a ’posterior’
distribution $\nu$ over $\mathcal{H}$ after training the model on some
training dataset. Note that the ’prior’ and ’posterior’ do not refer to the
corresponding terms in Bayesian probability in the strict sense. Here, the
posterior distribution $\nu$ can be thought of as being determined by the
training data and the training algorithm in general. The PAC-Bayes framework
studies an averaging of the generalization error according to the posterior
$\nu$ (the so-called Gibbs classifier). The bounds thus obtained are valid for
all posterior distributions $\nu$ over the hypothesis space $\mathcal{H}$, but
the bound itself depends on the posterior $\nu$. Thus, these bounds are _non-
uniform_ in the hypothesis space as well. In other words, the precise
numerical bounds are valid for parameters in the neighborhood of those learned
by the actual learning algorithm. Thus, these have the potential to provide
tight and useful bounds for complex models like DNNs.

A point to note is that the bulk of the work in statistical learning theory
within each of the above approaches has dealt with providing bounds for the
_inductive_ learning setting. In this case, the goal is to learn a labeling
rule from a finite set of labeled training examples, and this rule is then
used to label new (test) examples. In many realistic situations, however, one
has a _full sample_ , which consists of a training set of labeled examples,
together with an unlabeled set of points which needs to be labeled.This is the
_transductive_ setting, where one is interested only in labeling the unlabeled
set as accurately as possible. In this case, it is plausible that if one
leverages the information of the test data points appropriately, one may be
able to derive tighter risk bounds. Transductive learning bounds are valid for
the particular dataset at hand. If the dataset is augmented with new data
points, the bounds have to be computed again in general. In this work, we show
that it is possible to get tight bounds on the full sample risk by using the
method of Hilbert coresets, a coreset defined in terms of a Hilbert norm on
the space of loss functions (a more precise description will follow later).
This work builds on the work of Campbell and Broderick [9, 10] and applies it
to the transductive learning framework. Our principal contributions in this
paper are as follows:

  * •

We compute non-vacuous and tight upper bounds on the _expected full sample
risk_ for multi-class classification tasks for a variety of state-of-the-art
hypothesis classes (DNN architectures) and datasets, and compare these results
with that obtained from explicit training algorithms. Our bounds are valid for
any data distribution that lies in a class of distributions $\mathcal{D}$. In
particular, we study the class of distributions $\mathcal{D}$: _Convex hull of
a probability simplex_ $\mathcal{D}_{CH}^{(K)}$. The calculation of the bounds
for the expected full sample risk requires knowledge of the full dataset along
with the labels.

  * •

The developed bounds have nice properties. First, these bounds are _non-
uniform_ in the sense that although the bounds are valid for all distributions
over $\mathcal{H}$, the precise bound depends on particular distribution over
$\mathcal{H}$, i.e. the precise neighborhood of the parameters
$\theta_{\star}$ learned by the training algorithm.

  * •

Second, although the bounds we develop technically apply to a _stochastic_
DNN, appropriate choice of data-independent priors and data-dependent
posteriors in many practical examples gives rise to Gibbs classifiers that are
highly clustered around the classifier learned by the original deterministic
DNN (with parameters $\theta_{\star}$ learned by the training algorithm, see
results in Section IV-E). Thus, although there is no guarantee, the bounds can
become effectively _deterministic_ in many practical cases.

  * •

Finally, the numerical bounds _decrease_ with an increase in the number of
training samples. This behavior has been tested for a variety of DNNs and
holds up across changes in depth and width of the DNN as well as different
mini-batch size used for training. This is in contrast with certain existing
bounds which may even _increase_ with an increase in the number of training
samples [3, 11, 12]; see [13] for a discussion.

Transductive bounds have received much less attention than inductive bounds in
general. Generalization bounds for the transductive setting have been
developed within the PAC-Bayesian framework in [14, 15]. Our bounds are not
PAC-Bayes bounds in the strict sense. Our methods do provide a bound on the
expected full sample risk as in the PAC-Bayes case, however this bound is
computed using Hilbert coresets which is a novel approach. Also, the bounds
developed in [14, 15] are applicable for a uniform distribution on the full
sample and have only been applied to simple hypothesis classes. On the other
hand, the bounds developed in this paper are applicable to more general data
distributions, such as those that lie in the convex hull of the probability
simplex of $K$ probability vectors. Finally, to our knowledge our paper is the
first one to develop meaningful and tight upper bounds within a transductive
framework for state-of-the-art DNNs.

##  II Problem Formulation

We start with considering an arbitrary input space $\mathcal{X}$ and binary
output space $\mathcal{Y}=\\{-1,+1\\}$. For simplicity, we focus on the binary
classification case; however our results are easily generalizable to the
multi-class classification case as well. Consistent with the original
transductive setting of Vapnik [16], we consider a _full sample_ $Z$ of $N$
examples (input-output pairs) that is assumed to be drawn i.i.d from an
unknown distribution $D$ on $\mathcal{X}\times\mathcal{Y}$:

$Z:=\\{(x_{n},y_{n})\\}_{n=1}^{N}$ (1)

For simplicity we will refer to the probability distribution as
$D\equiv\\{p_{n}\\}=\\{P(x_{n},y_{n})\\};n=1,2,..,N$. $D$ is assumed to lie in
a class of data distributions $\mathcal{D}$. As mentioned previously, we
consider the class of distributions where $\mathcal{D}$ is the _convex hull of
a probability simplex_ $\mathcal{D}_{CH}^{(K)}(\\{\hat{p}^{(i)}\\})_{i\in[K]}$
defined by the $K$ probability vectors $(\\{\hat{p}^{(i)}\\});i\in[K]$. A
training dataset
$Z^{(S)}=\\{(x_{i},y_{i})_{i=1}^{N_{S}}\\}\in\mathcal{X}\times\mathcal{Y}$,
that contains $N_{S}$ input-output pairs is obtained from the full sample by
sampling without replacement. Given $Z^{(S)}$ and the unlabeled set $U_{S}$
consisting of $N-N_{S}$ inputs, the goal in transductive learning is to learn
a classifier that correctly classifies the unlabeled inputs of the set
$U_{S}$, i.e. minimize the risk of the classifier on $U_{S}$.

We consider in this work general loss functions (including unbounded loss
functions) $l:\mathcal{Y}\times\mathcal{Y}\rightarrow\mathbf{R}$, and a
(discrete or continuous) set $\mathcal{H}$ of hypotheses
$h:\Theta\times\mathcal{X}\rightarrow\mathcal{Y}$ where $\theta$ lies in some
parameter space $\Theta$. Although the bounds we develop apply to general loss
functions, most of the relevant applications of our bounds for classification
apply to bounded, and in particular _zero-one_ , loss functions. We provide a
bound on the _full sample risk_ of a hypothesis $h\in\mathcal{H}$, which is
given by:

$R(h;\theta)=\mathbf{E}_{Z\sim D}\;l(h(\theta,x),y)$ (2)

Inspired by the PAC-Bayesian approach, we study an averaging of the full
sample risk according to a posterior distribution $\nu$ over $\mathcal{H}$
that is dependent on the data and the training algorithm in general. Since
$\mathcal{H}$ is parametrized by $\theta\in\Theta$, by a slight abuse of
notation we denote the distribution on the parameters $\theta$ by
$\nu(\theta)$ as well. So, we aim to provide bounds on the _expected full
sample risk_ $R(\nu)$:

$R(\nu)=\mathbf{E}_{\theta\sim\nu}\,R(h;\theta)=\mathbf{E}_{\theta\sim\nu}\;\mathbf{E}_{Z\sim
D}\;l(h(\theta,x),y)=\mathbf{E}_{\theta\sim\nu}\left(\sum_{n=1}^{N}\,p_{n}\,l_{n}(\theta)\right):=\mathbf{E}_{\theta\sim\nu}\,L(\theta)$
(3)

where the individual loss for a data point is defined as
$l_{n}(\theta):=l(h(\theta,x_{n}),y_{n})$ and
$L(\theta):=\sum_{n=1}^{N}\,p_{n}\,l_{n}(\theta)$ is the total loss function
of the full sample. Both $L(\theta)$ and $l_{n}(\theta),n\in[N]$ are functions
of the parameters $\theta$. We aim to develop upper bounds on $R(\nu)$ using
the Hilbert coresets approach, which we summarize next.

###  II-A Summary of Hilbert Coresets Approach

As described above, we consider a dataset $Z=\\{(x_{n},y_{n})\\}_{n=1}^{N}\sim
D$, parameters $\theta\in\Theta$, and a _cost function_
$L:\Theta\rightarrow\mathbf{R}$ that is additively decomposable into a set of
functions $L_{n}(\theta)$, i.e. $L(\theta):=\sum_{n=1}^{N}\,L_{n}(\theta)$. By
comparing to (3) and the discussion following that, we can make the
identification:

$L_{n}(\theta):=p_{n}\,l_{n}(\theta)$ (4)

Now, an _$\epsilon$ -coreset_ is given by a weighted dataset with nonnegative
weights $w_{n},n=1,2...,N$, only a small number $m$ of which are non-zero,
such that the weighted cost function
$L(w,\theta):=\sum_{n=1}^{N}\,w_{n}\,L_{n}(\theta)$ satisfies:

$|L(w,\theta)-L(\theta)|\leq\epsilon\,L(\theta)\hskip
14.22636pt\forall\theta\in\Theta$ (5)

for $\epsilon>0$. If we consider a bounded uniform norm for functions of the
form $f:\Theta\rightarrow\mathbf{R}$, which is weighted by the total cost
$L(\theta)$ as in

$||f||:=\sup_{\theta\in\Theta}\;\left|\frac{f(\theta)}{L(\theta)}\right|,$ (6)

we can rewrite the problem of finding the best coreset of size $m$ as:

$\min_{w}\;||L(w)-L||^{2}\hskip
8.5359pt\mathrm{s.t.}\sum_{n=1}^{N}\,\mathbf{I}\,[w_{n}>0]\leq m;\;\;w\geq 0$
(7)

The norm $||f||$ of any function $f:\Theta\rightarrow\mathbf{R}$ gets rid of
its dependence on $\theta$, thus there is no dependence on $\theta$ in the
norm in (7). The key insight of Campbell _et al_ [9, 10] was to approach the
problem from a geometrical point of view. In particular, the idea is to
construct coresets for cost functions in a _Hilbert space_ , i.e. in which the
cost functions $L_{n}(\theta)$ (defined in the beginning of this subsection)
are endowed with a norm corresponding to an inner product. The inner product
provides a notion of directionality to the various $L_{n}$ and allows us to
use more general norms than the uniform norm in (6) (which we will see in
section II-B). This in turn provides us with two distinct benefits:

  * •

It allows one to choose coreset points intelligently based on the residual
cost approximation error.

  * •

Theoretical guarantees on approximation quality incorporate the error in the
approximation via the (mis)alignment of $L_{n}$ relative to the total cost
$L$.

###  II-B Hilbert Coresets via Frank-Wolfe

Campbell et al [9, 10] developed various algorithms for computing Hilbert
coresets; here we focus on the Frank-Wolfe algorithm since it produces
exponential convergence for moderate values of the coreset size $m$.
Furthermore, it can be elegantly tied to risk bounds and computed explicitly
in many interesting cases. In this case, the problem in (7) is reformulated by
replacing the cardinality constraint on the weights with a polytope constraint
on $w$:

$\displaystyle\min_{w}\;||L(w)-L||^{2}$ $\displaystyle s.t.$
$\displaystyle\sum_{n=1}^{N}w_{n}||L_{n}||=\sum_{n=1}^{N}||L_{n}||;\;\;w\geq
0$

One can rewrite the objective function $||L(w)-L||^{2}$ as
$J(w):=(w-1)^{T}\,K\,(w-1)$ with the matrix $\hat{K}$ having elements defined
in terms of the inner product, $\hat{K}_{n,m}:=\langle L_{n},L_{m}\rangle$
(which will be defined shortly). This optimization problem is then solved by
applying the Frank-Wolfe algorithm for $m$ iterations to find the coresets
$w$[9]. The size of the coreset (number of non-zero coreset weights) is
smaller than $m$ (typically much smaller), see details in Algorithm 2 of [9].
Finally, it can be shown that the Frank-Wolfe algorithm satisfies:

$||L(w)-L||\leq\frac{\sigma\,\eta\,\bar{\eta}\,\beta}{\sqrt{\eta^{2}\,(m-1)+\bar{\eta}^{2}\,\beta^{-2(m-2)}}},$
(9)

where $\sigma:=\sum_{n=1}^{N}\,||L_{n}||$. The remaining quantities
$\\{\eta,\,\bar{\eta},\beta\\}$ can be computed from all the data
($\\{x_{n},y_{n}\\};n=1,2...,N$) and the $\nu$-weighted $L2$-norm, see [9].
Variants of these quantities (which will be relevant for us later) will be
defined precisely in section III. In particular $0<\beta<1$, and for
moderately large $m$, the second term in the denominator dominates giving rise
to an exponential dependence on $m$ ($\sim\beta^{m}$) .

The above bound on the coreset approximation error for the cost $L$ is valid
for any choice of norm. However, for our purposes we will be concerned with
the $\nu$-weighted $L2$-norm defined by:

$||L(w)-L||^{2}_{\nu,2}:=\mathbf{E}_{\theta\sim\nu}\,[L(w,\theta)-L(\theta)]^{2},$
(10)

with an induced inner product:

$\langle
L_{n},L_{m}\rangle:=\mathbf{E}_{\theta\sim\nu}[L_{n}(\theta)\,L_{m}(\theta)]=K_{n,m},$
(11)

where $K_{n,m}$ was defined below (II-B). In the following, we will see that
an upper bound on a variant of $||L-L(w)||$ can be connected to an upper bound
on the full sample risk.

##  III Risk Bounds via Frank-Wolfe Hilbert Coresets

Leveraging the Frank-Wolfe Hilbert coreset construction scheme, we provide an
upper bound on the expected full sample risk which is valid for all posteriors
$\nu\in\mathcal{H}$. However, the precise bound is dependent on the posterior
$\nu$. Furthermore, the bound is independent of the data-distribution over the
set of data-distributions $\mathcal{D}$. The result can be stated as follows:

###### Theorem III.1

Given a full sample $Z$ of size $N$, an unknown distribution on the data
belonging to $\mathcal{D}^{(K)}_{CH}(\\{\hat{p}^{(i)}\\})_{i\in[K]}$, a
hypothesis space $\mathcal{H}:\Theta\times\mathcal{X}\rightarrow\mathcal{Y}$,
a loss function $l:\mathcal{Y}\times\mathcal{Y}\rightarrow\mathbf{R}$, for any
prior distribution $\pi_{0}$ on $\mathcal{H}$, then the output
$\\{\tilde{w},p\\}$ obtained after $m$ iterations of applying _Algorithm 1_
satisfies:

$\displaystyle\forall\nu\,\mathrm{on}\,\mathcal{H},\,\forall
D\in\mathcal{D}_{CH^{(K)}},\hskip 14.22636pt|R(\nu)|$ $\displaystyle\leq$
$\displaystyle||L(\tilde{w})||_{\nu,2}+||L(\tilde{w})-L||_{\nu,2}$ (12)
$\displaystyle\leq$
$\displaystyle||L(\tilde{w})||_{\nu,2}+\frac{\hat{\sigma}\,\hat{\eta}\,\hat{\bar{\eta}}\,\hat{\beta}}{\sqrt{\hat{\bar{\eta}}^{2}\,\hat{\beta}^{-2(m-1)}+\hat{\eta}^{2}\,(m-1)}},$
(13)

where $L(\theta)$ is defined as in (3),
$L(\tilde{w},\theta):=\sum_{n=1}^{N}\,\tilde{w}_{n}\,l_{n}(\theta);\tilde{w}_{n}:=p_{n}\,w_{n}$,
and the quantities
$\\{\hat{\sigma},\hat{\eta},\,\hat{\bar{\eta}},\hat{\beta}\\}$ are defined as
follows:

$\displaystyle\hat{\sigma}$ $\displaystyle:=$
$\displaystyle\max_{i\in[K]}\,\sigma^{(i)};\;\hat{\eta}:=\max_{i\in[K]}\,\eta^{(i)};\;\hat{\beta}:=\max_{i\in[K]}\,\beta^{(i)}$
$\displaystyle\sigma^{(i)}$ $\displaystyle:=$
$\displaystyle\sum_{n=1}^{N}\,\hat{p}^{(i)}_{n}\,||l_{n}||;\;{\eta^{(i)}}^{2}=1-\frac{||L^{(i)}||^{2}}{(\sigma^{(i)})^{2}};\;{\beta^{(i)}}^{2}=1-\frac{(r^{(i)})^{2}}{(\sigma^{(i)})^{2}\hat{\bar{\eta}}^{2}};\;\;i\in[K]$
$\displaystyle\hat{\bar{\eta}}^{2}$ $\displaystyle:=$
$\displaystyle\max_{n,m}\left\lVert\frac{l_{n}}{||l_{n}||}-\frac{l_{m}}{||l_{m}||}\right\rVert;\;L^{(i)}(\theta):=\sum_{n=1}^{N}\hat{p}^{(i)}_{n}l_{n}(\theta),$
(14)

with $r^{(i)}$ corresponding to the shortest distance from the probability
vector $\hat{p}^{(i)}$ ($i^{th}$ vertex of the convex hull
$\mathcal{D}_{CH}^{(K)}$) to the relative boundary of the feasible region of
$\tilde{w}$. All of the quantities above as well as the quantities needed in
_Algorithm 1_ can be computed from the full sample
($\\{x_{n},y_{n}\\};n=1,2...,N$), the $\nu$-weighted $L2$-norm (10), and inner
product (11).

Algorithm 1 (FW)
---
Given:
$\\{l_{n}(\theta)\\}_{n=1}^{N},\,m,\langle.\rangle,\,\mathcal{D}_{CH}^{(K)}(\\{\hat{p}^{(i)}\\})_{i\in[K]}$
for $i$ in $\\{1,2,...,K\\}$ do
$\quad p_{0}\leftarrow\hat{p}^{(i)}\qquad\qquad\qquad\qquad\qquad\qquad$
$\triangleleft$ initialize $p$ to $i^{th}$ vertex in
$\mathcal{D}_{CH}^{(K)}(\\{\hat{p}^{(i)}\\})_{i\in[K]}$
$\quad\forall n\in[N]$, compute $||l_{n}||$ and $\sigma^{(i)}$ from (III)
$\triangleleft$ compute norms
$\quad f_{0}\leftarrow\operatorname*{arg\,max}_{n\in[N]}\langle
L^{(i)},\frac{l_{n}}{||l_{n}||}\rangle$ $\triangleleft$ select vertex in
$\tilde{w}$-polytope greedily
$\quad\tilde{w}_{0}^{(i)}\leftarrow\frac{\sigma^{(i)}}{||l_{f_{0}}||}\mathbf{1}_{f_{0}}$
$\triangleleft$ initialize $\tilde{w}$ (for $p=\hat{p}^{(i)}$) with full
weight on $f_{0}$
for $t$ in $\\{1,2,...,m-1\\}$ do $\triangleleft$ carry out the following for
$m-1$ iterations for each $\hat{p}^{(i)}$
$\quad\quad f_{t}\leftarrow\operatorname*{arg\,max}_{n\in[N]}\langle
L^{(i)}-L(\tilde{w}^{(i)}_{t-1}),\frac{l_{n}}{||l_{n}||}\rangle$
$\triangleleft$ find FW index at $t^{th}$ iteration
$\quad\quad\gamma_{t}\leftarrow\frac{\langle\frac{\sigma^{(i)}}{||l_{f_{t}}||}l_{f_{t}}-L(\tilde{w}^{(i)}_{t-1}),L^{(i)}-L(\tilde{w}^{(i)}_{t-1})\rangle}{\left\lVert\frac{\sigma^{(i)}}{||l_{f_{t}}||}l_{f_{t}}-L(\tilde{w}^{(i)}_{t-1})\right\rVert^{2}}$
$\triangleleft$ closed from line search for step-size $\gamma$ at $t^{th}$
iteration
$\quad\quad\tilde{w}^{(i)}_{t}\leftarrow(1-\gamma_{t})\tilde{w}^{(i)}_{t-1}+\gamma_{t}\frac{\sigma^{(i)}}{\tilde{\sigma}_{f_{t}}}\mathbf{1}_{f_{t}}$
$\triangleleft$ add or reweight $f_{t}^{th}$ data point
end for
end for
return
$\\{\tilde{w},p\\}=\\{\tilde{w}^{(I)}_{m-1},\hat{p}^{(I)}\\}=\operatorname*{arg\,max}_{i\in[K]}\,\tilde{J}(\tilde{w}^{(i)}_{m-1},\hat{p}^{(i)})$

The proof of the first inequality in (III.1) can be obtained in a
straightforward manner from Jensen’s inequality. In particular, applying
Jensen’s inequality to the convex function $\phi(.):=(.)^{2}$, one has:

$\phi(\mathbf{E}_{\theta\sim\nu}(.))\leq\mathbf{E}_{\theta\sim\nu}\,[\phi(.)]$
(15)

Thus, choosing, $(.)$ as $L(\tilde{w},\theta)-L(\theta)$ and
$L(\tilde{w},\theta)$ respectively, one gets:

$\displaystyle|\left[\mathbf{E}_{\theta\sim\nu}\,(L(\tilde{w},\theta)-L(\theta))\right]|$
$\displaystyle\leq$
$\displaystyle(\mathbf{E}_{\theta\sim\nu}\,(L(\tilde{w},\theta)-L(\theta))^{2})^{1/2}=||L(\tilde{w},\theta)-L(\theta)||_{\nu,2},$
$\displaystyle|\left[\mathbf{E}_{\theta\sim\nu}\,(L(\tilde{w},\theta)\right]|$
$\displaystyle\leq$
$\displaystyle(\mathbf{E}_{\theta\sim\nu}\,(L(\tilde{w},\theta))^{2})^{1/2}=||L(\tilde{w},\theta)-L(\theta)||_{\nu,2}$
(16) $\displaystyle\implies|R(\nu)|$ $\displaystyle=$
$\displaystyle|\mathbf{E}_{\theta\sim\nu}\,L(\theta)|\leq|\left[\mathbf{E}_{\theta\sim\nu}\,(L(\tilde{w},\theta)\right]|+|\left[\mathbf{E}_{\theta\sim\nu}\,(L(\tilde{w},\theta)-L(\theta))\right]|$
(17) $\displaystyle\leq$
$\displaystyle||L(\tilde{w},\theta)||_{\nu,2}+||L(\tilde{w},\theta)-L(\theta)||_{\nu,2}$

where we have used the definition of the expected full sample risk from (3).

The details of the proof of the second inequality in (III.1) are much more
involved and are provided in the appendix. Here we point out some salient
features. First we can write

$||L(\tilde{w},\theta)-L(\theta)||^{2}_{\nu,2}=(\tilde{w}-p)^{T}\tilde{K}(\tilde{w}-p)\,\;\;\mathrm{with}\;\;\tilde{K}_{n,m}:=\langle
l_{n},l_{m}\rangle,$ (18)

For notational simplicity we denote the RHS above as $\tilde{J}(\tilde{w},p)$,
which we interpret as an objective function to optimize. The above problem is
similar to the Hilbert coreset problem posed in (II-B) with $J(w)$ in (II-B)
replaced by a slightly different quantity $\tilde{J}(\tilde{w},p)$. In
contrast to (II-B), the objective function $\tilde{J}(\tilde{w},p)$ in (18)
depends on two vectors $\\{\tilde{w},p\\}$. Since the goal is to provide a
bound for all $D\in\mathcal{D}_{CH}^{(K)}$, one can consider the following
optimization problem by adapting the problem in (II-B) to this case:

$\displaystyle\max_{D=\\{p_{n}\\}\in\mathcal{D}_{CH}^{(K)}}\,\min_{\tilde{w}}\,(\tilde{w}-p)^{T}\tilde{K}(\tilde{w}-p)$
$\displaystyle s.t.$
$\displaystyle\sum_{n=1}^{N}\tilde{w}_{n}||l_{n}||=\sum_{n=1}^{N}p_{n}||l_{n}||,\,\tilde{w}\geq
0$ $\displaystyle s.t.$ $\displaystyle
p\in\mathcal{D}_{CH}^{(K)}(\\{\hat{p}^{(i)}\\})_{i\in[K]},\;\sum_{n=1}^{N}p_{n}=1,p\geq
0.$

Please refer to the appendix for the full proof. An important point to note
about the result in Theorem III.1 is that the upper bound on the full sample
risk denoted by the RHS in (13) approximately goes like $\hat{\beta}^{m}$ even
for moderate number of Frank-Wolfe iterations $m$.

Theorem (III.1) is valid for the $\nu$-weighted $L2$-norm. However, this norm
and the associated inner product involves computing an expectation over the
parameters $\theta$ (see (10) and (11)) which does not admit a closed form
evaluation in general. The expectation can be estimated in an unbiased manner
by the method of random projections, as in [9, 10, 17]. In particular, we
construct a random projection of each vector $l_{n}(\theta);n\in[N]$ into a
finite ($J$) dimensional vector space using $J$ samples
$\theta\sim\nu(\theta)$:

$u_{n}=\sqrt{\frac{1}{J}}[l_{n}(\theta_{1}),l_{n}(\theta_{2}),...,l_{n}(\theta_{J})]^{T},n\in[N],$
(20)

which serves as finite dimensional approximation of $L_{n}(\theta)$. In
particular, for all $n,m\in[N]$, we have:

$\langle l_{n},l_{m}\rangle\approx u_{n}^{T}\,u_{m}.$ (21)

Leveraging similar techniques as in [9], we can construct _Algorithm 2_ , and
arrive at the following result.

Algorithm 2 (RP-FW)
---
Given:
$\\{l_{n}(\theta)\\}_{n=1}^{N},\nu(\theta),J,m,\mathcal{D}_{CH}^{(K)}(\\{\hat{p}^{(i)}\\})_{i\in[K]}$
$\cdot$ $\nu(\theta)$: $\\{\theta_{j}\\}^{J}_{j=1}\sim\nu(\theta)$
$\triangleleft$ take $J$ samples from posterior $\nu(\theta)$
$\cdot$
$u_{n}=\sqrt{\frac{1}{J}}[l_{n}(\theta_{1}),...,l_{n}(\theta_{J})]^{T}$
$\forall n\in[N]$ $\triangleleft$ construct random projection vector for
$l_{n}(\theta)$ using (20)
$\cdot$ return
$\\{\tilde{w},p\\}=FW[u_{n},m,(.)^{T}(.),\mathcal{D}^{(K)}_{CH}(\\{\hat{p}^{(i)}\\})_{i\in[K]}]$
$\triangleleft$ construct coreset with random projection vectors $u_{n}$
by applying _Algorithm 1_ with inner product $\langle
u_{n},u_{m}\rangle=u_{n}^{T}\,u_{m}$

###### Theorem III.2

Given the general setup of theorem III.1, for loss functions $l$ such that
$l_{n}(\mu)\,l_{m}(\mu)$ (with $\mu\sim\nu(\theta)$) is sub-Gaussian with
parameter $\xi^{2}$, for $0<\delta<1$, for
$\mathcal{D}=\mathcal{D}^{(K)}_{CH}(\\{\hat{p}^{(i)}\\})_{i\in[K]}$, the
output of _Algorithm 2_ satisfies with probability $>1-\delta$:

$\displaystyle|R(\nu)|$ $\displaystyle\leq$
$\displaystyle||L(\tilde{w})||_{\nu,2}+||L(\tilde{w})-L||_{\nu,2}$ (22)
$\displaystyle\leq$
$\displaystyle\left(||u(\tilde{w})||^{2}_{\nu,2}+||\tilde{w}||_{1}^{2}\,\sqrt{\frac{2\xi^{2}}{J}\log{\frac{2N^{2}}{\delta}}}\right)^{1/2}+\left(||u(\tilde{w})-u||^{2}_{\nu,2}+||\tilde{w}-p||_{1}^{2}\,\sqrt{\frac{2\xi^{2}}{J}\log{\frac{2N^{2}}{\delta}}}\right)^{1/2},$
$\displaystyle\leq$
$\displaystyle\left(||u(\tilde{w})||^{2}_{\nu,2}+||\tilde{w}||_{1}^{2}\,\sqrt{\frac{2\xi^{2}}{J}\log{\frac{2N^{2}}{\delta}}}\right)^{1/2}+\left(\frac{\hat{\sigma}\,\hat{\eta}\,\hat{\bar{\eta}}\,\hat{\beta}}{\sqrt{\hat{\bar{\eta}}^{2}\,\hat{\beta}^{-2(m-1)}+\hat{\eta}^{2}\,(m-1)}}+||\tilde{w}-p||_{1}^{2}\,\sqrt{\frac{2\xi^{2}}{J}\log{\frac{2N^{2}}{\delta}}}\right)^{1/2},$
(23)

The quantities $\\{\hat{\sigma},\hat{\eta},\,\hat{\bar{\eta}},\hat{\beta}\\}$
are defined as before. The sub-Gaussianity parameter $\xi$ can be estimated
from the loss functions. For the $0-1$ loss function, a simple bound is
$\xi\leq 1/2$. Also note that
$\\{\tilde{w},p\\}=\\{\tilde{w}_{m-1}^{(I)},\hat{p}^{(I)}\\}$ where $I\in[K]$
is the index that maximizes $\tilde{J}(\tilde{w}^{(i)}_{M-1},\hat{p}^{(i)})$
(see _Algorithm 1_). The details of the proof can be found in the appendix.

##  IV Experiments and Instantiation of Bounds

In this section, we present our experimental setup for computing Hilbert
coreset bounds for three datasets: i) Binary MNIST, and ii) CIFAR-100. MNIST
is a very well known image classification dataset and is usually treated as a
multiclass problem. We follow the results of [4] and create a binary
classification dataset out of it by labeling the digits 0-4 with class ‘-1’
and digits 5-9 with class ‘1’. The CIFAR-100 datasets is also a well known
dataset with hundred different classes of images.

The distribution on the data in each case is assumed to be a uniform
distribution for simplicity, with $p_{n}=\frac{1}{N}$ for a full sample of
size $N$. In other words, the class of data distributions, e.g. the convex
hull $\mathcal{D}_{CH}^{(K)}$ is simply a point. This can be extended in a
straightforward manner to cases where $\mathcal{D}_{CH}^{(K)}$ using the
methods described in the appendix. Finally, the loss function $l$ used for the
computation of bounds is the standard 0-1 loss used in classification. Further
experimental details for each of these datasets are provided below.

###  IV-A Binary MNIST

The MNIST dataset is obtained from Keras. This dataset is split into the
training set (60000 images) and test set (10000 images). We convert this into
a binary classification dataset as discussed above. The images are
preprocessed by scaling them to float values in the range $[0,1]$. Keras is
used to construct a simple neural network binary classifier, similar to that
in [4]. In particular, the hypothesis class $\mathcal{H}$ for this dataset is
taken to be a feed-forward neural network which has three hidden layers each
with 600 nodes and a ReLu activation function. The output layer of the network
uses a sigmoid activation function. The weights at each node are initialized
according to a truncated normal distribution with mean $0$ and standard
deviation $0.04$. The biases of the first hidden layer are 0.1, whereas the
biases of the second hidden layer and output layer are initialized to 0.

The model is trained on different random subsets of the MNIST dataset using
the cross-entropy loss (however the bounds will be provided with respect to
the 0-1 loss).The first checkpoint uses 500 data points from the training set.
The second checkpoint used 3000 data points, including the 500 from the first
checkpoint.The third checkpoint uses 10000 data points, including the 3000
already used, and the final checkpoint uses 30000 data points. The models were
trained for 120 epochs for each checkpoint.

###  IV-B CIFAR-100

For this dataset, we use a state-of-the art neural network architecture as the
hypothesis class and a training algorithm is described in Xie et. al. [18].
The reason for doing this is the following. As a performer on DARPA Learning
with Less Labeling (LwLL) program, the goal of our team was to develop novel
theoretical bounds for other performers who are working on the model-building
and training of new algorithms for efficient learning of classifiers with
fewer labels. As part of this exercise, we approached the authors in [18] (who
are part of the Information Sciences Institute (ISI) team) for the training
models developed by them and developed theoretical bounds for those models
(the associated hypothesis classes $\mathcal{H}$). As we will show below, we
obtained tight non-vacuous upper bounds for such hypothesis classes as well,
indicating the power and diversity of the proposed approach.

We summarize the training algorithm used in [18] for completeness. [18] uses
mutual-information based unsupervised $\&$ semi-supervised Concurrent learning
(MUSCLE), a hybrid learning approach that uses mutual information to combine
both unsupervised and semisupervised learning. MUSCLE can be used as a
standalone training scheme for neural networks, and can also be incorporated
into other learning approaches. It is shown that the proposed hybrid model
outperforms state of the art for the CIFAR-100 dataset as well as others. The
basic idea of MUSCLE is to combine combine the loss from (semi)-supervised
learning methods including real and pseudo labels with a loss term arising
from the mutual information between two different transformed versions of the
input sample ($x_{i}$). Furthermore, a consistency loss that measures the
consistency of the prediction on an input sample and an augmented sample, can
be added as a regularization term. Please refer to [18] and references therein
for further details about the training algorithm and neural network models
(hypothesis classes).

The model is trained with different random subsets of the data, at the 500,
2500, 4000 and 10000 checkpoints.

###  IV-C Choice of Prior and Posterior Distributions

As explained in section I, the theoretical bounds we develop are stochastic
bounds. So, one has to choose an appropriate prior distribution (any
distribution independent of the data) and a posterior distribution (that
depends on the training data and the training algorithm in general). In our
analysis, we choose the prior distribution to be
$\pi_{0}=N(0,\sigma^{2}\cdot\mathbf{I})$. Inspired by the analysis of Banerjee
et al [19], we use a training data-dependent posterior
$\nu_{S}=N(\theta^{{\dagger}},\Sigma^{S}_{\theta^{{\dagger}}})$, where
$\theta^{{\dagger}}$ denotes the parameters _learned_ by training the models
(in the previous subsections) on the training set $S$, and the covariance
$\Sigma_{\theta^{{\dagger}}}$ is defined as follows:

$\displaystyle\Sigma^{S}_{\theta^{{\dagger}}}$ $\displaystyle=$ $\displaystyle
diag(\sigma_{j}^{2});\;\sigma_{j}^{-2}:=\max\,[H^{S}_{j,j}(\theta^{{\dagger}}),\sigma^{-2}],\;\;j=1,2,..,dim(\Theta)$
$\displaystyle\mathrm{where}\;\;H^{S}(\theta^{{\dagger}})$ $\displaystyle:=$
$\displaystyle\frac{1}{N_{S}}\sum_{i=1}^{N_{S}}\,\nabla^{2}\,l(y_{i},\phi(\theta^{{\dagger}},x_{i}))\;\;\mathrm{(Diagonal\;Hessian)}$
(24)

The anisotropy in the posterior $\nu_{S}$ can be understood as follows. For
directions $j$ in parameter space where $H^{S}_{j,j}(\theta^{{\dagger}})$ is
_larger_ than $\sigma^{-2}$, the posterior variance is determined by the
Hessian, and will be small enough such that $\theta_{j}$ will not deviate far
from the learned parameters $\theta^{{\dagger}}_{j}$ when sampling from the
posterior $\nu_{S}$. On the other hand, for directions $j$ in parameter space
where $H^{S}_{j,j}(\theta^{{\dagger}})$ is _smaller_ than $\sigma^{-2}$, the
posterior variance $\sigma_{j}^{2}$ is the same as that of the prior
$\sigma^{2}$. Thus, the posterior standard deviation is upper bounded by the
prior standard deviation and different weights and biases are assumed to be
independent. One can optimize over the prior standard deviation $\sigma$ to
get a good bound.

###  IV-D Bound Calculation Workflow

Given the data points $(x_{n},y_{n}),n=1,2,..,N_{S}$ of the training set
checkpoint $S$ of a particular dataset, a trained model (with the learned
parameters $\theta^{{\dagger}}$) for that checkpoint, and a prior distribution
$\pi_{0}=N(0,\sigma^{2}\cdot\mathbf{I})$ with an appropriately chosen
$\sigma$, the workflow for computing the coreset bounds corresponding to the
training set $S$ is as follows:

  * •

We compute the diagonal Hessian of the (cross-entropy) loss for the trained
model at the trained parameters $H^{S}_{j,j}(\theta^{{\dagger}})$ for the
training set $S$ via the methods described in [20]. Then we compute the
posterior variance($\sigma_{j}^{2}$) as described in section IV-C. This gives
the data dependent posterior $\nu_{S}$ for the training checkpoint $S$.

  * •

Given $\nu_{S}$, we compute an approximation to the inner product
$\langle\cdot\rangle$ by the method of _random features projection_ [17, 9] as
in _Algorithm 2_ which converge to the true inner products with high
probability as the number of samples from the posterior ($J$) increases (see
Theorem III.2).

  * •

From the inner products computed with the posterior for the training set $S$,
and using the full dataset of size $N$, we compute all the constants
$\\{\hat{\sigma},\hat{\eta},\,\hat{\bar{\eta}},\hat{\beta}\\}$ in Theorems
III.1 and III.2 and then implement _Algorithm 1_ with number of iterations $m$
equal to training set checkpoint size $N_{S}$ to find the coreset weights
$\tilde{w}^{(I)}_{m-1}$ and the appropriate probability vector
$\hat{p}^{(I)}$.

  * •

Finally, using the coreset weights and the inner product $(.)^{T}\,(.)$, we
compute _coreset upper bound-I_ (first term in the RHS of (22)) as well as
_coreset upper bound-II_ (RHS of (23)).

The above workflow was carried out for all different training checkpoints (of
increasing sizes) for both datasets mentioned in sections IV-A and IV-B.

###  IV-E Experimental Results

|
---|---
Figure 1: Left: _Coreset Upper Bound - I_ (RHS of (22), shown in blue) and
_Coreset Upper Bound -II_ (RHS of (23), shown in red) on the expected full
sample risk for the three-layer ReLU DNN model for four training checkpoints
of the Binary MNIST dataset with prior standard deviation $\sigma=10^{-4}$,
based on Theorem III.2. The expected empirical risk (EER) of the trained model
on the full sample is shown in green. Right: _Coreset Upper Bound - I_ (RHS of
(22), shown in blue) and _Coreset Upper Bound -II_ (RHS of (23), shown in red)
on the expected full sample risk for the Xie _et al_ [18] model for four
training checkpoints of the CIFAR-100 dataset with prior standard deviation
$\sigma=10^{-4}$, based on Theorem III.2. The expected empirical risk (EER) of
the trained model on the full sample is shown in green.

For our numerical results, we choose $J$, the number of samples from the
posterior (see section IV-D), as 7000. We experimented with other values of
$J$ and found that the numerical results for the bounds were stable with
$J\gtrsim 7000$. We also experimented with the choice of the prior standard
deviation $\sigma$. As discussed in section IV-C, the posterior standard
deviation is upper bounded by $\sigma$. So, we would like to choose a low-
value of $\sigma$ to come up with an _effectively deterministic_ bound. After
some experimentation, we chose a value of $\sigma=10^{-4}$ that seemed to
provide good bounds which were effectively deterministic in the sense that the
expected empirical risk (EER) is very close to the empirical risk. For each
dataset, we use the number of Frank-Wolfe iterations $m$ (in _Algorithm 1_)
equal to the size of the training set checkpoint $S$. We find that the size of
the coreset (the number of non-zero weights) is _much smaller_ than the number
of Frank-Wolfe iterations or the size of the training set $S$. For example,
for the CIFAR-100 dataset, the coreset size is only $\sim 150$, around two
orders of magnitude smaller than the size of the largest training set
checkpoint $S=10000$. This is also consistent with results in [9, 10].

The left plot in Figure 1 shows that the bounds for the MNIST model are very
tight. For example, the empirical expected full sample risk at the first
checkpoint with 500 samples (0.49) is quite close to the _coreset upper bound-
I_ (0.51) and the _coreset upper bound-II_ (0.56). We get tight non-vacuous
bounds for the DNN hypothesis class used in section IV-A for all training
checkpoints. Also, the coreset upper bounds obtained for the training
checkpoint with 30000 samples ($0.145\,\&\,0.154$) are better than that
obtained in [4] for 55000 samples (0.207). As discussed earlier, we use the
same model and training procedure as [4] to allow a fair comparison.

The right plot in Figure 1 shows the results for the CIFAR-100 model with the
state-of-the-art neural network architecture and training algorithm as in
[18]. For this model, the _coreset upper bound-I_ is non-vacuous at all
training checkpoints, while the _coreset upper bound-II_ is non-vacuous for
training set sizes larger than around 8000. Thus we are able to compute
meaningful and effectively deterministic upper bounds for the full sample risk
for complex datasets (like CIFAR-100) and state-of-the-art hypothesis classes
and training algorithms, as in [18] in section IV-B.

The constants in the bounds in Theorems III.1 and III.2 depend on the
posterior distribution $\nu$. As discussed in section IV-D, the posterior
distribution $\nu_{S}$ is computed for each training checkpoint $S$.
Therefore, the constants in the bounds plotted in Figure 1 also change with
each training checkpoint $S$. An important point to note is that the bounds
for both datasets/models decrease significantly with the size of the training
set $S$. This is in contrast to some other recent bounds for DNNs that may
even increase with the size of the training set [3, 11, 12], see discussion in
[13].

Finally, as pointed out in the beginning of section IV, we have assumed a
uniform distribution on the data for simplicity implying that the class of
data distributions $\mathcal{D}$ is trivial. However, it is straightforward to
extend these results to non-trivial classes like $\mathcal{D}_{CH}^{(K)}$,
using the methods described in the appendix. The expected empirical full
sample risk will still be upper bounded by the bounds computed for these data
distribution classes as long as the uniform distribution belongs to
$\mathcal{D}_{CH}^{(K)}$.

##  V Conclusions

In this paper, we have developed the formalism for computing upper bounds on
the _expected full sample risk_ for supervised classification tasks within a
transductive setting using the method of Hilbert coresets. We compute explicit
numerical bounds for complex datasets and complex hypothesis classes (such as
state-of-the-art DNNs) and compare these to the performance of modern training
algorithms. The bounds obtained are tight and meaningful. The bounds we
develop are stochastic, however they are effectively deterministic for
appropriate choice of priors. The bounds are _non-uniform_ in the hypothesis
space, since the precise numerical values of the bound depend on the learned
weights of the training algorithm ($\theta_{\star}$). Finally, the bounds we
develop decrease significantly with the size of the training set ($S$) as can
be seen from section IV-E.

There are a couple of avenues for future research. First, it should be
possible to derive generalization bounds for other classes of data
distributions $\mathcal{D}$, such as in which the probability vector for the
full sample lies inside an ellipsoid. Second, as mentioned in section IV-E,
the coreset size is much smaller than the size of the training set $S$ for a
given training checkpoint. This suggests that it might be possible to develop
a model-free (agnostic) active learning approach leveraging the coresets
framework. For example, given a prior and an initial (presumably small) set of
labeled data, one can use the coresets approach to construct an initial
posterior on the hypothesis space. Using this posterior, one can make
predictions on the unlabeled data and compute a _pseudo risk_ for the entire
dataset. Then one could use this _pseudo-risk_ to construct the next set of
coreset weights and update the coreset weighted risk. At each step, one should
also update the posterior appropriately. The goal would be to show that this
procedure converges and provides an upper bound on the expected true risk. The
hope is that this should provide much tighter generalization and sample
complexity bounds.

##  VI Acknowledgement

The material in this paper is based upon work supported by the United States
Air Force and DARPA under Contract No. FA8750-19-C-0205 for the DARPA I20
_Learning with Less Labeling (LwLL)_ program. S.D., P.K., and R.K.P are
grateful to Prof. T. Broderick for providing inspiration to the ideas in this
paper as well as providing feedback on the manuscript. The authors would like
to particularly thank Prof. V. Chandrasekaran and Prof. A. Willsky for
insightful discussions and comments on the manuscript.

###  -A Proof of Theorem III.1 for
$\mathcal{D}=\mathcal{D}_{CH}^{(K)}(\hat{p}^{(I)}),I=1,2..,K$

In this case, we define $\mathcal{D}$ as the convex hull of $K$ points in the
probability simplex characterized by $N$ points.

###### Lemma .1

Let $\mathcal{D}_{CH}^{(K)}(\hat{p}^{(I)})$ be the convex hull of
$\\{\hat{p}^{(I)};1\leq I\leq K\\}$, where each $\hat{p}^{(I}$ is an
$N$-dimensional probability vector with
$\sum_{n=1}^{N}\,\hat{p}^{(I)}_{n}=1;\hat{p}^{(I)}_{n}>0\;\forall n\in[N]$.
Then the optimal solution to
$\max_{p\in\mathcal{D}_{CH}^{(K)}}\,\tilde{J}(\tilde{w},p)$ for a fixed
$\tilde{w}$ is $\hat{p}^{(k)}$ for some $k$ in $\\{1,2,..,K\\}$.

###### Proof:

Let $p_{0}$ be any point which is _not_ a vertex in the convex hull and fix
$\tilde{w}$. Because $\tilde{J}$ is quadratic, we can represent the function
as

$\tilde{J}(\tilde{w},p)=\tilde{J}(\tilde{w},p_{0})+\nabla\tilde{J}(\tilde{w},p_{0})^{T}(p-p_{0})+(p-p_{0})^{T}\tilde{K}(p-p_{0})$

Since $p_{0}$ is not an extreme point, there is a vector $d$ such that both
$p_{0}+d$ and $p_{0}-d$ are in the convex hull. Clearly, either
$\nabla\tilde{J}(\tilde{w},p_{0})^{T}d\geq 0$ or
$-\nabla\tilde{J}(\tilde{w},p_{0})^{T}d\geq 0$. So without loss of generality,
we can claim there is a vector $d$ such that
$\nabla\tilde{J}(\tilde{W},p_{0})^{T}d\geq 0$. Furthermore, because
$\tilde{J}$ is quadratic, $\tilde{K}$ is positive definite. So we conclude
$p_{0}$ is suboptimal because

$\tilde{J}(\tilde{w},p_{0})+\nabla\tilde{J}(\tilde{w},p_{0})^{T}(p-p_{0})+(p-p_{0})^{T}\tilde{K}(p-p_{0})>\tilde{J}(\tilde{w},p_{0})$

∎

We now provide proofs pertaining to _Algorithm 1_ for
$\mathcal{D}=\mathcal{D}_{CH}^{(K)}$.

###### Lemma .2

_Algorithm 1_ is equivalent to the Frank-Wolfe coreset construction of [9] for
$K=1$.

###### Proof:

For $K=1$, the outer maximization is trivial. Since all elements of
$\hat{p}^{(1)}$ are non-zero, so we can convert back to $w$. Notice that the
objective then becomes

$\tilde{J}(\tilde{w},p)=(\tilde{w}-\hat{p}^{(1)})^{T}\tilde{K}(\tilde{w}-\hat{p}^{(1)})=(w-1)^{T}\hat{K}(w-1)$
(25)

after writing $w=\tilde{w}\,\hat{p}^{(1)}$ and identifying
$\hat{p}^{(1)}_{n}\,l_{n}$ with $L_{n}$. Thus, we have:

$\hat{p}^{(1)}_{n}||l_{n}||=||L_{n}||,$ (26)

and the solution to the original Frank-Wolfe algorithm in [9] is also a
solution to _Algorithm 1_. ∎

Proof of Theorem III.1 for $\mathcal{D}_{CH}^{(K)};K>1$:

Let us first define a number of constants that would be helpful later:

$\displaystyle\tilde{\sigma}_{n}$ $\displaystyle=$
$\displaystyle||l_{n}||;\;\sigma^{(i)}=\sum_{n=1}^{N}\hat{p}^{(i)}_{n}\tilde{\sigma}_{n};\;L^{(i)}=\sum_{n=1}^{N}\hat{p}^{(i)}_{n}l_{n}$
$\displaystyle\hat{\bar{\eta}}^{2}$ $\displaystyle=$
$\displaystyle\max_{n,m}\left\lVert\frac{l_{n}}{||l_{n}||}-\frac{l_{m}}{||l_{m}||}\right\rVert;\;{\eta^{(i)}}^{2}=1-\frac{||L^{(i)}||^{2}}{(\sigma^{(i)})^{2}};\;{\beta^{(i)}}^{2}=1-\frac{(r^{(i)})^{2}}{(\sigma^{(i)})^{2}\bar{\eta}^{2}},$
(27)

where $r^{(i)}$ is the shortest distance from the solution $\hat{p}^{(i)}$ to
the relative boundary of the feasible region for $\tilde{w}$. The first
inequality in Theorem III.1 was proven in the main text. We now prove the
second inequality (see (13)) in Theorem III.1.

###### Theorem .1

Let $\hat{\sigma}=\max_{i\in[K]}\sigma^{(i)}$,
$\hat{\eta}=\max_{i\in[K]}\eta^{(i)}$ , and
$\hat{\beta}=\max_{i\in[K]}\beta^{(i)}$. Then, with
$L=\sum_{n=1}^{N}\,p_{n}l_{n}$ and
$L(\tilde{w})=\sum_{n=1}^{N}\,\tilde{w}_{n}\,l_{n}=\sum_{n-1}^{N}\,w_{n}p_{n}l_{n}$,
the coreset constructed with _Algorithm 1_ gives rise to a Hilbert coreset of
size at most $m$ with

$||L-L(\tilde{w})||\leq\frac{\hat{\sigma}\hat{\eta}\hat{\bar{\eta}}\hat{\beta}}{\sqrt{\hat{\bar{\eta}}^{2}\hat{\beta}^{-2(m-1)}+\hat{\eta}^{2}(m-1)}}$
(28)

###### Proof:

In order to prove theorem (.1), we first borrow two lemmas from [9] which will
be useful shortly.

###### Lemma .3

If $\hat{p}^{(i)}$ has no non-zero entries, then
$L^{(i)}\equiv\sum_{n=1}^{N}\,\hat{p}^{(i)}_{n}\,l_{n}$ is in the relative
interior of the convex hull:
$Conv\\{\frac{\sigma^{(i)}}{\tilde{\sigma}_{n}}l_{n}\\}_{n=1}^{N}$

The proof the above is a simple consequence of Lemma A.5 in [9] and by the
following identification:

$\frac{\sigma^{(i)}}{\tilde{\sigma}_{n}}l_{n}=\frac{\sigma^{(i)}}{\hat{p}^{(i)}_{n}||l_{n}||}\hat{p}^{(i)}_{n}l_{n}=\frac{\sigma^{(i)}}{||L_{n}||}L^{(i)}_{n}$
(29)

###### Lemma .4

Let $\\{x_{n}\\}$ be a sequence of real numbers and $0\leq\alpha\leq 1$ which
satisfy $x_{k+1}\leq\alpha x_{k}(1-x_{k})$ for all $k$. Then for all $n$,

$x_{n}\leq\frac{x_{0}}{\alpha^{-n}+nx_{0}}$ (30)

This is listed as lemma A.6 in [9].

From Lemma .1, we know the outer maximization is solved by one of the vertices
of the convex hull $\mathcal{D}_{CH}^{(K)}$. Therefore, following _Algorithm
1_ in Table 1, we first select a vertex $\hat{p}^{(i)}$ from the $K$ vertices
in $\mathcal{D}_{CH}^{(K)}$ at random. Then for that $\hat{p}^{(i)}$, we use
the Frank-Wolfe algorithm to minimize $\tilde{w}$ subject to a convex polytope
constraint on $\tilde{w}$ as in (III). This gives a value of
$\tilde{w}_{t}^{(i)}$ after $t=m-1$ iterations corresponding to
$\hat{p}^{(i)}$. The same procedure is then carried out for all $K$ vertices
of $\mathcal{D}_{CH}^{(K)}$, and we choose the index $(I)$ that maximizes the
objective $\tilde{J}(\tilde{w}_{t}^{(i)},\hat{p}^{(i)});i=1,2,..,K$.

We now show that the $\tilde{w}_{t}^{(I)}$ obtained for the index $(I)$
satisfies Theorem III.1. For simplicity, we omit the index $(I)$ from
$\tilde{w}_{t}^{(I)}$ and just use $\tilde{w}_{t}$. We first initialize the
Frank-Wolfe algorithm by choosing an appropriate $\tilde{w}_{0}$ as follows:

$\displaystyle f_{0}$ $\displaystyle=$
$\displaystyle\operatorname*{arg\,max}_{n\in[N]}\,\langle\frac{l_{n}}{||l_{n}||},\sum_{n=1}^{N}\,{\hat{p}_{n}^{(I)}}\,l_{n}\rangle$
$\displaystyle\tilde{w}_{0}$ $\displaystyle=$
$\displaystyle\frac{\sigma^{(I_{0})}}{||l_{f_{0}}||}\mathbf{1}_{f_{0}}=\frac{\sigma^{(I)}}{||\tilde{\sigma}_{f_{0}}||}\mathbf{1}_{f_{0}}$
(31)

We now show that
$\tilde{J}(\tilde{w}_{0},\hat{p}^{(I)})\leq\sigma^{2}\eta^{2}$. Notice that
for any $\xi\in\mathbb{R}^{N}$ such that $\sum_{n=1}^{N}\xi_{n}=1$,

$\displaystyle\frac{\tilde{J}(\tilde{w}_{0},\hat{p}^{(I)})}{{\sigma^{(I)}}^{2}}$
$\displaystyle=$ $\displaystyle
1-2<\frac{l_{f_{0}}}{\tilde{\sigma_{f_{0}}}},\frac{L^{(I)}}{\sigma^{(I)}}>+\frac{||L^{(I)}||}{{\sigma^{(I)}}^{2}}$
(32) $\displaystyle\leq$ $\displaystyle
1-2\sum_{n=1}^{N}\xi_{n}<\frac{l_{n}}{\sigma_{n}},\frac{L^{(I)}}{\sigma^{(I)}}>+\frac{||L^{(I)}||}{{\sigma^{(I)}}^{2}},$

since $l_{f_{0}}$ maximizes the inner product above. Choosing
$\xi_{n}=\frac{\hat{p}^{(I)}\tilde{\sigma}_{n}}{\sigma^{(I)}}$ yields

$\displaystyle\frac{\tilde{J}(\tilde{w}_{0},\hat{p}^{(I)})}{{\sigma^{(I)}}^{2}}$
$\displaystyle\leq$ $\displaystyle
1-2\sum_{n=1}^{N}<\frac{\hat{p}_{n}^{(I)}l_{n}}{\sigma^{(I)}},\frac{L^{(I)}}{\sigma^{(I)}}>+\frac{||L^{(I)}||}{{\sigma^{(I)}}^{2}}$
(33) $\displaystyle=$ $\displaystyle
1-\frac{||L^{(I)}||}{{\sigma^{(I)}}^{2}}:={\eta^{(I)}}^{2}$

Thus,
$\tilde{J}(\tilde{w}_{0},\hat{p}^{(I)})\leq{\sigma^{(I)}}^{2}{\eta^{(I)}}^{2}\leq\hat{\sigma}^{2}\hat{\eta}^{2}$.

Next, we will look at the sequence of values,
$\tilde{J}(\tilde{w}_{t},\hat{p}^{(I)})$. The goal is to use Lemma .4 on this
sequence to get the desired bound. The Frank-Wolfe algorithm uses the update
$\tilde{w}_{t+1}=\tilde{w}_{t}+\gamma^{(I)}_{t}\,d_{(I)}$, where

$\displaystyle d_{(I)}$ $\displaystyle=$
$\displaystyle\operatorname*{arg\,min}_{s}[\nabla
J(\tilde{w}_{t},\hat{p}^{(I)})^{T}s]-\tilde{w}_{t}$
$\displaystyle\mathrm{or},\;d_{(I)}$ $\displaystyle=$
$\displaystyle\operatorname*{arg\,min}_{s}[2(\tilde{w}_{t}-\hat{p}^{(I)})^{T}\,\tilde{K}\,s]-\tilde{w}_{t}=\operatorname*{arg\,min}_{s}[2(\tilde{w}_{t}-\hat{p}^{(I)})^{T}\,\tilde{K}\,\frac{\sigma^{(I)}}{\tilde{\sigma}_{n}}\,\mathbf{1}_{n}]-\tilde{w}_{t},$
(34)

where we have used the fact that the minimizer of a linear function in $s$
over a convex polytope will lie at one of the vertices
$\frac{\sigma^{(I)}}{\tilde{\sigma}_{n}}\,\mathbf{1}_{n};\;n\in[N]$. Then it
is straightforward to show that:

$d_{(I)}=\frac{\sigma^{(I)}}{\tilde{\sigma}_{f_{t}}}\mathbf{1}_{f_{t}}-\tilde{w}_{t};\;\;f_{t}=\operatorname*{arg\,max}_{n\in[N]}\langle
L-L(\tilde{w}_{t}),\frac{l_{n}}{\tilde{\sigma}_{n}}\rangle$ (35)

This gives:

$\tilde{J}(\tilde{w}_{t+1},\hat{p}^{(I)})=\tilde{J}(\tilde{w}_{t},\hat{p}^{(I)})+2\gamma^{(I)}_{t}d_{(I)}\tilde{K}(\tilde{w}_{t}-\hat{p}^{(I)})+{\gamma_{t}^{(I)}}^{2}(d_{(I)})^{T}\tilde{K}d_{(I)}$
(36)

According to the Frank-Wolfe algorithm, we choose $\gamma_{t}^{(I)}$ which
minimizes this quantity in the domain [0,1], which gives rise to:

$\gamma_{t}^{(I)}=\frac{d_{(I)}\tilde{K}(\hat{p}^{(I)}-\tilde{w}_{t})}{d_{(I)}\tilde{K}d_{(I)}}$
(37)

Substituting this in (36) gives:

$\displaystyle\tilde{J}(\tilde{w}_{t+1},\hat{p}^{(I)})$ $\displaystyle=$
$\displaystyle\tilde{J}(\tilde{w}_{t},\hat{p}^{(I)})-\frac{((d_{t}^{(I)})^{T}\tilde{K}(\hat{p}^{(I)}-\tilde{w}_{t}))^{2}}{(d_{t}^{(I)})^{T}\tilde{K}d_{t}^{(I)}}$
(38) $\displaystyle=$
$\displaystyle\tilde{J}(\tilde{w}_{t},\hat{p}^{(I)})\left(1-\langle\frac{\frac{\sigma^{(I)}}{\tilde{\sigma}_{f_{t}}}l_{f_{t}}-L(\tilde{w}_{t})}{\left\lVert\frac{\sigma^{(I)}}{\tilde{\sigma}_{f_{t}}}l_{f_{t}}-L(\tilde{w}_{t})\right\rVert},\frac{L^{(I)}-L({\tilde{w}_{t})}}{\left\lVert
L^{(I)}-L({\tilde{w}_{t})}\right\rVert}\rangle\right)$

Now, notice that since Lemma .3 says that $L^{(I)}$ is in the relative
interior of the convex polytope (in $\tilde{w}$ space) for a given
$\hat{p}^{(I)}$, there exists an $r^{(I)}>0$ such that all $\tilde{w}$ that
are feasible,

$L(\tilde{w})+(\left\lVert
L^{(I)}-L(\tilde{w})\right\rVert+r^{(I)})\frac{L^{(I)}-L(\tilde{w})}{\left\lVert
L^{(I)}-L(\tilde{w})\right\rVert}$ (39)

is also feasible. Furthermore, the Frank-Wolfe vertex
$\frac{\sigma^{(I)}}{\tilde{\sigma}_{f_{t}}}l_{f_{t}}$ maximizes the inner
product $\langle
L^{(I)}-L(\tilde{w}_{t}),L(\tilde{w})-L(\tilde{w}_{t})\rangle$ over all
feasible $\tilde{w}$. This implies that

$\displaystyle\langle\frac{\frac{\sigma^{(I)}}{\tilde{\sigma}_{f_{t}}}l_{f_{t}}-L(\tilde{w}_{t})}{\left\lVert\frac{\sigma^{(I)}}{\tilde{\sigma}_{f_{t}}}l_{f_{t}}-L(\tilde{w}_{t})\right\rVert},\frac{L^{(i)}-L(\tilde{w}_{t})}{\left\lVert
L^{(i)}-L(\tilde{w}_{t})\right\rVert}\rangle$ $\displaystyle\geq$
$\displaystyle\langle\frac{(\left\lVert
L^{(I)}-L(\tilde{w}_{t})\right\rVert+r^{(I)})\frac{L^{(I)}-L(\tilde{w}_{t})}{\left\lVert
L^{(I)}-L(\tilde{w}_{t})\right\rVert}}{\left\lVert\frac{\sigma^{(I)}}{\tilde{\sigma}_{f_{t}}}l_{f_{t}}-L(\tilde{w}_{t})\right\rVert},\frac{L^{(I)}-L(\tilde{w}_{t})}{\left\lVert
L^{(I)}-L(\tilde{w}_{t})\right\rVert}\rangle$ $\displaystyle=$
$\displaystyle\frac{\sqrt{\tilde{J}(\tilde{w}_{t},\hat{p}^{(I)}})+r^{(I)}}{\left\lVert\frac{\sigma^{(I)}}{\tilde{\sigma}_{f_{t}}}l_{f_{t}}-L(\tilde{w}_{t})\right\rVert}$
$\displaystyle\geq$
$\displaystyle\frac{\sqrt{\tilde{J}(\tilde{w}_{t},\hat{p}^{(I)}})+r^{(I)}}{\sigma^{(I)}\bar{\eta}^{(I)}}$

Combining equations (36) and (-A) gives:

$\displaystyle\tilde{J}(\tilde{w}_{t+1},\hat{p}^{(I)})$ $\displaystyle\leq$
$\displaystyle\tilde{J}(\tilde{w}_{t},\hat{p}^{(I)})\left(1-\left(\frac{\sqrt{\tilde{J}(\tilde{w}_{t},\hat{p}^{(I)}})+r^{(I)}}{\sigma^{(I)}\bar{\eta}^{(I)}}\right)^{2}\right)$
(41) $\displaystyle\leq$
$\displaystyle{\beta^{(I)}}^{2}\tilde{J}(\tilde{w}_{t},\hat{p}^{(I)})\left(1-\frac{\tilde{J}(\tilde{w}_{t},\hat{p}^{(I)})}{(\sigma^{(I)})^{2}(\bar{\eta}^{(I)})^{2}(\beta^{(I)})^{2}}\right)$
$\displaystyle\leq$
$\displaystyle\hat{\beta}^{2}\tilde{J}(\tilde{w}_{t},\hat{p}^{(I)})\left(1-\frac{\tilde{J}(\tilde{w}_{t},\hat{p}^{(I)})}{\hat{\sigma}^{2}\hat{\bar{\eta}}^{2}\hat{\beta}^{2}}\right)$

Setting $\alpha=\hat{\beta}^{2}$ and
$x_{t}:=\frac{\tilde{J}(\tilde{w}_{t},\hat{p}^{(I)})}{\hat{\sigma}^{2}\hat{\bar{\eta}}^{2}\hat{\beta}^{2}}$.
Then by Lemma .4,

$\frac{\tilde{J}(\tilde{w}_{t},\hat{p}^{(i)})}{\hat{\sigma}^{2}\hat{\bar{\eta}}^{2}\hat{\beta}^{2}}\leq\frac{\frac{\tilde{J}(\tilde{w}_{0},\hat{p}^{(i)})}{\hat{\sigma}^{2}\hat{\bar{\eta}}^{2}\hat{\beta}^{2}}}{\hat{\beta}^{(-2t)}+t\frac{\tilde{J}(\tilde{w}_{0},\hat{p}^{(i)})}{\hat{\sigma}^{2}\hat{\bar{\eta}}^{2}\hat{\beta}^{2}}}$
(42)

Using the fact that $\frac{a}{a+b}$ is monotonically increasing in $a$ when
$a,b>0$ and
$\tilde{J}(\tilde{w}_{0},\hat{p}^{(i)})\leq\hat{\sigma}^{2}\hat{\eta}^{2}$,
one gets:

$\displaystyle\tilde{J}(\tilde{w}_{t},\hat{p}^{(I)})$ $\displaystyle\leq$
$\displaystyle\frac{\hat{\sigma}^{2}\hat{\bar{\eta}}^{2}\hat{\beta}^{2}}{\hat{\eta}^{2}\,t+\hat{\beta}^{-2(t-1)}\hat{\bar{\eta}}^{2}}$
$\displaystyle\mathrm{or},\;||L-L(\tilde{w}_{t}||$ $\displaystyle\leq$
$\displaystyle\frac{\hat{\sigma}^{2}\hat{\bar{\eta}}^{2}\hat{\beta}^{2}}{\sqrt{\hat{\eta}^{2}\,(m-1)+\hat{\beta}^{-2(m-2)}\hat{\bar{\eta}}^{2}}}$
(43)

after $t=m-1$ iterations. ∎

###  -B Proof of Theorem III.2 for
$\mathcal{D}=\mathcal{D}^{(K)}_{CH}(\\{\hat{p}^{(i)}\\})_{i\in[K]}$

###### Proof:

The proof of this closely follows Theorem 5.2 in [9]. Given a posterior
distribution $\nu(\theta)$, we have:

$\langle
l_{m},l_{n}\rangle:=\mathbb{E}_{\theta\sim\nu}\,[l_{n}(\theta)l_{m}(\theta)].$
(44)

Given the true vector as $l_{n}$ and it’s random projection as $u_{n}$, one
can write:

$\displaystyle\mathbb{P}[\max_{m,n}|\langle l_{m},l_{n}\rangle-
u_{m}^{T}\,u_{n}|\geq\epsilon]$ $\displaystyle\leq$
$\displaystyle[\sum_{m,n}\,\mathbb{P}[|\langle l_{m},l_{n}\rangle-
u_{m}^{T}\,u_{n}|\geq\epsilon]$ (45) $\displaystyle\leq$ $\displaystyle
N^{2}\,\max_{m,n}\mathbb{P}[|\langle l_{m},l_{n}\rangle-
u_{m}^{T}\,u_{n}|\geq\epsilon]=N^{2}\,\max_{m,n}\mathbb{P}[|\langle
l_{m},l_{n}\rangle-\frac{1}{J}\sum_{j=1}^{J}\,u_{mj}^{T}\,u_{nj}|\geq\epsilon]$
$\displaystyle\leq$ $\displaystyle
2N^{2}\,e^{-\frac{J\epsilon^{2}}{2\xi^{2}}},$

where in the final inequality we have used Hoeffding’s Lemma for the gaussian
random variable $x:=l_{n}(\theta)\,l_{m}(\theta)$. Thus, with probability
$\geq 1-\delta$, we have:

$\max_{m,n}|\langle l_{m},l_{n}\rangle-
u_{m}^{T}\,u_{n}|\leq\sqrt{\frac{2\xi^{2}}{J}\log{\frac{2N^{2}}{\delta}}}$
(46)

The above can then be used to upper bound the difference between the true
expectation $||L-L(\tilde{w})||_{\nu,2}$ and the approximate expectation
constructed from random projections $||u-u(\tilde{w})||_{\nu,2}$, where
$\\{\tilde{w},p\\}=\\{\tilde{w}_{m-1}^{(I)},\hat{p}^{(I)}\\}$ is the output of
Algorithm 2:

$\displaystyle||L-L(\tilde{w})||^{2}_{\nu,2}-||u-u(\tilde{w})||^{2}_{\nu,2}$
$\displaystyle=$ $\displaystyle\sum_{m,n}\,[(\tilde{w}_{m}-p)(\langle
l_{m},l_{n}\rangle-u_{m}^{T}u_{n})(\tilde{w}_{n}-p)]$ (47) $\displaystyle\leq$
$\displaystyle\sum_{m,n}\,|\tilde{w}_{m}-p||\langle l_{m},l_{n}\rangle-
u_{m}^{T}u_{n}||\tilde{w}_{n}-p|$ $\displaystyle\leq$
$\displaystyle||\tilde{w}-p||_{1}^{2}\,\max_{m,n}|\langle l_{m},l_{n}\rangle-
u_{m}^{T}\,u_{n}|$

Therefore, using (46), with probability greater than $1-\delta$, one has:

$||L-L(\tilde{w})||^{2}_{\nu,2}\leq||u-u(\tilde{w})||^{2}_{\nu,2}+||\tilde{w}-p||_{1}^{2}\,\sqrt{\frac{2\xi^{2}}{J}\log{\frac{2N^{2}}{\delta}}}$
(48)

∎

## References

  * [1] V. N. Vapnik and A. Y. Chervonenkis, “On the uniform convergence of relative frequencies of events to their probabilities,” Theory of Probability and its Applications, vol. 16, no. 2, pp. 264–280, 1971.
  * [2] B. Neyshabur, R. Tomioka, and N. Srebro, “Norm-based capacity control in neural networks,” in Proceedings of The 28th Conference on Learning Theory, COLT 2015, Paris, France, July 3-6, 2015 (P. G. ”u nwald, E. Hazan, and S. Kale, eds.), vol. 40 of JMLR Workshop and Conference Proceedings, pp. 1376–1401, JMLR.org, 2015.
  * [3] P. L. Bartlett, D. J. Foster, and M. Telgarsky, “Spectrally-normalized margin bounds for neural networks,” CoRR, vol. abs / 1706.08498, 2017.
  * [4] G. K. Dziugaite and D. M. Roy, “Computing nonvacuous generalization bounds for deep (stochastic) neural networks with many more parameters than training data,” 2017.
  * [5] D. A. McAllester, “Some pac-bayesian theorems,” in Machine Learning, pp. 230–234, ACM Press, 1998.
  * [6] J. Langford and J. Shawe-Taylor, “Pac-bayes & margins,” in Advances in Neural Information Processing Systems 15 (S. Becker, S. Thrun, and K. Obermayer, eds.), pp. 439–446, MIT Press, 2003.
  * [7] M. Seeger, Bayesian Gaussian Process Models: PAC-Bayesian Generalisation Error Bounds and Sparse Approximation.  PhD thesis, University of Edinburgh, December 2003.
  * [8] O. Catoni, Chapter 1: Inductive PAC-Bayesian learning, vol. Volume 56 of Lecture Notes–Monograph Series, pp. 1–50.  Beachwood, Ohio, USA: Institute of Mathematical Statistics, 2007.
  * [9] T. Campbell and T. Broderick, “Automated scalable bayesian inference via hilbert coresets,” Journal of Machine Learning Research, vol. 20, no. 15, pp. 1–38, 2019.
  * [10] T. Campbell and T. Broderick, “Bayesian coreset construction via greedy iterative geodesic ascent,” JMLR, 2018.
  * [11] N. Golowich, A. Rakhlin, and O. Shamir, “Size-independent sample complexity of neural networks,” CoRR, vol. abs/1712.06541, 2017.
  * [12] B. Neyshabur, Z. Li, S. Bhojanapalli, Y. LeCun, and N. Srebro, “Towards understanding the role of over-parametrization in generalization of neural networks,” CoRR, vol. abs/1805.12076, 2018.
  * [13] V. Nagarajan and J. Z. Kolter, “Uniform convergence may be unable to explain generalization in deep learning,” CoRR, vol. abs/1902.04742, 2019.
  * [14] P. Derbeko, R. El-Yaniv, and R. Meir, “Explicit learning curves for transduction and application to clustering and compression algorithms,” CoRR, vol. abs/1107.0046, 2011.
  * [15] L. Bégin, P. Germain, F. Laviolette, and J.-F. Roy, “PAC-Bayesian Theory for Transductive Learning,” in Proceedings of the Seventeenth International Conference on Artificial Intelligence and Statistics (S. Kaski and J. Corander, eds.), vol. 33 of Proceedings of Machine Learning Research, (Reykjavik, Iceland), pp. 105–113, PMLR, 22–25 Apr 2014.
  * [16] V. N. Vapnik, Statistical Learning Theory.  Wiley-Interscience, 1998.
  * [17] A. Rahimi and B. Recht, “Random features for large-scale kernel machines,” in Advances in Neural Information Processing Systems 20 (J. C. Platt, D. Koller, Y. Singer, and S. T. Roweis, eds.), pp. 1177–1184, Curran Associates, Inc., 2008.
  * [18] H. Xie, M. E. Hussein, A. Galstyan, and W. Abd-Almageed, “MUSCLE: Strengthening Semi-Supervised Learning Via Concurrent Unsupervised Learning Using Mutual Information Maximization,” arXiv e-prints, p. arXiv:2012.00150, Nov. 2020.
  * [19] A. Banerjee, T. Chen, and Y. Zhou, “De-randomized pac-bayes margin bounds: Applications to non-convex and non-smooth predictors,” CoRR, vol. abs / 2002.09956, 2020.
  * [20] L. Sagun, U. Evci, V. U. Guney, Y. N. Dauphin, and L. Bottou, “Empirical analysis of the hessian of over-parametrized neural networks,” CoRR, vol. abs / 1706.04454, 2017.


# Developing a Monolithic Silicon Sensor in
a 65 nm CMOS Imaging Technology
for Future Lepton Collider Vertex Detectors

A. Simancas1,2, J. Braach3,4, E. Buschmann3, A. Chauhan1, D. Dannheim3, M. Del
Rio Viera1,2, K. Dort3,5, D. Eckstein1, F. Feindt1, I.M. Gregor1, K. Hansen1,
L. Huth1, L. Mendes1,6, B. Mulyanto1, D. Rastorguev1,7, C. Reckleben1, S. Ruiz
Daza1,2, P. Schütze1, W. Snoeys3, S. Spannagel1, M. Stanitzki1, A. Velyka1, G.
Vignola1,2, H. Wennlöf1 1Deutsches Elektronen-Synchrotron DESY, Germany
(corresponding author email: adriana.simancas@desy.de)2University of Bonn,
Germany3 Conseil Européen pour la Recherche Nucléaire CERN, Switzerland4
University of Hamburg, Germany5 University of Giessen, Germany6University of
Campinas, Brazil7University of Wuppertal, GermanyManuscript received December
9, 2022.The Tangerine project is funded by the Helmholtz Innovation Pool, 2021
- 2023.This work has been sponsored by the Wolfgang Gentner Programme of the
German Federal Ministry of Education and Research (grant no. 13E18CHA).This
project has received funding from the European Union’s Horizon 2020 Research
and Innovation programme under GA no 101004761.

###### Abstract

Monolithic CMOS sensors in a 65 nm imaging technology are being investigated
by the CERN EP Strategic R&D Programme on Technologies for Future Experiments
for an application in particle physics. The appeal of monolithic detectors
lies in the fact that both sensor volume and readout electronics are
integrated in the same silicon wafer, providing a reduction in production
effort, costs and scattering material. The Tangerine Project WP1 at DESY
participates in the Strategic R&D Programme and is focused on the development
of a monolithic active pixel sensor with a time and spatial resolution
compatible with the requirements for a future lepton collider vertex detector.
By fulfilling these requirements, the Tangerine detector is suitable as well
to be used as telescope planes for the DESY-II Test Beam facility. The project
comprises all aspects of sensor development, from the electronics engineering
and the sensor design using simulations, to laboratory and test beam
investigations of prototypes. Generic TCAD Device and Monte-Carlo simulations
are used to establish an understanding of the technology and provide important
insight into performance parameters of the sensor. Testing prototypes in
laboratory and test beam facilities allows for the characterization of their
response to different conditions. By combining results from all these studies
it is possible to optimize the sensor layout. This contribution presents
results from generic TCAD and Monte-Carlo simulations, and measurements
performed with test chips of the first sensor submission.

##  I Introduction

Monolithic CMOS sensors have found their way through imaging technologies into
High Energy Physics thanks to their good performance in particle detection.
The main characteristic is the integration of an active sensor and readout in
the same silicon wafer, which provides a reduction in production effort, costs
and scattering material.

The Tangerine Project (Towards Next Generation Silicon Detectors) [1] aims to
develop the next generation of low-material and high-precision pixel detectors
using novel technologies, targeting the requirements of a vertex detector in a
future lepton collider experiment. Work Package 1 at DESY collaborates closely
with the CERN EP Strategic R&D Programme [2] Work Package 1.2 as well as the
ALICE ITS3 Upgrade project [3] in designing, simulating and testing prototypes
in a 65 nm CMOS Imaging Technology.

In the context of particle physics applications, this contribution presents
the development of a monolithic active pixel sensor (MAPS) with a small
collection electrode in the aforementioned technology. Compared to the
previously investigated 180 nm process, the 65 nm technology offers a
significant improvement in the logic density of the pixels and a reduction of
the material thickness. The small collection electrode sensor is characterized
by a low input capacitance (in the order of fF), granting a high signal to
noise ratio (S/N) and a low power consumption.

The goal of the project is to develop a detector suitable for vertex-finding
at future lepton colliders. Therefore, the performance targets are a position
resolution of $\sim$3 µm, a time resolution of the order of nanoseconds and a
material budget of maximally 50 µm of silicon. Each pixel will provide a
charge measurement making it possible to perform charge-weighted position
measurements. Due to similar performance requirements for telescope planes at
test beam facilities, this detector will be suited to replace the current
telescope planes at the DESY-II Test Beam facility [4].

In the first foundry submission of the 65 nm technology collaboration, test
chips have been designed and produced. Joint contributions have taken place
for their characterization, as well as generic simulations validating process
modifications for enhanced sensor performance. A test chip featuring a
Krummenacher type charge sensitive amplifier designed at DESY was
characterized in [5]. The ALICE Analog Pixel Test Structure (APTS) prototype
[6] is an analog prototype designed to characterize different sensor layouts
and has been designed as part of the ALICE ITS3 upgrade efforts.

In this work, first sensor simulations as well as preliminary measurement
results of the ALICE APTS are presented.

##  II SENSOR TECHNOLOGY

Currently, three different sensor layouts are studied: standard [7], n-blanket
[8] and n-gap [9]. These designs were originally developed in a 180 nm CMOS
imaging technology to enhance depletion, timing performance and radiation
tolerance of small collection electrode MAPS.

Fig. 1 is a schematic representation of a detector with the n-gap layout. The
cross section shows half of a pixel on each side of the schematic and the
pixel edge in the center. A thin epitaxial p-doped layer is grown on a low
resistivity p-doped substrate. The p-well is the structure that hosts the in-
pixel electronics and shields them from the electric field of the active
sensor region. The n-blanket is a low dose n-doped layer implemented to create
a planar pn-junction, enlarging the depleted volume of the sensor.
Furthermore, a gap in the n-blanket produces a vertical pn-junction that
generates a lateral electric field in the farthest position from the readout
electrodes (pixel boundaries).

Figure 1: n-gap layout example of the Tangerine sensor.

To understand the 65 nm CMOS Imaging Technology and optimize the sensor
design, generic simulations and prototype testing are carried out
simultaneously. The goal of the design optimization is to customize the
electric field inside the sensor to improve efficiency and signal formation
time. The simulation cycle is of the utmost importance to reduce time and
costs invested in producing and testing prototypes.

##  III SENSOR SIMULATIONS

Since the electric fields in small collection electrode MAPS are complex,
device simulations are needed to provide insight into performance parameters
of the sensor. Studies using generic doping profiles were performed for the
three different layouts and for pitches between 10 µm and 35 µm. These
simulations include only the epitaxial layer, since no significant
contribution from the substrate to the electric field modelling is expected.
Additionally, to bias the substrate in the TCAD simulations a metallic contact
is added in the backside, as opposed to the the real detector that does not
include a backside metalization.

###  III-A Technology Computer-Aided Design (TCAD) Simulations

TCAD is a powerful tool to simulate electrical properties of semiconductors
and can be used to optimize the sensor layout and other features, such as bias
voltage configuration, to achieve the desired performance goals. TCAD contains
a finite-element simulation tool that constructs a mesh over the studied
structure and solves Poisson’s equation and carrier continuity equations to
model the electrostatic potential and other properties in each node of the
mesh. The software used for the simulations in this work is Sentaurus TCAD
from Synopsys [10].

3D quasi-stationary simulations using generic doping profiles were performed
to model the electric fields of the studied layouts, since 2D simulations
cannot resolve some of the effects introduced by the lateral electric fields.
An example of a generic 3D TCAD simulation for the n-gap design is shown in
Fig. 2. The TCAD studies carried out so far include scans over different
geometrical and operational parameters of the sensor, such as p-well opening
and bias voltage, and observing the behavior of the electric field, the
lateral electric field strength, as well as the depleted volume. The results
are reported in the following.

Figure 2: Generic 3D TCAD simulation of the n-gap layout.

The p-well opening is the distance between the edge of the collection implant
and the edge of the p-well, as shown in Fig. 1. It was varied from 1 µm to 4
µm and from the results of the standard layout it was observed that increasing
the p-well opening provided a larger depleted volume and stronger lateral
electric field. This translates into a larger signal and a faster signal
formation, but reduces the space for readout electronics. However, this effect
was much less prominent once the n-blanket modification was added.

The gap size in the n-gap layout was varied from 1 µm to 4 µm. The gap creates
a lateral electric field at the pixel boundaries, which increases with the gap
size. The gap is introduced to allow for the free charges to avoid the field
minimum reported in [9] and drift with a shorter mean free path towards the
collection electrodes, thus it improves the signal formation time.

(a) Standard

(b) n-blanket

(c) n-gap

Figure 3: Generic TCAD Simulation: current density for the three layouts with
p-well bias of -1.2 V and substrate bias of -4.8 V. Pixel pitch of 25 µm.
Brown line indicates the position of the pn-junction and white line
corresponds to the depleted volume. The high current density in the edge of
the pixels for the standard and the n-gap design indicates breakdown for these
biasing conditions.

A breakdown of the sensor was observed when fixing the p-well bias at -1.2 V,
and applying a higher bias to the substrate. The breakdown was detected early
on for the standard layout, at -2.4 V, for the n-blanket layout it was reached
at a bias of -11 V, and for the n-gap layout it was observed at -4.8V. This
behavior is in agreement with observations made during monitoring of the
leakage current in experimental measurements. Fig. 3 shows how the breakdown
is modelled in the generic TCAD simulation. The color scale corresponds to the
current density, while the white line indicates the limits of the depleted
volume. The breakdown produces a high current density in the edge of the
pixels, and the depleted volume is deformed.

The bias voltage applied to the p-well and the substrate was scanned from 0 V
to -20 V, while fixing the bias of the readout electrodes to 1.2 V. When the
bias voltage was simultaneously increased for both electrodes, an increased
depleted volume was observed. However, high values of the electric field were
detected inside the p-well structure, compromising the shielding of the
electronics. When the bias was increased only for the substrate while fixing
the p-well bias, a similar behaviour was observed for the n-blanket and n-gap
layouts, but the p-well integrity was preserved.

Fig. 4 shows the electric field obtained from TCAD simulations of the three
sensor layouts with generic doping profiles, where p-well and substrate bias
voltage were set to -4.8 V. The brown line indicates the position of the pn-
junction, the white line delimits the depleted region and the streamlines
(black arrows) indicate the instantaneous tangent to the velocity vector of
the moving charges.

Comparing between the different layouts with a pixel pitch of 25 µm, the
following can be concluded:

  * •

The standard layout (Fig. 4(a)) has a small depleted volume. The electron-hole
pairs produced outside the depleted volume will move predominantly by
diffusion in random directions and some might not reach the readout
electrodes. The expected effect on charges produced at the edge of the pixel
is a low efficiency, but a high charge-sharing between pixels, which improves
the spatial resolution.

  * •

The n-blanket layout (Fig. 4(b)) shows a larger depleted volume. The electron-
hole pairs produced in the active volume will move predominantly by drift
towards the readout electrodes. This is foreseen to produce an improvement in
efficiency, but with an impairment in spatial resolution due to lower charge-
sharing.

  * •

The n-gap layout (Fig. 4(c)) has a higher lateral electric field in the pixel
edges. The electron-hole pairs produced in the edge of the pixel will drift
with a shorter mean free path towards the readout electrodes. As a
consequence, an improvement is expected in efficiency as well as the signal
formation time, but with a further impairment in spatial resolution due to
even lower charge-sharing.

In order to quantify the effects discussed here, Monte-Carlo simulations are
required as in Section III-B.

(a) Standard

(b) n-blanket

(c) n-gap

Figure 4: Generic TCAD Simulation: electric field for the three layouts with
p-well and substrate bias of -4.8 V. Pixel pitch of 25 µm. Brown line
indicates the position of the pn-junction, black arrows represent the
streamlines and the white line delimits the depleted volume.

TCAD is also capable of simulating current pulses produced by the interaction
of a charged particle with the sensor. This is carried out with transient
simulations and can be used to estimate signal characteristics, such as signal
magnitude and time evolution.

A transient simulation for a minimum ionizing particle (MIP) traversing the
corner of a pixel with an electron/hole pair production of 63 eh/µm was
performed for the studied layouts. The result is shown in Fig. 5, confirming
that the time evolution of the signal is improved by the modifications in the
sensor, and particularly for the n-gap layout.

Figure 5: Generic TCAD Simulation: current pulse produced by a MIP traversing
the corner of a pixel, for the three different layouts.

For a standard layout sensor, a similar simulation was performed with a MIP
traversing the center of the pixel. By integrating the induced charge for the
duration of the signal, a total of $\sim$750 electrons was obtained.
Furthermore, generic TCAD transient simulations have provided valuable
feedback for the ASIC design, such as expected signal magnitudes to define a
reasonable threshold in the readout electronics.

###  III-B Allpix2 Simulations

TCAD transient simulations are time consuming and have a high computational
budget, thus quasi-stationary simulations from TCAD are combined with Monte-
Carlo simulations to obtain high statistics data and calculate the performance
parameters of the sensor [11]. Within the Tangerine project, this combination
of simulations is used to quantify the effects reported in Section III-A.

Allpix2 [12] is a modular framework developed for Monte-Carlo simulations of
semiconductor radiation detectors. It provides the possibility to build a
matrix of pixels by replicating the TCAD electric fields simulated on a single
cell.

First results from Monte-Carlo simulations using generic TCAD fields from this
work have been reported in [1]. The performance parameters of interest are
detection efficiency, cluster size, and spatial resolution. Results confirm
that the modifications to the sensor layout are valuable as they increase the
efficient operating margin of the sensor. The observed trends are equivalent
for all tested bias voltages.

##  IV PROTOTYPE TESTING

The purpose of the ALICE APTS [6] prototype is to characterize different
sensor designs. Several of these chips are being studied and the results will
allow for direct comparisons with simulations to be performed. This section
presents the activities within the Tangerine project of system integration and
testing of the ALICE APTS in laboratories and test beam campaigns with
preliminary results.

###  IV-A APTS and Data-Acquisition (DAQ) System

The ALICE APTS contains a matrix of 4 × 4 square pixels and has been produced
in the three studied sensor designs. The tested devices have a pixel pitch of
25 µm, the pixels are DC coupled to the front-end electronics and each pixel
contains a source follower as buffered analog output. The latter makes the
APTS a structure suitable for comparisons against simulations.

The Caribou System [13, 14] is used as data-acquisition system for the ALICE
APTS. It is an open-source set of hardware, firmware and software for
laboratory and beam tests. The modular hardware consists of three boards: the
application-specific chip board containing the sensor, the periphery CaR board
which provides current and voltage sources together with a physical interface
between System-on-Chip (SoC) and the detector, and the evaluation board that
contains the Xilinx Zynq SoC, which runs the detector control and the data-
processing firmware.

For the ALICE APTS, a custom chip board was designed and produced with
amplification and signal shaping at each pixel output. The readout was
performed with two 8-channel 65 MS/s ADCs on the CaR board, using a custom
firmware block. The firmware provides a configurable trigger logic, where
either an external or an internal trigger can be selected with programmable
thresholds and adjustable latency. The firmware is compatible with the AIDA
Trigger Logic Unit [15] (TLU).

The data-acquisition framework employed in the project is EUDAQ2 [16]. It is a
generic data-acquisition software for use in conjunction with beam telescopes
at charged particle beam lines. It allows for storage and synchronization of
data from several systems. A detector-specific decoder is used to interpret
the raw data for further analysis.

The recorded data with this DAQ system consists of waveforms sampled at 65 MHz
with a configurable number of samples. Internal triggers were used for
laboratory tests, while for the test beam data acquisition an external trigger
was provided by the TLU. Fig. 6 shows a waveform example produced by a MIP in
an APTS with the n-blanket layout. The shape of the waveform is dominated by
the amplifiers employed in the chip board of the DAQ system.

Figure 6: Waveform of a MIP measured with the ALICE APTS in the n-blanket
layout at -3.6 V p-well and substrate bias.

###  IV-B Laboratory Characterization

The activities performed in the laboratory involved the optimization of the
front-end operation parameters and studies with charge injection and
radioactive sources for gain calibration.

Since the acquired signals are represented in ADC units, a calibration is
required to relate these values to number of electrons. This way, it is
possible to match the thresholds of data and simulation and quantify the
agreement between them. The calibration of the tested devices was performed
with X-ray fluorescence, an 55Fe radioactive source and test pulses.

The decay of 55Fe produces two X-ray emissions that are considered
monochromatic: K-alpha of 5.9 keV and K-beta of 6.5 keV. For the interaction
of these X-rays with silicon, a production of approximately 1640 and 1800
electron-hole pairs, respectively, is expected. The calibration uses the
K-alpha and K-beta peaks shown in the 55Fe spectrum in Fig. 7, which was
measured with an APTS in the n-blanket layout. The result is a combination of
the individual pixel charge spectra, so the tail at the left of the spectrum
is a combination of Compton scattering, charge-sharing and threshold effects
that differ per pixel. The calibration was validated by applying the
calibration on the spectrum from the X-ray fluorescence of titanium.

Figure 7: Calibrated spectrum of an 55Fe source acquired with an ALICE APTS in
the n-blanket layout at -3.6 V p-well and substrate bias.

###  IV-C test beam Measurements

The ALICE APTS was tested in the DESY-II Test Beam Facility [4] to
characterize its performance with MIPs. A MIMOSA26 telescope [17] was used as
a reference system to reconstruct the individual particle tracks of the beam.
The APTS acted as a Device Under Test (DUT) and was placed in the middle of
the 6 telescope planes, orthogonal to the beam. Additionally, a TelePix [18]
plane was used as trigger in the data acquisition. It provides time stamps
with a precision below 5 ns, and a fast digital hit output signal which allows
for triggering in a configurable region of interest. This plane was placed
downstream with respect to the telescope. With this setup it is possible to
determine the effect of particles impinging on different positions in the
evaluated prototype and obtain the performance parameters mentioned in Section
III-B.

For the data analysis, a modular framework called Corryvreckan [19] is used.
It allows for online monitoring and offline event building in complex data-
taking environments combining detectors with different readout architectures.
With this software, the particle tracks are reconstructed from the information
provided by the telescope planes. One particle interaction can produce a
signal in different pixels. The one that registers the higher signal is called
seed pixel, and together with the surrounding pixels they constitute a
cluster. After correlating the tracks with the DUT information, clusters in
the DUT are associated to the reconstructed tracks and it is possible to
perform efficiency and resolution studies.

The devices were tested under different operational conditions. Given the
maximum beam rate of 5 kHz, a beam cross section of the order of centimeters
and the active area of the 25 µm pitch ALICE APTS of 100 × 100 µm, it was
necessary to record data for several hours for each investigated setting. It
was discovered that, due to the long acquisition time and temperature changes
in the test beam area, the relative position of the DUT with respect to the
reference was changing. This required to develop a new alignment method that
corrected the relative drift within the same acquisition run. This method has
shown promising results, but it is still under evaluation. For this reason,
the quantitative results are still preliminary, but it has been possible to
observe encouraging qualitative results that are reported in the following.

After applying the calibration to test beam data of an n-blanket design, the
charge distribution of the seed pixels was obtained as shown in Fig. 8. The
distribution corresponds approximately to a Landau distribution convolved with
a Gaussian. A most probable value of around 600 electrons is observed. The
peak at the end of the distribution is due to saturation of the ADCs.

Figure 8: Charge distribution of seed pixels for associated clusters acquired
with an ALICE APTS in the n-blanket layout at -3.6 V p-well and substrate
bias.

The cluster size depends on the charge-sharing between pixel cells and the
chosen thresholds. The detectors with sensor modifications (n-blanket and
n-gap layouts) exhibited a lower cluster size in comparison to the standard
layout, indicating a similar trend to the observations made in simulations.

The detection efficiency of the detector was measured by relating the particle
hits in the DUT to the reconstructed tracks. As expected from simulations and
previous works [8, 20], an improvement in efficiency was observed for the
n-blanket and n-gap layouts.

##  V Conclusion

The Tangerine project participates in the investigation of monolithic pixel
detectors in 65 nm CMOS Imaging Technology lead by the CERN EP Strategic R&D
Programme on Technologies for Future Experiments.

A monolithic active pixel sensor in 65 nm CMOS imaging technology is being
developed within the Tangerine project. Simulations and prototype testing are
complementary activities carried out within the project and the collaboration
institutes.

Device simulations using generic doping profiles provide valuable insight for
sensor optimization and understanding sensor behaviour for novel technologies.
The combination with Monte-Carlo simulations produces results which can be
directly compared with experimental tests.

Generic TCAD simulations have confirmed that the principles to modify the
sensor in the 180 nm and 65 nm are generally applicable. Generic simulations
have allowed for establishing the parameters of relevance for the sensor
optimization and their effect on the operation of the detector, such as p-well
opening, n-gap size and sensor bias voltage. The results show an overall
agreement with what has been observed in other technologies investigated for
MAPS. Efficiency and resolution studies using generic TCAD simulations
combined with Monte-Carlo simulations are ongoing.

The ALICE APTS has been evaluated in laboratory and test beams. A custom DAQ
system based on Caribou was developed and integrated for these measurements.
From the experimental activities, the functionality of the test setup has been
demonstrated. Preliminary results on charge distribution, cluster size and
detection efficiency have shown a qualitative agreement with simulations. More
detailed studies including spatial resolution and timing are continuing.

Further test campaigns are planned for the near future, as well as dedicated
simulations on timing performance. A fully integrated chip with a larger pixel
matrix, designed jointly by CERN, DESY and IFAE, has been submitted to the
foundry. This will allow for recording high-statistic data samples and for
further improvement of the comparison with simulations.

## Acknowledgments

The measurements presented have been performed at the test beam Facility at
DESY Hamburg (Germany), a member of the Helmholtz Association (HGF).

The authors wish to express their gratitude to the CERN EP R&D WP 1.2 and
especially to the designers of the APTS and the ALICE ITS3 measurement team
for their support.

(c) All figures and pictures by the author(s) under a CC BY 4.0 license,
unless otherwise stated.

## References

  * [1] H. Wennlöf, A. Chauhan, M. Del Rio Viera, D. Eckstein, F. Feindt, I. M. Gregor, K. Hansen, L. Huth, L. Mendes, B. Mulyanto, D. Rastorguev, C. Reckleben, S. Ruiz Daza, P. Schütze, A. Simancas, S. Spannagel, M. Stanitzki, A. Velyka, and G. Vignola, “The Tangerine project: Development of high-resolution 65 nm silicon MAPS,” _Nuclear Instruments and Methods in Physics Research, Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ , vol. 1039, no. March, p. 167025, 2022. [Online]. Available: https://doi.org/10.1016/j.nima.2022.167025
  * [2] CERN Collaboration, “Strategic R&D Programme on Technologies for Future Experiments - Annual report 2020,” _CERN Document Server_ , no. April, 2020\. [Online]. Available: http://cds.cern.ch/record/2764386?ln=en
  * [3] ALICE Collaboration, “Letter of Intent for an ALICE ITS Upgrade in LS3,” 2019\.
  * [4] R. Diener, J. Dreyling-Eschweiler, H. Ehrlichmann, I. M. Gregor, U. Kötz, U. Krämer, N. Meyners, N. Potylitsina-Kube, A. Schütz, P. Schütze, and M. Stanitzki, “The DESY II test beam facility,” _Nuclear Instruments and Methods in Physics Research, Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ , vol. 922, no. November 2018, pp. 265–286, 2019. [Online]. Available: https://doi.org/10.1016/j.nima.2018.11.133
  * [5] A. Chauhan, M. Del Rio Viera, D. Eckstein, F. Feindt, I.-M. Gregor, K. Hansen, L. Huth, L. Mendes, B. Mulyanto, D. Rastorguev, C. Reckleben, S. R. Daza, P. Schütze, A. Simancas, S. Spannagel, M. Stanitzki, A. Velyka, G. Vignola, and H. Wennlöf, “Towards a new generation of Monolithic Active Pixel Sensors,” _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ , vol. 1047, p. 167821, February 2023. [Online]. Available: https://doi.org/10.1016/j.nima.2022.167821
  * [6] W. Deng, G. A. Rinella, M. Aresti, J. Baudot, F. Benotto, S. Beole, W. Bialas, G. Borghello, S. Bugiel, M. Campbell, F. Carnesecchi, E. Charbon, L. Cecconi, A. D. Martin, A. Dorokhov, C. Ferrero, J. Hasenbichler, G. H. Hong, C. Hu, G. Huang, A. Kluge, T. Kugathasan, P. L. Rocca, M. Mager, P. Martinengo, M. Munker, S. Perciballi, F. Piro, F. Prino, K. Rebane, I. Sanna, V. Sarritzu, U. Savino, W. Snoeys, M. Suljic, G. Termo, A. Triffiro, A. Turcato, G. Usai, Y. Wu, and S. Korea, “Design of an analog monolithic pixel sensor prototype in TPSCo 65 nm CMOS imaging technology,” _Submitted to Journal of Instrumentation_ , 2022.
  * [7] S. Senyukov, J. Baudot, A. Besson, G. Claus, L. Cousin, A. Dorokhov, W. Dulinski, M. Goffe, C. Hu-Guo, and M. Winter, “Charged particle detection performances of CMOS pixel sensors produced in a 0.18µm process with a high resistivity epitaxial layer,” _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ , vol. 730, no. 2013, pp. 115–118, December 2013. [Online]. Available: http://dx.doi.org/10.1016/j.nima.2013.03.017
  * [8] W. Snoeys, G. Aglieri Rinella, H. Hillemanns, T. Kugathasan, M. Mager, L. Musa, P. Riedler, F. Reidt, J. Van Hoorne, A. Fenigstein, and T. Leitner, “A process modification for CMOS monolithic active pixel sensors for enhanced depletion, timing performance and radiation tolerance,” _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ , vol. 871, no. June, pp. 90–96, November 2017. [Online]. Available: http://dx.doi.org/10.1016/j.nima.2017.07.046
  * [9] M. Munker, M. Benoit, D. Dannheim, A. Fenigstein, T. Kugathasan, T. Leitner, H. Pernegger, P. Riedler, and W. Snoeys, “Simulations of CMOS pixel sensors with a small collection electrode, improved for a faster charge collection and increased radiation tolerance,” _Journal of Instrumentation_ , vol. 14, no. 5, 2019.
  * [10] Synopsys, “TCAD - Technology Computer Aided Design (TCAD).” [Online]. Available: https://www.synopsys.com/silicon/tcad.html
  * [11] D. Dannheim, K. Dort, D. Hynds, M. Munker, A. Nürnberg, W. Snoeys, and S. Spannagel, “Combining TCAD and Monte Carlo methods to simulate CMOS pixel sensors with a small collection electrode using the Allpix2 framework,” _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ , vol. 964, no. December 2019, p. 163784, June 2020. [Online]. Available: https://doi.org/10.1016/j.nima.2020.163784
  * [12] S. Spannagel, K. Wolters, D. Hynds, N. Alipour Tehrani, M. Benoit, D. Dannheim, N. Gauvin, A. Nürnberg, P. Schütze, and M. Vicente, “Allpix2: A modular simulation framework for silicon detectors,” _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ , vol. 901, no. June, pp. 164–172, September 2018. [Online]. Available: https://doi.org/10.1016/j.nima.2018.06.020
  * [13] H. Liu, M. Benoit, H. Chen, K. Chen, F. A. Di Bello, G. Iacobucci, F. Lanni, I. Peric, B. Ristic, M. V. B. Pinto, W. Wu, L. Xu, and G. Jin, “Development of a modular test system for the silicon sensor R&D of the ATLAS Upgrade,” _Journal of Instrumentation_ , vol. 12, no. 1, 2017.
  * [14] T. Vanat, “Caribou – A versatile data acquisition system,” in _Proceedings of Topical Workshop on Electronics for Particle Physics — PoS(TWEPP2019)_ , vol. 370, no. September. Trieste, Italy: Sissa Medialab, March 2020, p. 100. [Online]. Available: https://pos.sissa.it/370/100
  * [15] P. Baesso, D. Cussans, and J. Goldstein, “The AIDA-2020 TLU: A flexible trigger logic unit for test beam facilities,” _Journal of Instrumentation_ , vol. 14, no. 9, 2019.
  * [16] Y. Liu, M. Amjad, P. Baesso, D. Cussans, J. Dreyling-Eschweiler, R. Ete, I. Gregor, L. Huth, A. Irles, H. Jansen, K. Krueger, J. Kvasnicka, R. Peschke, E. Rossi, A. Rummler, F. Sefkow, M. Stanitzki, M. Wing, and M. Wu, “EUDAQ2—A flexible data acquisition software framework for common test beams,” _Journal of Instrumentation_ , vol. 14, no. 10, pp. P10 033–P10 033, oct 2019. [Online]. Available: https://iopscience.iop.org/article/10.1088/1748-0221/14/10/P10033
  * [17] H. Jansen, S. Spannagel, J. Behr, A. Bulgheroni, G. Claus, E. Corrin, D. Cussans, J. Dreyling-Eschweiler, D. Eckstein, T. Eichhorn, M. Goffe, I. M. Gregor, D. Haas, C. Muhl, H. Perrey, R. Peschke, P. Roloff, I. Rubinskiy, and M. Winter, “Performance of the EUDET-type beam telescopes,” _EPJ Techniques and Instrumentation_ , vol. 3, no. 1, pp. 1–18, 2016.
  * [18] H. Augustin, R. Diener, S. Dittmeier, P. M. Freeman, J. Hammerich, A. Herkert, L. Huth, D. Immig, U. Krämer, N. Meyners, I. Perić, O. Schäfer, A. Schöning, A. Simancas, M. Stanitzki, D. Stuart, and B. Weinläder, “Upgrading the beam telescopes at the DESY II Test Beam Facility,” _Nuclear Instruments and Methods in Physics Research, Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ , vol. 1040, no. July, p. 167183, 2022. [Online]. Available: https://doi.org/10.1016/j.nima.2022.167183
  * [19] D. Dannheim, K. Dort, L. Huth, D. Hynds, I. Kremastiotis, J. Kröger, M. Munker, F. Pitters, P. Schütze, S. Spannagel, T. Vanat, and M. Williams, “Corryvreckan: a modular 4D track reconstruction and analysis software for test beam data,” _Journal of Instrumentation_ , vol. 16, no. 03, p. P03008, March 2021. [Online]. Available: https://dx.doi.org/10.1088/1748-0221/16/03/P03008
  * [20] M. Dyndal, V. Dao, P. Allport, I. A. Tortajada, M. Barbero, S. Bhat, D. Bortoletto, I. Berdalovic, C. Bespin, C. Buttar, I. Caicedo, R. Cardella, F. Dachs, Y. Degerli, H. Denizli, L. F. S. de Acedo, P. Freeman, L. Gonella, A. Habib, T. Hemperek, T. Hirono, B. Hiti, T. Kugathasan, I. Mandić, D. Maneuski, M. Mikuž, K. Moustakas, M. Munker, K. Oyulmaz, P. Pangaud, H. Pernegger, F. Piro, P. Riedler, H. Sandaker, E. Schioppa, P. Schwemling, A. Sharma, L. S. Argemi, C. S. Sanchez, W. Snoeys, T. Suligoj, T. Wang, N. Wermes, and S. Worm, “Mini-MALTA: radiation hard pixel designs for small-electrode monolithic CMOS sensors for the High Luminosity LHC,” _Journal of Instrumentation_ , vol. 15, no. 02, pp. P02 005–P02 005, February 2020. [Online]. Available: https://iopscience.iop.org/article/10.1088/1748-0221/15/02/P02005