# Multi-layered Cepstrum for Instantaneous Frequency Estimation

###### Abstract

We propose the multi-layered cepstrum (MLC) method to estimate multiple
fundamental frequencies (MF0) of a signal under challenging contamination such
as high-pass filter noise. Taking the operation of cepstrum (i.e., Fourier
transform, filtering, and nonlinear activation) recursively, MLC is shown as
an efficient method to enhance MF0 saliency in a step-by-step manner.
Evaluation on a real-world polyphonic music dataset under both normal and low-
fidelity conditions demonstrates the potential of MLC.

Index Terms—  Cepstrum, multiple fundamental frequency estimation, time-
frequency analysis.

##  1 Introduction

Multiple fundamental frequency (MF0) estimation is an essential problem in
analyzing various multi-component signals including multi-talker speech [1,
2], polyphonic music [3], fetal electrocardiogram (FECG),
photoplethysmographic (PPG) signals [4], and IoT sensor data [5], to name but
a few. Numerous methods have been proposed to enhance the saliency of F0s in
the signal representation, such as the autocorrelation function (ACF) [6],
generalized cepstrum [7], nonnegative matrix factorization (NMF), MUltiple
SIgnal Classification (MUSIC) [8], combined frequency and periodicity (CFP)
[9, 10, 11], and deep learning approaches [12]. The generalized cepstrum,
being one of the earliest proposed MF0 estimation methods [7], still gains
attention in recent years, by showing its effectiveness in suppressing the
unwanted harmonics and in localizing the F0 components by combining it with a
spectrum [9, 10, 13]. A recent study further proved that a cepstrum can
perfectly localize the instantaneous F0s of a signal with the intrinsic mode
function (IMF) type [11], a general type of non-stationary oscillatory
signals.

One issue which is relatively less discussed in MF0 estimation is the so-
called _convolutional_ noise in signals [14]. Convolutional noise is
ubiquitous; it occurs when signals are transmitted through an intermediate
channel or device. An MF0 algorithm usually fails in the presence of
convolutional noise, for example, in the processing of audio signals played by
a smartphone, as the response of a smartphone speaker typically behaves like a
high-pass filter whose cutoff frequency is around 500 Hz [15, 9]. Such a high-
pass filter suppresses the true F0 peaks, as similar to the case of the
“missing fundamental” effects in low-pitched musical signals [16]. In such
cases, to identify the true F0 peaks merely from the interleaved harmonics of
multiple components is challenging.

In this paper, we generalize the notion of cepstrum to resolve the above-
mentioned issue in MF0 estimation. The proposed _multi-layered cepstrum_ (MLC)
refines the saliency of the F0s in a recursive manner. In a nutshell, the MLC
takes the operation of cepstrum (i.e., Fourier transform, filtering and
nonlinear activation) repeatedly: for a time-domain input signal
$\mathbf{x}\in\mathbb{R}^{N}$, $l\geq 0$, the $l^{\textrm{th}}$-layer output
of MLC is:

$\mathbf{z}^{(0)}=\sigma^{(0)}\left(|\mathbf{F}\mathbf{x}|\right)\,,\quad\mathbf{z}^{(l)}=\sigma^{(l)}\left(\mathbf{W}^{(l)}\mathbf{F}\mathbf{z}^{(l-1)}\right)\,,$

where $\mathbf{F}\in\mathbb{C}^{N\times N}$ is the $N$-point discrete Fourier
transform (DFT) matrix, $\mathbf{W}\in\mathbb{R}^{N\times N}$ is a high-pass
filter, and $\sigma(\cdot)$ is an element-wise power-scaled rectification unit
(see Section 3). In the following sections, we will explain how the idea of
MLC is motivated, how it works in MF0 estimation, and its structural
resemblance to multi-layered perceptrons (MLP). We will then focus on
combining two MLCs lying in both the frequency and time domains to achieve MF0
estimation by employing the CFP approach. Experiments on both synthetic and
real data demonstrates the effectiveness of MLC in MF0 estimation in both
normal and contaminated conditions, and indicate positively the fact that MLC
performs better than other methods based on a ‘shallow’ DFT.

##  2 Motivation of MLC

The _cepstrum_ $\hat{\mathbf{z}}:=\mathbf{F}^{-1}\sigma(\mathbf{F}\mathbf{x})$
of $\mathbf{x}$ [17, 18] is a classic method to estimate the F0 of
$\mathbf{x}$. A logarithm cepstrum employing $\sigma(x):=\log x$ is commonly
used in single F0 estimation, while the _generalized cepstrum_ (GC) employing
$\sigma(x):=x^{\gamma}$, $0<\gamma\leq 2$, is found more applicable in MF0
estimation, probably because it is numerically more stable than using the
logarithm function and is more robust to noise [19, 20]. Cepstrum-based MF0
estimation is based on an assumption that the fast-varying components in a
spectrum (e.g., harmonic sequences) are relevant to F0s, while the slow-
varying ones (e.g., spectral envelope) are irrelevant to F0s. Saliency of F0s
can be enhanced in the _quefrency_ 111Since the quefrency has the same unit as
time, in this paper, the terms quefrency and time are used interchangeably.
domain by taking another Fourier transform and a high-pass filter with a
cutoff quefrency at $q_{c}$ on a nonlinearly-stretched spectrum, where the
nonlinear operation is to fit the magnitude of the spectrum to humans’
perception scale.

To motivate the MLC, we revise the above-mentioned assumption to that 1) the
fast-varying and _periodic_ components in a spectrum are relevant to F0s while
others are not, and that 2) such periodic components still remain periodic
(and therefore fast-varying) after a Fourier transform, which is a basic
property of DFT. These assumptions allow us to consider repeating the
operation of cepstrum on $\mathbf{x}$ to iteratively discard the slow-varying
or aperiodic components and to purify the F0-relevant components, which should
remain after a number of filtering processes. This idea has been preliminarily
verified in [21], where the _generalized cepstrum of spectrum_ (GCoS), being
the $2^{\text{nd}}$layer output of the MLC, behaves as a more succinct feature
indicating the F0 peaks than a spectrum. This paper is an extension of [21] to
a multi-layered time-frequency analysis tool.

Admittedly, the above-mentioned idea is also motivated by deep learning, which
catches wide attention recently in MF0 estimation [22, 12, 23, 24, 25]. Not
only achieving state-of-the-art results, deep learning also provides inspiring
observation when solving the MF0 estimation problem. For example, results in
[12] show that the learned networks for melody tracking resemble traditional
pitch detection methods: the first layer behaves like a spectral analyzer,
while the second layer behaves like a comb filter. It is again intriguing to
consider a multi-layered network for MF0 estimation, where the well-known
pitch detection blocks such as the Fourier transform and filters are taken
directly as subnetworks.

##  3 The MLC-CFP method

We formulate the MLC of a non-stationary signal $\mathbf{x}$ in the time-
frequency or time-quefrency domain. Given a window function
$\mathbf{h}\in\mathbb{R}^{N}$ and hop size $H$, the short-time Fourier
transform (STFT) of $\mathbf{x}$, $\mathbf{X}:=\mathbf{X}[k,n]$, is
represented as
$\mathbf{X}:=\mathbf{X}[k,n]=\sum^{N-1}_{m=0}\mathbf{x}[m+nH]\mathbf{h}[m]e^{-\frac{j2\pi
km}{N}}$, where $\mathbf{X}\in\mathbb{C}^{N\times M}$, $M$ is the number of
frames, $k$ is the frequency index, and $n$ is the time index. Denote
$|\mathbf{X}|$ the element-wise absolute value of $\mathbf{X}$. The
$l^{\textrm{th}}$-layer short-time MLC is:

$\mathbf{Z}^{(0)}:=\sigma^{(0)}\left(\left|\mathbf{X}\right|\right)\,,\quad\mathbf{Z}^{(l)}:=\sigma^{(l)}\left(\mathbf{W}^{(l)}\mathbf{F}\mathbf{Z}^{(l-1)}\right)\,,$
(1)

where $\mathbf{Z}^{(l)}\in\mathbb{R}^{N\times M}$,
$\mathbf{Z}^{(l)}:=\mathbf{Z}^{(l)}[k,n]$ for even $l$,
$\mathbf{Z}^{(l)}:=\mathbf{Z}^{(l)}[q,n]$ for odd $l$, and $q$ is the
quefrency index. $\sigma^{(l)}$ is a nonlinear, element-wise power function
such that for $\gamma_{l}>0$, $\sigma^{(l)}\left(x\right)=x^{\gamma_{l}}$ for
$x>0$, and $\sigma^{(l)}\left(x\right)=0$ for $x\leq 0$. With the $N$-point
DFT matrix $\mathbf{F}$, the frequency resolution of $\mathbf{Z}^{(l)}[k,n]$
is $f_{s}/N$ per bin, and the time resolution of $\mathbf{Z}^{(l)}[q,n]$ is
then $1/f_{s}$ per bin. $\mathbf{W}^{(i)}\in\mathbb{R}^{N\times N}$ represents
a high-pass filter; it filters out the elements lying below a cutoff frequency
(quefrency) $i_{c}$. More precisely, $\mathbf{W}^{(l)}$ is a diagonal matrix
such that $\mathbf{W}^{(l)}[i,i]=1$ when $i>i_{c}$ and $i<N-i_{c}$, while
$\mathbf{W}^{(l)}[i,i]=0$ otherwise. For even $l$, $i_{c}:=k_{c}$ defines the
cutoff frequency at $f_{c}=k_{c}f_{s}/N$, and for odd $l$, $i_{c}:=n_{c}$
defines the cutoff quefrency at $q_{c}=n_{c}/f_{s}$. Also note that since
$\mathbf{Z}^{(l)}[:,n]$ is even for $l>2$, we use $\mathbf{F}$ instead of
$\mathbf{F}^{-1}$ when transforming from the frequency domain to the time
domain. Many of the basic pitch detection functions are essentially special
cases of (1). For example, $\mathbf{Z}^{(0)}$ is the STFT of $\mathbf{x}$, and
$\mathbf{Z}^{(1)}$ is known as the autocorrelation function (ACF) for
$\gamma_{0}=2$ and $\gamma_{1}=1$. (1) also resembles a DNN structure with
zero bias vector: $\mathbf{W}^{(l)}\mathbf{F}$ resembles the fully-connected
layer and $\sigma^{(l)}$ resembles the activation function.

To enhance the true F0 components and suppress the high-order harmonics,
consider the last two output layers, $\mathbf{Z}^{(l_{e})}$ and
$\mathbf{Z}^{(l_{o})}$, of an MLC, where $l_{e}$ is even and $l_{o}$ is odd.
$\mathbf{Z}^{(l_{e})}$ is in the frequency domain and $\mathbf{Z}^{(l_{o})}$
is in the quefrency domain. Based on the CFP approach, we utilize the duality
between them to enhance the saliency of F0 by employing a nonlinear mapping to
$\mathbf{Z}^{(l_{o})}$ from the quefrency domain to frequency domain, and then
multiplying it by $\mathbf{Z}^{(l_{e})}$ [13, 9, 10, 11]:

$\mathbf{Y}^{(l_{e},l_{o})}[k,n]=\mathbf{Z}^{(l_{e})}[k,n]\mathbf{Z}^{(l_{o})}\left[\lfloor
N/k\rceil,n\right]$ (2)

where the time-quefrency representation $\mathbf{Z}^{(l_{o})}[k,n]$ is
nonlinearly mapped into the frequency domain, and $\lceil\cdot\rfloor$ is the
rounding function. Since $\mathbf{Z}^{(l_{o})}$ suppresses the harmonic peaks
in $\mathbf{Z}^{(l_{e})}$ while $\mathbf{Z}^{(l_{e})}$ suppresses the sub-
harmonic peaks in $\mathbf{Z}^{(l_{o})}$, the resulting CFP representation
$\mathbf{Y}^{(l_{e},l_{o})}:=\mathbf{Y}^{(l_{e},l_{o})}[k,n]$ has only true F0
peaks enhanced. The nonlinear mapping of $\mathbf{Z}^{(l_{o})}$ in (2) can be
implemented by a filterbank following a custom-defined indexes of $k$; see
[10, 21] for details. We name the above process as the MLC-CFP algorithm
hereafter.

Fig. 1: Simulation result. Upper: $\mathbf{x}$; lower: $\hat{\mathbf{x}}$.
From left to right: $\mathbf{Z}^{(0)}$ (power-scale spectrogram),
$\mathbf{Z}^{(1)}$ (GC), $\mathbf{Z}^{(2)}$ (GCoS), and $\mathbf{Y}^{(2,1)}$
(CFP representation). The window size is 8 seconds and the hop size is 1
second.

##  4 Numerical results

We illustrate the effectiveness of the MLC-CFP method in emphasizing F0 peaks
of a simulated signal contaminated with both additive and convolutional noise.
Given two components, $\mathbf{x}_{1}$ and $\mathbf{x}_{2}$, both are sampled
at $f_{s}=1$ kHz. $\mathbf{x}_{1}$ is a square wave having 20% duty cycle, and
has a constant F0 at $2$ Hz. $\mathbf{x}_{2}$ is a frequency-modulated
sawtooth wave with positive ramp and vertical drop, and has time-varying F0
represented as $f_{0}[n]=2.5+\cos(2\pi n/10f_{s})$. To simulate the effect of
convolutional noise, $\mathbf{x}_{1}$ is passed through a 10th-order
Butterworth low-pass filter, and $\mathbf{x}_{2}$ is passed through a 10th-
order Butterworth high-pass filter. Both filters have a cutoff frequency at 10
Hz. We consider the 2-layer MLC of two sample signals: the first one is
$\mathbf{x}=\mathbf{x}_{1}+\mathbf{x}_{2}$, which is with convolutional noise
only, while the second one,
$\tilde{\mathbf{x}}=\mathbf{x}_{1}+\mathbf{x}_{2}+\mathbf{n}+\mathbf{d}$, is
further contaminated with additive pink noise and impulse noise, where
$\mathbf{n}$ is additive pink noise making the SNR of $\tilde{\mathbf{x}}$ be
10 dB and $\mathbf{d}$ is an impulse at the 80th second.

Fig. 1 shows the STFT ($\mathbf{Z}^{(0)}$), GC ($\mathbf{Z}^{(1)}$), GCoS
($\mathbf{Z}^{(2)}$), and CFP representation ($\mathbf{Y}^{(2,1)}$) of
$\mathbf{x}$ and $\hat{\mathbf{x}}$, respectively. The parameters
$[\gamma_{0},\gamma_{1},\gamma_{2}]=[0.24,0.6,1]$ are used. For the STFT of
both samples, the F0 peaks and the low-order harmonics of $\mathbf{x}_{2}$ are
all unseen, and only the first 5 multiples of $\mathbf{x}_{1}$ are visible due
to low-pass filtering. To capture the true F0 of $\mathbf{x}_{2}$ from the
STFT is very challenging. However, such information can be found in the GC:
although high-pass filtering noise attenuates the true F0 component in the
spectrum, it does not change the period of the signal. Such a phenomenon can
also be seen in the case with pink noise. The true F0 trajectory then appears
in the GCoS with high resolution; this is because the nonlinear scaling of the
GC provides the mechanism to _down-mixing_ the frequency when it is
transformed back to the frequency domain. Interestingly, we can see that
although the true F0 of $\mathbf{x}_{2}$ is missed, the resulting F0
trajectory of $\mathbf{x}_{2}$ has resolution even higher than the one of
$\mathbf{x}_{1}$. Since most of the harmonic peaks of $\mathbf{x}_{1}$ are
eliminated, its spectrum cannot support better resolution when transformed to
the quefrency domain. Finally, the CFP representation computed through (2)
clearly indicates the true F0 trajectories of $\mathbf{x}_{1}$ and
$\mathbf{x}_{2}$ with less interference.

##  5 Real-data Experiments

Table 1: Comparison of different layer combinations and testing method on the Bach10 dataset. Layers | Brute Force | Greedy | SGD (10-fold CV)
---|---|---|---
$\\{L,L+1\\}$ | $\gamma_{0}$ | $\gamma_{1}$ | $\gamma_{2}$ | $\gamma_{3}$ | $\gamma_{4}$ | $\gamma_{5}$ | P (%) | R (%) | F (%) | $\gamma_{i}$ | P (%) | R (%) | F (%) | P (%) | R (%) | F (%)
0th & 1st | 0.3 | 1 | x | x | x | x | 77.72 | 83.51 | 80.51 | 0.33 | 78.05 | 83.57 | 80.71 | 77.67 | 78.43 | 78.05
1st & 2nd | 0.3 | 0.5 | 1 | x | x | x | 80.55 | 87.20 | 83.74 | 0.48 | 81.05 | 85.23 | 83.09 | 79.35 | 88.10 | 83.50
2nd & 3rd | 0.2 | 0.6 | 0.9 | 1 | x | x | 79.33 | 89.93 | 84.30 | 0.71 | 80.28 | 86.66 | 83.35 | 82.61 | 85.99 | 84.26
3rd & 4th | 0.1 | 0.9 | 0.9 | 0.5 | 1 | x | 81.06 | 90.80 | 85.66 | 0.65 | 80.16 | 87.08 | 83.47 | 81.02 | 85.83 | 83.36
4th & 5th | 0.1 | 0.9 | 0.9 | 0.7 | 0.8 | 1 | 81.14 | 90.67 | 85.64 | 0.78 | 80.20 | 87.27 | 83.58 | 81.80 | 86.22 | 83.95
5th & 6th | 0.1 | 0.9 | 0.9 | 0.7 | 0.8 | 0.5 | 82.35 | 91.10 | 86.50 | 0.70 | 80.29 | 87.57 | 83.77 | 81.97 | 82.18 | 82.07

We apply the MLC-CFP method for MF0 estimation of polyphonic music data, and
investigate the effects of 1) the number of layers, 2) critical parameters
$\gamma_{i}$ of each layer, and 3) the presence of convolutional noise. We
used the Bach10 dataset, containing ten quartets of four different instruments
[26] for the evaluation of multi-pitch estimation (MPE). To see the behavior
of the algorithm in a straightforward way, we consider three methods to
optimize parameters $\\{\gamma_{i}\\}^{L}_{i=0}$ for an $L$-layer MLC-CFP
algorithm combining $\mathbf{Z}^{(L-1)}$ and $\mathbf{Z}^{(L)}$:

1) the _brute-force method_ that searches for the optimal parameters
exhaustively for each $\gamma_{i}$. The search range of every $\gamma_{i}$ is
from 0.1 to 0.9 with step size of 0.1.

2) the _greedy method_ that searches for the optimal value of $\gamma_{l}$
based on the optimal values of $\\{\gamma_{i}\\}^{l-1}_{i=0}$. The search
range is from 0.01 to 0.99 with step size of 0.01, since its light computation
allows us to use finer grids. For the post-processing of the brute-force and
the greedy methods, we use the CFP-based algorithm reported in [21], as a
simplified version of [9] without too detailed hand-crafted rules.

3) the _stochastic gradient descent (SGD) method_ that minimizes the binary
cross-entropy between $\mathbf{Y}^{(L)}$ (with 88 log-frequency filterbanks
and a sigmoid function as the output layer) and the ground-truth piano roll by
performing SGD over every $\gamma_{i}$. We apply non-repeated 10-fold cross-
validation on the dataset: for each time one music piece is chosen as the test
data and others as training data. The resulting F-scores are obtained by
summing all the counts of the true positives, false positives, and false
negatives. We initialize $\gamma_{0}$ to 0.24, $\gamma_{1}$ to 0.6, and all
other $\gamma_{i}$s to 1, and use the SGD optimizer with 0.1 learning rate, a
mini-batch size of 256, and train the network up to 40 epochs.

For other parameters, we compute 7939-point DFT using Blackman-Harris window,
the hop size is 10 ms, the cutoff frequency $f_{c}$ is 27.5 Hz (frequency of
A0), and the cutoff quefrency $q_{c}$ is 0.24 ms (period of C8). Since there
are many task-dependent parameters, all the source codes will be announce
after this paper is accepted for reproducibility.

Table 1 lists the optimal precision (P), recall (R), F-score (F), along with
$\gamma_{i}$ for every pairs
$\left[\mathbf{Z}^{(L-1)},\mathbf{Z}^{(L)}\right]$ for all optimization
approaches. When using the brute-force method, the resulting P, R and F all
increase consistently when the number of layers increases. Note that the
recall increases dramatically with higher-order layer because the weak F0
peaks are enhanced, similar to the simulation example in Fig. 1. When using
the greedy method, we observe that although the performance is sub-optimal, it
still increases as $L$ increases. When using gradient descent, there is no
consistent improvement in F-scores related to the number of layers, possibly
because the hand-crafted rules such as harmonic selection [9] could not be
applied to this model. However, it still gives competitive performance for
$L=3$. More specifically, when compared to other MF0 estimation methods, all
the three methods in our discussion outperform the state-of-the-art matrix
factorization approaches such as constrained NMF and probabilistic latent
component analysis (PLCA) (see Table III in [9] for more details). The brute-
force method also outperform the CFP result reported in [9], in which the
parameters for the hand-crafted rules are also based on brute-force searching.

Fig. 2 shows the resulting F-scores on the Bach10 dataset with high-pass
degradation, which is simulated by a 4th-order Butterworth high-pass filters
with cutoff frequency from 10 to 1000 Hz. The parameters used here are those
found from the brute force method; see Table 1. Fig. 2 clearly shows that with
more layers the results are better: even under a severe cutoff frequency at 1
kHz, a 6-layer MLC still gives 72.74% F-score, 25% higher than the one of a
1-layer MLC. Even more interestingly, results in some degradation conditions
are better than those in the normal condition. For example, pairs
$\mathbf{Z}^{(5)}$ and $\mathbf{Z}^{(6)}$ achieve 87.10% F-score with cutoff
frequency at 100 Hz; this is probably because the simulated degradation
suppresses the low-frequency noise incidentally.

To summarize, our experiments show that 1) in the MLC-CFP algorithm, a deep
DFT do performs better than a shallow DFT, and 2) MLC is highly robust to
high-pass noise.

$10^{1}$$10^{2}$$10^{3}$$0.5$$0.6$$0.7$$0.8$$0.9$Cutoff frequency
(Hz)F-score$\mathbf{Z}^{(0)}$ & $\mathbf{Z}^{(1)}$$\mathbf{Z}^{(1)}$ &
$\mathbf{Z}^{(2)}$$\mathbf{Z}^{(2)}$ & $\mathbf{Z}^{(3)}$$\mathbf{Z}^{(3)}$ &
$\mathbf{Z}^{(4)}$$\mathbf{Z}^{(4)}$ & $\mathbf{Z}^{(5)}$$\mathbf{Z}^{(5)}$ &
$\mathbf{Z}^{(6)}$ Fig. 2: F-scores of Bach10 dataset under high-pass noise.

##  6 Discussion and conclusion

We have introduced the MLC, explained its physical insight in signal
processing, and verified its potential in detecting multiple F0s through
simulation and real-world data under normal and noisy conditions. We also
emphasize that the MLC is not only useful in the signals contaminated with
high-pass convolutional noise, but also useful in all signals which inherently
have weak F0 magnitude. For example, the F0 of bassoon in the Bach10 dataset
is usually weaker than other instruments (violin, clarinet and saxophone), but
experiments show that this issue is addressed by the MLC. From the result in
the filter test we also suggest that in practice using MLC with a little high-
pass filter can get better result.

One risk of MLC is that it may generate fake F0 peaks. For example, two F0s at
$a$ Hz and $b$ Hz may produces a cross term at $|a-b|$ Hz due to effect of
nonlinear activation. This cross term cannot be suppressed and may be even
enhanced because it appears in both the time-domain and frequency-domain
representations. To restrict the strength of cross terms, the results indicate
a possible strategy: setting a small $\gamma_{0}$ and larger $\gamma_{i}$s at
the higher layers. Our pilot study also shows that even fake cross terms are
produced, they are mostly weak unless all $\gamma_{i}$s approach zero.

Developing a more efficient way to optimize the parameters efficiently on a
large-scale dataset with complicated sources and harmonic structures will be
another focus in the future. From the current finding, we summarize that MLC-
CFP works well in general and is relatively easy to optimize for $L\leq 3$,
but a great potential exists in a deeper MLC.

## References

  * [1] G. Pirker, M. Wohlmayr, S. Petrik, and F. Pernkopf,  “A pitch tracking corpus with evaluation on multipitch tracking scenario,”  in INTERSPEECH, 2011.
  * [2] M. Wu, D.L. Wang, and G. J. Brown,  “A multipitch tracking algorithm for noisy speech,”  IEEE Trans. Speech Audio Process., vol. 11, no. 3, pp. 229–241, 2003.
  * [3] E. Benetos, S. Dixon, D. Giannoulis, H. Kirchhoff, and A. Klapuri,  “Automatic music transcription: challenges and future directions,”  J. Intelligent Information Systems, vol. 41, no. 3, pp. 407–434, 2013.
  * [4] Z. Zhang, Z. Pi, and B. Liu,  “TROIKA: A general framework for heart rate monitoring using wrist-type photoplethysmographic signals during intensive physical exercise,”  IEEE Trans. Bio-Med. Eng., vol. 62, no. 2, pp. 522–531, 2015.
  * [5] Z. Li, J. Wang, and J. Han,  “ePeriodicity: Mining event periodicity from incomplete observations,”  IEEE Trans. Knowledge and Data Engineering, vol. 27, no. 5, pp. 1219–1232, 2015.
  * [6] L. Rabiner,  “On the use of autocorrelation analysis for pitch detection,”  IEEE Trans Acoust., Speech, Signal Process., vol. 25, no. 1, pp. 24–33, 1977.
  * [7] T. Kobayashi and S. Imai,  “Spectral analysis using generalized cepstrum,”  IEEE Trans. Acoust., Speech, Signal Process, vol. 32, no. 5, pp. 1087–1089, 1984.
  * [8] M. Græsbøll Christensen, P. Stoica, A. Jakobsson, and S. Holdt Jensen,  “Multi-pitch estimation,”  Signal Processing, vol. 88, no. 4, pp. 972–983, 2008.
  * [9] L. Su and Y.-H. Yang,  “Combining spectral and temporal representations for multipitch estimation of polyphonic music,”  IEEE/ACM Trans. Audio, Speech, Language Process., vol. 23, no. 10, pp. 1600–1612, 2015.
  * [10] L. Su, T.-Y. Chuang, and Y.-H. Yang,  “Exploiting frequency, periodicity and harmonicity using advanced time-frequency concentration techniques for multipitch estimation of choir and symphony,”  in Proc. ISMIR, 2016, pp. 393–399.
  * [11] C.-Y. Lin, L. Su, and H.-T. Wu,  “Wave-shape function analysis–when cepstrum meets time-frequency analysis,”  Journal of Fourier Analysis and Applications, pp. 1–55, 2017.
  * [12] P. Verma and R. W. Schafer,  “Frequency estimation from waveforms using multi-layered neural networks,”  in INTERSPEECH, 2016, pp. 2165–2169.
  * [13] G. Peeters,  “Music pitch representation by periodicity measures based on combined temporal and spectral representations,”  in Proc. ICASSP, 2006.
  * [14] H. Hermansky, N. Morgan, and H.-G. Hirsch,  “Recognition of speech in additive and convolutional noise based on rasta spectral processing,”  in Proc. ICASSP. IEEE, 1993, vol. 2, pp. 83–86.
  * [15] M. Mauch and S. Ewert,  “The audio degradation toolbox and its application to robustness evaluation.,”  in ISMIR, 2013, pp. 83–88,  https://code.soundsoftware.ac.uk/projects/audio-degradation-toolbox (v0.2).
  * [16] T. D. Rossing, F. R. Moore, and P. A. Wheeler,  The Science of Sound,  Addison Wesley San Francisco, 2002.
  * [17] B.P. Bogert, M.J.R. Healy, and J.W. Tukey,  “The quefrency analysis of time series for echoes: Cepstrum, pseudo-autocovariance, cross-cepstrum and shape cracking,”  in Time Series Analysis, 1963, vol. 15, pp. 209–243.
  * [18] A.V. Oppenheim and R.W. Schafer,  Discrete-Time Signal Processing,  Prentice Hall, 3rd edition, 2009.
  * [19] T. Tolonen and M. Karjalainen,  “A computationally efficient multipitch analysis model,”  IEEE Speech Audio Process., vol. 8, no. 6, pp. 708–716, 2000.
  * [20] A. Klapuri,  “Multipitch analysis of polyphonic music and speech signals using an auditory model,”  IEEE Trans. Audio, Speech, Language Process., vol. 16, no. 2, pp. 255–266, 2008.
  * [21] L. Su,  “Between homomorphic signal processing and deep neural networks: Constructing deep algorithms for polyphonic music transcription,”  in APSIPA ASC, 2017.
  * [22] S. Sigtia, E. Benetos, and S. Dixon,  “An end-to-end neural network for polyphonic piano music transcription,”  IEEE/ACM Trans. Audio, Speech, Language Process., vol. 24, no. 5, pp. 927–939, 2016.
  * [23] R. M. Bittner, B. McFee, J. Salamon, P. Li, and J. P. Bello,  “Deep salience representations for $f_{0}$ estimation in polyphonic music,”  in Proc. ISMIR, Oct. 2017.
  * [24] J. Thickstun, Z. Harchaoui, and S. Kakade,  “Learning features of music from scratch,”  arXiv preprint arXiv:1611.09827, 2016.
  * [25] A. Elowsson,  “Polyphonic pitch tracking with deep layered learning,”  arXiv preprint arXiv:1804.02918, 2018.
  * [26] Z. Duan, B. Pardo, and C. Zhang,  “Multiple fundamental frequency estimation by modeling spectral peaks and non-peak regions,”  IEEE/ACM Trans. Audio, Speech, Language Process., vol. 18, no. 8, pp. 2121–2133, 2010.


such a mixture, and the results we give in this section should extend to that
setting as well.

One starting point from (García Trillos et al., 2021) is the concrete
description that they give for a mixture model to be “well-separated” in an
appropriate sense.

######  Definition 8 (García Trillos et al. (2021))

Consider a probability density

$\rho(x)=\sum_{k=1}^{K}w_{k}\rho_{k}(x)$

on $\mathbb{R}^{d}$ with $w_{k}\geq 0$ and $\rho_{k}$ being probability
densities. We call such a density a _mixture model_ , and each of the
$\rho_{k}$ the mixture components. We then define the following parameters:

  1. 1.

Overlapping. The overlapping, $\mathcal{S}$, of a mixture model is defined to
be

$\mathcal{S}:=\max_{i\neq j}\int\frac{\rho_{i}\rho_{j}}{\rho}\,dx,$

  2. 2.

Coupling. The coupling, $\mathcal{C}$, of a mixture model is defined to be

$\mathcal{C}:=\max_{k}\int\left|\frac{\nabla\rho_{k}}{\rho_{k}}-\frac{\nabla\rho}{\rho}\right|\rho_{k}\,dx,$

  3. 3.

Indivisibility. The indivisibility, $\Theta$, of the mixture model is defined
to be

$\Theta:=\min_{k}\min_{u\perp 1,\int u^{2}\rho_{k}=1}\int|\nabla
u|^{2}\rho_{k}.$

In García Trillos et al. (2021) the authors’ definition also applies to
probability distributions defined on manifolds. We choose not to state our
results in terms of manifolds as it doesn’t align with other portions of this
work, but all of the results in this subsection would apply in the manifold
setting as well. In that work one of the main assumptions (Assumption 7 in
García Trillos et al. (2021)) is that

$\begin{cases}&\rho,\rho_{k}\in C^{1}\\\
&\Delta_{\rho},\Delta_{\rho_{k}}\text{ have discrete spectrum }\\\
&L^{2}(\rho),L^{2}(\rho_{k})\text{ each have an orthogonal
eigenbasis.}\end{cases}$ (S)

We will take this as a standing assumption throughout this section.

In the previous definition, the overlapping $\mathcal{S}$ describes the amount
to which pairs of the mixtures have coincident densities, relative to the
underlying density. In a case where the mixtures have disjoint supports, this
parameter will be zero. The indivisibility $\Theta$ is a measure of how
difficult it is, in terms of Dirichlet energy, to divide any one of the
clusters: this parameter would be large in the case of strongly log-concave
mixture components. The coupling parameter $\mathcal{C}$ is somewhat more
subtle, and measures a type of relative entropy between $\rho$ and $\rho_{k}$:
again this parameter would be zero for disjoint mixture components.

García Trillos et al. (2021) then define a well-clustered mixture model via
the relationships

$\mathcal{S}\ll 1\quad\text{and}\quad\frac{\mathcal{C}}{\Theta}\ll 1.$ (11)

These conditions, of course, require some knowledge of the underlying
distributions, and may be difficult to verify for a particular data set.
However, most natural models of well-clustered data would satisfy these
assumptions. For example, in the case where we have $K$ clusters with disjoint
support we will have $\mathcal{S}$ and $\mathcal{C}$ both equal to zero. They
would also hold for Gaussian mixtures which are sufficiently separated (and
$\mathcal{S}$ and $\mathcal{C}$ could be quantified in terms of the means and
variances). In this sense, we view these conditions as applicable to many
different models of separated data components.

With these definitions in hand, we now state our main exploration result,
which guarantees that clusters are explored efficiently: this proposition
largely turns out to be a direct consequence of the estimates given by García
Trillos et al. (2021).

###### Proposition 9

Consider a mixture model with parameters $(\mathcal{S},\mathcal{C},\Theta)$
which satisfies Assumption (S). Let $\mathcal{K}_{t}(x^{*},x)$ be the heat
kernel with parameter $t$, namely the solution to the heat equation satisfying
$\mathcal{K}_{t}(x^{*},x)=e^{t\Delta_{\rho}}\delta_{x^{*}}$ with Neumann
boundary conditions. Consider parameters $a,b,\sigma$ as in Proposition 11.
Then, under the assumptions on parameters stated in Proposition 11, the kernel
$\mathcal{K}_{t}$ will be a $(\delta,\zeta,\varepsilon)$ class separator with

$\delta=\tilde{\delta}+a^{2},\quad\zeta=\frac{(1-\sin(b))^{2}}{w_{max}}\cos(2(\sigma+b))-\xi,\quad\varepsilon=\frac{(1+\sin(b))^{2}}{w_{min}}\cos(\pi/2-2(\sigma+b))+\xi,$

where

$\xi:=|e^{\lambda_{K}t}-1|\frac{(1+\sin(b))^{2}}{w_{min}}\|\rho\|_{\infty}+\min_{t_{1},s>0,t_{1}+s<t}e^{\lambda_{K+1}(t-t_{1}-s)}\|e^{t_{1}\Delta}\|_{L^{1}(\rho)\to
L^{2}(\rho)}\|e^{s\Delta}\|_{L^{2}(\rho)\to L^{\infty}(\rho)}$

Proof  We consider a sequence of eigenvalue, eigenfunction pairs
$(\lambda_{k},e_{k})$ of the operator $\Delta_{\rho}$, normalized in the
$\rho$-weighted $L^{2}$ norm. We then define
$g(y)=e^{t_{1}\Delta}\delta_{x^{*}}-\sum_{k=1}^{K}e^{\lambda_{k}t_{1}}e_{k}(x^{*})\rho(x^{*})e_{k}(y)$,
with $t_{1}$ being a parameter that we choose later. We can interpret $g$ as
the projection of the solution to the heat equation onto the “high” Fourier
modes. We find that

$\mathcal{K}_{t}(x^{*},x)=\sum_{k=1}^{K}e_{k}(x^{*})\rho(x^{*})e_{k}(x)+\sum_{k=1}^{K}(e^{\lambda_{k}t}-1)e_{k}(x^{*})\rho(x^{*})e_{k}(x)+e^{s\Delta}\sum_{k=K+1}^{\infty}e^{\lambda_{k}(t-t_{1}-s)}e_{k}(x)\langle
e_{k},g\rangle_{\rho}.$ (12)

We can then bound

$\left|\sum_{k=1}^{K}(e^{\lambda_{k}t}-1)e_{k}(x)e_{k}(x^{*})\rho(x^{*})\right|\leq|e^{\lambda_{K}t}-1|\max_{k=1\dots
K}\|e_{k}\|_{L^{\infty}(E_{k})}^{2}\|\rho\|_{\infty},$

where we will only need to estimate this term in regions (i.e., a family of
disjoint sets $E_{k}$ that are defined in Proposition 11, for the purposes of
proving $(\delta,\zeta,\varepsilon)$ class separation) where the
eigenfunctions $e_{k}$ are not too large. We can also bound

$\displaystyle\left|e^{s\Delta}\sum_{k=K+1}^{\infty}e^{\lambda_{k}(t-t_{1})}e_{k}(x)\langle
e_{k},g\rangle_{\rho}\right|$
$\displaystyle\leq\|e^{s\Delta}\|_{L^{2}(\rho)\to
L^{\infty}(\rho)}\left\|\sum_{k=K+1}^{\infty}e^{\lambda_{k}(t-t_{1})}e_{k}(x)\langle
e_{k},g\rangle\right\|_{L^{2}(\rho)}$ $\displaystyle\leq
e^{\lambda_{K+1}(t-t_{1})}\|g\|_{L^{2}(\rho)}\|e^{s\Delta}\|_{L^{2}(\rho)\to
L^{\infty}(\rho)}$ $\displaystyle\leq
e^{\lambda_{K+1}(t-t_{1}-s)}\|e^{t_{1}\Delta}\|_{L^{1}(\rho)\to
L^{2}(\rho)}\|e^{s\Delta}\|_{L^{2}(\rho)\to L^{\infty}(\rho)}.$

Hence we define

$\xi:=\min_{t_{1},s>0,t_{1}+s<t}|e^{\lambda_{K}t}-1|\max_{k=1\dots
K}\|e_{k}\|_{L^{\infty}(E_{k})}^{2}\|\rho\|_{\infty}+e^{\lambda_{K+1}(t-t_{1}-s)}\|e^{t_{1}\Delta}\|_{L^{1}(\rho)\to
L^{2}(\rho)}\|e^{s\Delta}\|_{L^{2}(\rho)\to L^{\infty}(\rho)}.$

In turn, we only need to estimate the first term in Equation (12). This is
considered in Proposition 11. That Proposition also provides an upper bound
upon $\|e_{k}\|_{L^{\infty}(E_{k})}^{2}$, which completes the proof.

The quantity $\xi$ represents all of the “high Fourier mode” effects encoded
by the heat kernel. In the case of well-separated mixtures, one would expect
that $\lambda_{K}\sim 0$ and $\lambda_{K+1}$ would be much larger, and hence
one would anticipate that it whould be possible to choose time parameters so
that $\xi$ is small. This “spectral gap” is quantified in the following
proposition from (García Trillos et al., 2021).

######  Proposition 10 (Propositions 27 and 38, García Trillos et al. (2021))

Consider a mixture model with parameters $(\mathcal{S},\mathcal{C},\Theta)$
which satisfies Assumption (S), and suppose that $K\mathcal{S}<1$. Then the
following bounds hold:

$\displaystyle\left(\sqrt{\Theta(1-K\mathcal{S})}-\frac{\sqrt{\mathcal{C}K\mathcal{S}}}{1-\mathcal{S}}\right)^{2}\leq\lambda_{K+1}$
$\displaystyle\lambda_{K}\leq\frac{K\mathcal{C}}{1-K\mathcal{S}^{1/2}}.$

In turn, the main consideration in bounding $\xi$ relates to the choice of $t$
and the operator norm bounds on the heat kernel. The types of bounds that are
implicitly assumed in the definition of $\xi$ are called _ultracontractivity_
estimates for the heat kernel, and are known to be finite under various
classes of assumptions: a standard reference on the subject is (Davies, 1989).
In particular, these bounds will hold for densities which are compactly
supported and do not degenerate, or for mixtures whose densities are smooth
and (asymptotically) log-concave. While for empirical data it will be
challenging to verify bounds on $\xi$, there are many important models,
including those used previously in the active learning literature (Karzand and
Nowak, 2020; Cloninger and Mhaskar, 2021), where it is possible to bound $\xi$
(e.g., Gaussian mixtures, uniform distributions on sets with very small
overlap). Furthermore, if we restrict our “hypothesis class” to be densities
for which the heat kernel obeys desired bounds (for example by assuming
uniform log-concavity on the mixture components), then we can completely
control all of the remaining terms in Proposition 9 using the parameters
$\mathcal{S},\mathcal{C},\Theta$.

Now, we turn to a key estimate for the heat kernel, regarding the effect of
the low Fourier modes after neglecting decay. For convenience, and following
the notation in (García Trillos et al., 2021), we define the spectral
embedding $F$ to be the mapping from $\mathbb{R}^{d}$ to $\mathbb{R}^{K}$
given by

$F(x):=(e_{i}(x))_{i=1}^{K}.$

###### Proposition 11

Consider a mixture model with parameters $(\mathcal{S},\mathcal{C},\Theta)$.
We define the parameters

$\Xi:=K\frac{(\Lambda-\sqrt{\mathcal{S}})^{2}}{4}+4K^{3/2}\left(\frac{1}{\sqrt{1-K\Lambda}}-1\right),$

with

$\Lambda:=4\left(\sqrt{\frac{\Theta(1-K\mathcal{S})}{\mathcal{C}}-\frac{\sqrt{K\mathcal{S}}}{(1-\mathcal{S})}}\right)^{-1}+\sqrt{\mathcal{S}}.$

We assume that

$\displaystyle a,b>0$ (13)
$\displaystyle\frac{a\sin(b)}{\sqrt{w_{max}}}\geq\sqrt{\Xi}$ (14)
$\displaystyle\mathcal{S}<\frac{w_{min}(1-\cos^{2}(\sigma))}{w_{max}\cos^{2}(\sigma)K^{2}}$
(15) $\displaystyle\Lambda-\sqrt{\mathcal{S}}>0$ (16) $\displaystyle\Lambda
K<1$ (17) $\displaystyle b+\sigma<\frac{\pi}{4}$ (18)

and we let

$\tilde{\delta}:=\frac{w_{max}K^{2}\mathcal{S}}{w_{min}^{2}(1-\cos^{2}(\sigma))}+\frac{1-\cos^{2}(\sigma)}{w_{min}}$

Then there exists disjoint sets $E_{k}$ so that
$\rho_{k}(E_{k})>1-(\tilde{\delta}+a^{2})$, which satisfy the inequalities,
for $x\in E_{j},x^{*}\in E_{k}$

$\displaystyle\left|\sum_{k=1}^{K}e_{k}(x)e_{k}(x^{*})\right|$
$\displaystyle=|\langle
F(x),F(x^{*})\rangle|\leq\frac{(1+\sin(b))^{2}}{w_{min}}\cos\left(\frac{\pi}{2}-2(\sigma+b)\right)\qquad\text{
if }j\neq k$ (19)
$\displaystyle\left|\sum_{k=1}^{K}e_{k}(x)e_{k}(x^{*})\right|$
$\displaystyle=|\langle
F(x),F(x^{*})\rangle|\geq\frac{(1-\sin(b))^{2}}{w_{max}}\cos(2(\sigma+b))\qquad\text{
if }j=k$ (20) $\displaystyle\left|\sum_{k=1}^{K}e_{k}(x)e_{k}(x^{*})\right|$
$\displaystyle=|\langle
F(x),F(x^{*})\rangle|\leq\frac{(1+\sin(b))^{2}}{w_{min}}\qquad\text{ if }j=k$
(21)

Proof  The ingredients for this proof are all present in the proof of Theorem
10 in García Trillos et al. (2021). However, their statements are all made in
terms of the overall probability of points which are in the support of a set
of approximately orthogonal cones: in other words, they estimate probabilities
in terms of $\rho$. This was natural in their unsupervised setting, whereas in
ours we are closer to the supervised setting and care about the labels
associated with each mixture component; that is, we want to estimate the size
of sets in terms of the $\rho_{k}$. We simply sketch how the estimates we use
can be obtained by small modifications of their proofs.

In one of the first main steps in their proof, they identify the functions
$\tilde{q}_{k}:=\sqrt{\frac{\rho_{k}}{\rho}}$ as a system of approximate
eigenfunctions, which are almost orthogonal in a way that can be quantified by
$\mathcal{S}$. Specifically, in the proof of Proposition 29, given an angle
$0<\sigma<\frac{\pi}{4}$ they define the disjoint sets

$A_{k}:=\left\\{x:\rho_{k}(x)>\cos^{2}(\sigma)\sum_{j=1}^{K}\rho_{j}(x)\right\\}$

They then prove that

$\rho\left(\bigcup_{k=1}^{K}A_{k}\right)\geq
1-\frac{w_{max}K^{2}\mathcal{S}}{w_{min}(1-\cos^{2}(\sigma))}.$

We now extend their estimate to apply to the conditional probabilities
$\rho_{k}$, as needed in estimating the class separator parameters in
Definition 4. We first notice that the previous inequality immediately implies
that

$\rho_{k}\left(\left(\bigcup_{j=1}^{K}A_{j}\right)^{c}\right)\leq\frac{w_{max}K^{2}\mathcal{S}}{w_{min}^{2}(1-\cos^{2}(\sigma))}.$

On the other hand, by the definition of the $A_{k}$ we have that for $j\neq k$

$w_{j}\rho_{k}(A_{j})\leq(1-\cos^{2}(\sigma))w_{j}\rho_{j}(A_{j})\leq(1-\cos^{2}(\sigma))\rho(A_{j}),$

which after summing over $j\neq k$ gives

$w_{min}\rho_{k}\left(\bigcup_{j\neq k}A_{j}\right)\leq(1-\cos^{2}(\sigma)).$

This then implies that

$\rho_{k}(A_{k}^{c})\leq\frac{w_{max}K^{2}\mathcal{S}}{w_{min}^{2}(1-\cos^{2}(\sigma))}+\frac{1-\cos^{2}(\sigma)}{w_{min}}=:\tilde{\delta}.$

In terms of embeddings, they define the mapping
$F^{Q}:=(\tilde{q}_{k})_{k=1}^{K}$. The proof of Theorem 10 in their paper
shows that there exists an orthonormal matrix $O$ so that

$W_{2}^{2}(OF_{\sharp}\rho,F_{\sharp}^{Q}\rho)\leq
K\frac{(\Lambda-\sqrt{\mathcal{S}})^{2}}{4}+4K^{3/2}\left(\frac{1}{\sqrt{1-K\Lambda}}-1\right)=:\Xi,$

with

$\Lambda:=4\left(\sqrt{\frac{\Theta(1-K\mathcal{S})}{\mathcal{C}}-\frac{\sqrt{K\mathcal{S}}}{(1-\mathcal{S})}}\right)^{-1}+\sqrt{\mathcal{S}}$

Here we have used the notation $F_{\sharp}\rho$ to denote the standard push-
forward measure, which is defined to be the measure so that
$F_{\sharp}\rho(B):=\rho(F^{-1}(B))$ for Borel sets $B$.

Using the definition of the Wasserstein distance, we can also infer that

$W_{2}^{2}(OF_{\sharp}\rho_{k},F_{\sharp}^{Q}\rho_{k})\leq\frac{\Xi}{w_{min}}.$

We also notice that

$F^{Q}(A_{k})\subset\left\\{z:\frac{z_{k}}{|z|}>\cos(\sigma),\frac{1}{\sqrt{w_{max}}}\leq|z|\leq\frac{1}{\sqrt{w_{min}}}\right\\}.$

In words, the embedding $F^{Q}$ separates the $A_{k}$ into orthogonal cones
with angle $\sigma$ and places the mass in $A_{k}$ at least distance
$\frac{1}{\sqrt{w_{max}}}$ from the origin.

Finally, in Proposition 22 of their paper they show stability bounds of cones
under the Wasserstein distance. In particular, given a vector $v\neq 0$ and an
angle parameter $0<\sigma<\frac{\pi}{4}$, if we let the set $C$ be a cone
defined by

$C:=\left\\{x:\frac{\langle v,x\rangle}{|v||x|}\geq\cos(\sigma)\right\\}$

then we can define the set $C_{r}:=C\cap(B(0,r))^{c}$. By following the proof
of Proposition 22 by García Trillos et al. (2021), as applied to one cone as
opposed to an orthogonal system of cones, one has that, for a probability
measure $\mu_{1}$ with $\mu_{1}(C_{r})>1-\delta$, then for any $a,b>0$
satisfying

$ra\sin(b)\geq W_{2}(\mu_{1},\mu_{2})$

and assuming that $0<\sigma+b<\pi/4$, we will also have that
$\mu_{2}(\tilde{C}\cap(B(0,\tilde{r}))^{c})>1-\hat{\delta}$, with $\tilde{C}$
being the cone centered on $v$ with angle smaller than $\sigma+b$,
$\tilde{r}=r(1-\sin(b))$ and $\hat{\delta}=\delta+a^{2}$. We furthermore note
that their proof could be extended so that if we have
$\mu_{1}(C\cap(B(0,R)\setminus B(0,r)))$ with $R>r$ then we will have
$\mu_{2}(\tilde{C}\cap(B(0,\tilde{R})\setminus B(0,\tilde{r})))$ with
$\tilde{R}=R(1+\sin(b))$.

Putting these facts together, we then have that if our parameters satisfy
equations (14)-(18) and if we define a cone centered on the $k$-th coordinate
vector with angle $\sigma+b$ to be $C_{k}$, then for
$\tilde{r}=\frac{1-\sin(b)}{\sqrt{w_{max}}}$ and
$\tilde{R}=\frac{1+\sin(b)}{\sqrt{w_{min}}}$ and letting
$\tilde{C}_{k}=C_{k}\cap(B(0,\tilde{R})\setminus B(0,\tilde{r}))$ we have

$OF_{\sharp}\rho_{k}(\tilde{C}_{k})\geq 1-(\tilde{\delta}+a^{2}).$

We notice that if $z\in\tilde{C}_{k}$ and $z^{\prime}\in\tilde{C}_{j}$ with
$j\neq k$ then we have that

$\langle z,z^{\prime}\rangle\leq\tilde{R}^{2}\cos(\pi/2-2(\sigma+b)).$

On the other hand, if $j=k$ we have that

$\tilde{r}^{2}\cos(2(\sigma+b))\leq\langle
z,z^{\prime}\rangle\leq\tilde{R}^{2}.$

This then concludes the proof.

The parameters in the previous propositions are somewhat involved, and for
convenience we provide the following simplified result.

###### Corollary 12

Suppose that $\mathcal{C}\leq\kappa,\mathcal{S}\leq\kappa$ and that
$\Theta\geq\sqrt{\kappa}$. Let $\|e^{\Delta}\|_{L^{1}(\rho)\to
L^{\infty}(\rho)}\leq\sqrt{M}$. Then for $\kappa$ sufficiently small and for
$t=-\log(\kappa)\kappa^{1/2}$ we have that the heat kernel $\mathcal{K}_{t}$
is a $(\delta,\zeta,\varepsilon)$ class separator with $\delta\leq
C\kappa^{1/4},\zeta\geq w_{max}^{-1}-C(|\log(\kappa)|\kappa^{1/2}+M\kappa)$
and $\varepsilon\leq C(\kappa^{1/4}+M\kappa)$, where $C$ is a constant that
depends only upon $K,w$ and $\|\rho\|_{\infty}$.

Proof  This simply amounts to choosing the parameters $a,b,\sigma$
appropriately, and then bounding all of the terms. Throughout this proof, we
let $C$ be a constant that varies line by line, and may depend upon
$K,w_{min},w_{max},\|\rho\|_{\infty}$, but does not depend upon the other
parameters. First, we notice that

$\left(\sqrt{\frac{\Theta(1-K\mathcal{S})}{\mathcal{C}}-\frac{\sqrt{K\mathcal{S}}}{(1-\mathcal{S})}}\right)^{-1}\leq
C\kappa^{1/4}$

In turn

$\Xi\leq C\kappa^{1/2}$

Now let $a=\kappa^{1/4}$ and $b=\kappa^{1/4}$, and $\sigma=\kappa^{1/4}$. Then
$\tilde{\delta}\leq C\kappa^{1/2}$.

We also note that using our assumptions we have that $|\lambda_{K}|\leq\kappa$
and $|\lambda_{K+1}|\geq C\kappa^{1/2}$. Therefore by setting
$t=-\log(\kappa)\kappa^{-1/2}$ we have that $|e^{\lambda_{K}t}-1|\leq
C|\log(\kappa)|\kappa^{1/2}$ and $|e^{\lambda_{K+1}(t-2)}|\leq C\kappa$. In
turn this implies that $\xi\leq C(|\log(\kappa)|\kappa^{1/2}+M\kappa)$ Putting
this all together, then gives that our kernel is a
$(\delta,\zeta,\varepsilon)$ class separator with $\delta\leq C\kappa^{1/2}$,
$\zeta\geq w_{max}^{-1}-C(|\log(\kappa)|\kappa^{1/2}+M\kappa)$ and
$\varepsilon\leq C(\kappa^{1/4}+M\kappa)$.

###  6.3 Derivation of asymptotic exploitation

In this section, we identify a formal continuum limit associated with
Dirichlet Learning, and identify its steady states. For simplicity, we focus
our attention on the case where $\mathcal{X}$ is a bounded open set in
$\mathbb{R}^{d}$, and that the underlying features are associated with a
smooth and bounded density $\rho(x)=\sum_{k=1}^{K}w_{k}\rho_{k}(x)$ supported
on $\mathcal{X}$ with smooth and bounded class-conditional densities
$\rho_{k}(x)$ for $k\in[K]$. Recalling (8), we furthermore define
$V(\alpha(x)):=\mathcal{A}_{var}(x)$ and assume that our acquisition function
samples proportional to $e^{\lambda V(\alpha(x))}$ (i.e., the Proportional
Sampling policy introduced in 4). Furthermore, we will assume that the value
of $\lambda$ is allowed to vary during the sampling process.

In this context, we can view the expected change in the pseudo-labels from the
$\ell$-th to $(\ell+1)$-th observation by the relation

$\mathbb{E}_{X_{\ell+1},Y_{\ell+1}}[\alpha^{\ell+1}_{k}(x)-\alpha^{\ell}_{k}(x)|(X_{1},Y_{1}),\ldots,(X_{\ell},Y_{\ell})]=\int_{\mathcal{X}}\frac{w_{k}\rho_{k}(z)}{\rho(z)}\mathcal{K}(z,x)\frac{e^{\lambda(\ell)V(\alpha^{\ell}(z))}}{\int_{\mathcal{X}}e^{\lambda(\ell)V(\alpha^{\ell}(\tilde{z}))}\,d\tilde{z}}\,dz,$

where $\alpha_{k}^{\ell}(x)\geq 0$ is the $k$-th entry of the concentration
parameter vector $\alpha^{\ell}(x)\in\mathbb{R}_{+}^{K}$ at the $\ell$-th
observation.

By assuming that we take a large number of discrete steps in one unit of
“continuous time”, the law of large numbers then formally leads to the
evolution equation

$\displaystyle\partial_{t}\alpha(x,t)$
$\displaystyle=\int_{\mathcal{X}}\frac{w_{k}\rho_{k}(z)}{\rho(z)}\mathcal{K}(z,x)\frac{e^{\lambda(t)V(\alpha(z,t))}}{\int_{\mathcal{X}}e^{\lambda(t)V(\alpha(\tilde{z},t))}\,d\tilde{z}}\,dz$
$\displaystyle=:\int_{\mathcal{X}}\eta^{K}(z)\mathcal{K}(z,x)q(z,t)\,dz,$ (22)

where we have introduced the notation $\eta^{K}(x)\in\mathbb{R}_{+}^{K}$ as
the vector of weighted class-conditional densities concatenated together and
$q(x,t)$ as the sampling distribution according to Dirichlet variance. This
evolution equation provides a convenient means of understanding the effect of
newly observed labels on our future acquisitions in the limit of a large
number of observations. Furthermore, it is natural to consider kernels that
are increasingly localized when analyzing the large-sample limits of kernel
methods (Devroye et al., 1996). With that in mind, it is instructive to
consider a limiting case where $\mathcal{K}(z,x)$ is given by the Dirac mass
$\delta_{z}(x)$. In that case, we can write (6.3) as simply

$\partial_{t}\alpha(x,t)=q(x,t)\eta^{K}(x),$ (23)

from which we see that the trajectories of the
$\alpha(x,t)\in\mathbb{R}_{+}^{K}$ lie along the corresponding rays
$R_{x}=\\{a\eta^{K}(x):a\geq 0\\}$, whose speed is determined by the sampling
$q(x,t)$.

We can then write

$\alpha(x,t)=\left(\int_{0}^{t}q(x,s)\,ds\right)\eta^{K}(x),$ (24)

from which we can define the amount of sampling that has occurred at $x$ up to
time $t\geq 0$ as

$\displaystyle\beta_{x}(t):=\sum_{k=1}^{K}\alpha_{k}(x,t)=\int_{0}^{t}q(x,s)\,ds,$

since the entries of $\eta^{K}(x)$ sum to $1$. By appealing to (24), we can
now write the variance at $x$ as

$\displaystyle V(\alpha(x,t))$ $\displaystyle=\frac{\sum_{i\neq
j}\alpha_{i}(x,t)\alpha_{j}(x,t)}{\beta_{x}^{2}(t)(\beta_{x}(t)+1)}=\frac{\sum_{i\neq
j}\eta_{i}(x)\eta_{j}(x)}{\beta_{x}(t)+1}=\frac{\sum_{k=1}^{K}\eta_{k}(x)(1-\eta_{k}(x))}{\beta_{x}(t)+1}$
$\displaystyle=:\frac{G(x;\eta)}{\beta_{x}(t)+1},$

where this function

$G(x;\eta)=\sum_{k=1}^{K}\eta_{k}(x)(1-\eta_{k}(x))\in\left[0,1-\frac{1}{K}\right]$
(25)

is a measure of the population level uncertainty as in in Example 5.

We now assume that the distribution of points sampled by the acquisition
function, namely
$q(x,t)=\frac{e^{\lambda(t)V(\alpha(x,t))}}{\int_{\mathcal{X}}e^{\lambda(t)V(\alpha(z,t))}\,dz}$
is asymptotically convergent as $t\to\infty$. The question of whether this
actually always occurs for solutions of the evolution equation (23) is not
simple, but as $q$ will be weakly compact in $t$ it is natural to guess that
it will approach some steady state $\bar{q}(x)$. If one were to use this time-
independent distribution throughout the entire stochastic process, we would
obtain the simple sampling relationship

$\bar{\beta}_{x}(t)=\int_{0}^{t}\bar{q}(x)=t\bar{q}(x),$

with corresponding pseudo-label densities
$\bar{\alpha}(x,t)=t\bar{q}(x)\eta^{K}(x)$.

The assumption that $q(x,t)\to\bar{q}(x)$ as $t\to\infty$ suggests then that
$\partial_{t}\alpha(x,t)\approx\partial_{t}\bar{\alpha}(x,t)$ for large enough
$t$. From this, we can approximate
$\beta_{x}(t)\approx\bar{\beta}_{x}(t)=t\bar{q}(x)$ for $t$ large. Thus, we
expect that

$q(x,t)\propto\exp\left(\frac{\lambda(t)G(x;\eta)}{t\bar{q}(x)+1}\right)\approx\exp\left(\frac{\lambda(t)G(x;\eta)}{t\bar{q}(x)}\right)$
(26)

as $t\to\infty$. This gives a non-linear algebraic relation that allows us to
narrow our search for possible limiting sampling distributions. At this point,
we consider two cases for the scaling of $\lambda(t)$: when (i)
$\lambda(t)=\lambda_{0}t^{p}=o(t)$ for $p<1$ and (ii)
$\lambda(t)=\lambda_{0}t$ for $\lambda_{0}>0$.

Case 1: When $\lambda(t)=\lambda_{0}t^{p}$ for $p<1$, then the the exponent in
the last expression of (26) decreases like $t^{-(1-p)}$ for every
$x\in\mathcal{X}$ and the result is that the limiting distribution is uniform
over the domain,
$q(x,t)\to\bar{q}(x)=\left(\int_{\mathcal{X}}\,dz\right)^{-1}$. Roughly, we
have that the equivalent expression yields

$q(x,t)=\left(\int_{\mathcal{X}}\exp\left\\{\frac{\lambda_{0}}{t^{1-p}}\left(\frac{G(z;\eta)}{\bar{q}(z)}-\frac{G(x;\eta)}{\bar{q}(x)}\right)\right\\}\,dz\right)^{-1}\to\left(\int_{\mathcal{X}}\,dz\right)^{-1}$

as $t\to\infty$ implying that the limiting distribution $\bar{q}(x)$ is
uniform over $\mathcal{X}$ when $\lambda(t)=o(t)$. We can interpret this as
asymptotic exploration that samples throughout the domain irrespective of the
corresponding marginal density $\rho(x)$ in contrast to passive sampling via
said marginal distribution.

Case 2: Now consider the case when $\lambda(t)=\lambda_{0}t$. The expression
(26) the becomes autonomous in the large $t$ limit, and suggests that the
limiting distribution in this case satisfies

$\bar{q}(x)\propto\exp\left(\frac{\lambda_{0}G(x;\eta)}{\bar{q}(x)}\right).$
(27)

In the next section, we provide numerical simulations that suggest that for
moderate values of $\lambda_{0}$ this $\bar{q}(x)$ emphasizes regions of
greater population-level uncertainty; namely, where multiple conditional
probabilities ($\eta_{i}(x)=\mathbb{P}(Y=i|X=x)$) are simultaneously large. In
other words, this limiting distribution asymptotically exploits in regions
around the true decision boundaries. For $\lambda_{0}\ll 1$, however, we
expect (27) to imply that $\bar{q}(x)$ behaves like the uniform distribution
on $\mathcal{X}$.

We notice that the two cases here both give asymptotic sampling densities
which are independent of $\rho$. In large sample regimes, we posit that
independence on $\rho$ ought to be viewed positively, meaning that we will
explore (Case 1), or exploit (Case 2), in a fashion that is unbiased by the
density $\rho$.

Of course, several steps in this derivation are formal. The passage to
“continuous time” is likely justifiable using stochastic approximation
techniques under appropriate scalings, but the details would be significant.
The use of non-Dirac kernels certainly will increase the number of possible
steady-state sampling densities, but given a density with compact support it
also seems plausible that one could show that sampling density steady states
need to be “close” to the one given by the Dirac mass kernel. Finally, the
existence of steady-state densities is not completely obvious, especially in
light of the $\alpha$ large approximations that we made. On the other hand, we
do anticipate that it should be possible to demonstrate that sampling
densities associated with the evolution equation (23) do converge to the type
of steady states we have formally described here.

Although the work here is informal, we think that it gives valuable insight
into the behavior observed in Figure 6, where DiAL effectively transitions
from exploratory to exploitative behavior. We leave more rigorous proof of the
types of formulas given in this section for future work. We now present some
numerics to illustrate the behavior of (26) in the two cases of $\lambda(t)$
that we’ve considered here.

####  6.3.1 Numerical Demonstration

We return to the mixture of Gaussians setup of Example 5 and numerically solve
the evolution equation (23) for $\alpha(x,t)$ up to $t=10^{10}$ in order to
simulate convergence of $q(x,t)\to\bar{q}(x)$ for various scalings of
$\lambda(t)$. In panel (b) of Figure 16, we plot the corresponding
distributions $\bar{q}(x)$ for algebraic scalings of $\lambda(t)=t^{p}$ for
some values of $p\in[0,1]$; note that the case of $p=0$ corresponds to
$\lambda(t)=\lambda_{0}=1$. In panel (c), we plot $\bar{q}(x)$ for
$\lambda(t)=\lambda_{0}t$ for a range of values of $\lambda_{0}$.

As suggested by our analysis above, we see in panel (b) that for $\lambda(t)$
constant (and for sufficiently small powers $p<1$), the limiting distribution
$\bar{q}(x)$ corresponds to uniform sampling over the whole domain. We can
interpret this as asymptotic exploration that samples throughout the domain
irrespective of the corresponding marginal density $\rho(x)$ in contrast to
passive sampling via said marginal distribution. As $p\to 1$, however, we
observe that $\bar{q}(x)$ focuses on the regions of the domain that lie near
the true decision boundary between the classes; we can interpret this as
asymptotic exploitation. Similarly, in panel (c) of Figure 16, we observe the
influence of $\lambda_{0}$ on $\bar{q}(x)$, including the expected behavior as
$\lambda_{0}\ll 1$ corresponding to the uniform distribution.

These numerical experiments simply demonstrate the intimately coupled nature
of the scaling of $\lambda(t)$ and the limiting distribution
$q(x,t)\to\bar{q}(x)$ of this stochastic process that represents Dirichlet
Active Learning (DiAL). In addition to further work on analyzing the steady
states of the stochastic process associated with DiAL, we suggest that an in-
depth numerical exploration into this process is warranted pursuant our
findings here.

(a) Setup of 2 classes

(b) $\bar{q}(x)$ when $\lambda(t)=t^{p}$

(c) $\bar{q}(x)$ when $\lambda(t)=\lambda_{0}t$

Figure 16: Numerical demonstration of the effect of $\lambda(t)$ on the steady
state distribution $\bar{q}(x)$ in a 1D setting. With class-conditional
distributions $p(x|y=0),p(x|y=1)$ show respectively in gray and purple panel
(a), we numerically solve (23) up to $t=10^{10}$ for different scalings of
$\lambda(t)$ to simulate convergence of $q(x,t)\to\bar{q}(x)$. Panels (b) and
(c) suggest that the limiting behavior when $\lambda(t)=\lambda_{0}t$ for
$\lambda_{0}\gg 0$ corresponds to asymptotic exploitation, whereas when
$\lambda(t)=\lambda_{0}t^{p}$ for $p<1$ or $\lambda_{0}\ll 1$ corresponds to
asymptotic exploration.

###  6.4 Consistency of Dirichlet Learning classifier

We now turn to establishing the asymptotic consistency of the underlying
classifier (3) associated with Dirichlet Learning in the binary case ($K=2$)

$\hat{y}(x)=\operatorname*{arg\\!max}_{k=1,2}\alpha_{k}(x)=\operatorname*{arg\\!max}_{k=1,2}\sum_{x^{\ell}\in\mathcal{L}_{k}}\mathcal{K}(x^{\ell},x).$

We consider an open, connected, smooth, and bounded domain
$\mathcal{X}\subset\mathbb{R}^{d}$ with $n$ observed labeled data pairs
$(X_{i},Y_{i})\in\mathbb{R}^{d}\times\\{0,1\\},i=1,2,\ldots,n$ drawn
independently and identically according to a joint distribution with density
$\nu$ with marginal over the inputs given by the density $\rho(x)$. Notice
that we have shifted the labels to be $Y_{i}\in\\{0,1\\}$ as is common
practice in analyzing methods in the binary classification case. Consider the
density-dependent heat kernel propagation
$\mathcal{K}(z,x)\equiv\mathcal{K}_{t}(z,x)$ with source $z\in\mathcal{X}$
that solves

$\begin{dcases}\begin{aligned}
\partial_{t}\mathcal{K}_{t}(z,x)-\Delta_{\rho}\mathcal{K}_{t}(z,x)&=0&x\in\mathcal{X},t>0\\\
\mathcal{K}_{0}(z,x)&=\delta_{z}(x)&\\\ \end{aligned}\end{dcases}$ (28)

where
$\Delta_{\rho}=\frac{1}{\rho}\operatorname{div}\left(\rho^{2}\nabla\cdot\right)$
is a self-adjoint diffusion operator with respect to the $\rho$-weighted inner
product. As mentioned previously, please note that the differential operators
in (28) are with respect to the variable $x$, while the variable $z$ is the
source. Furthermore, recall that depending on the domain $\mathcal{X}$, one
must assume boundary conditions to make (28) well-defined, but we will instead
state assumptions on the heat kernel, $\mathcal{K}_{t}(z,x)$, and provide
reasonable example situations in which these assumptions hold.

###### Remark 13

Asymptotic consistency of classifiers is nearly always studied under the
setting of labeled data drawn i.i.d. from the underlying data-generating
distribution. However, in the case of active learning it is important to note
that the labeled data is not drawn i.i.d. from this distribution; rather, the
sequence of observed data points is highly dependent on the previously labeled
points. Furthermore, when the active learning policy selects query points
randomly according to the distribution of acquisition function values (e.g.,
Prop. Sampling from Table 1), this evolving distribution does not align with
the underlying marginal density $\rho$. As such, we emphasize that the result
of this section is to establish a classical type of consistency property of
the Dirichlet Learning classification rule in this continuum limit setting
(i.e., infinite unlabeled data) when the labeled data is drawn i.i.d. from the
joint distribution.

The establishment of consistency of the Dirichlet Learning classifier with the
biased sampling of query points according to Dirichlet Variance is left as a
future direction of work. We posit, however, that a reasonable simplification
can come from assuming that the distribution of proportional sampling has
converged to the steady-state “asymptotic exploitation” distribution that we
derived in the previous section.

For a data density $\rho(x)$ that corresponds to the marginal distribution of
inputs for the joint distribution of $(X,Y)\in\mathcal{X}\times\\{0,1\\}$,
define the true trend function $\eta(x):=\mathbb{E}[Y|X=x]$. Showing
consistency of the Dirichlet Learning classifier can be related to the
convergence of an estimated trend function related to (30). We define the
Dirichlet Learning estimated trend function to be

$\hat{\eta}(x):=\frac{1}{n\mathbb{E}_{\rho}[\mathcal{K}_{t}(Z,x)]}\sum_{i=1}^{n}Y_{i}\mathcal{K}_{t}(X_{i},x),$
(29)

which gives that $\hat{y}(x)=\mathbbm{1}\\{\hat{\eta}(x)\geq\frac{1}{2}\\}$.

We make the following assumptions that will simplify our consistency proof:

###### Assumption 1

The true trend function, $\eta(x)$, is uniformly continuous on $\mathcal{X}$.

###### Assumption 2

Given open and connected domain $\mathcal{X}\subset\mathbb{R}^{d}$ and density
$\rho:\mathcal{X}\rightarrow\mathbb{R}_{+}$, then there exists $p\geq 2$ such
that

$\|\mathcal{K}_{t}(\cdot,x)\|_{L^{p}(\rho)}\leq t^{-\frac{d}{2}}$

for each $x\in\mathcal{X}$ and $t\leq 1$.

###### Assumption 3

Given open and connected domain $\mathcal{X}\subset\mathbb{R}^{d}$ and density
$\rho:\mathcal{X}\rightarrow\mathbb{R}_{+}$, then there exists a universal
constant $\theta>0$ such that

$\sup_{z\in\mathcal{X}}\left(\int_{\mathcal{X}}\mathcal{K}_{t}(z,x)dx\right)\leq\theta<\infty$

for all $t\leq 1$.

###### Remark 14

The question of uniform heat kernel estimates has received significant
attention in the mathematical community, and classical references on the topic
include (Davies, 1989; Grigoryan, 2009; Ouhabaz, 2009). These types of
estimates are intimately linked with log-Sobolev and isoperimetric
inequalities, which are also often linked with quantifiable sampling
estimates. One notable class where these types of bounds hold are for log-
concave distributions. Another is for densities supported on a bounded, smooth
domain, which are bounded uniformly away from zero on the entire domain. We
have assumed heat kernel estimates of this form for convenience, and certainly
this assumption could be relaxed somewhat at the cost of added complexity in
the proofs.

We now state the following theorem that implies consistency of the Dirichlet
Learning classifier decision rule (30).

###### Theorem 1

Consider the binary case of the Dirichlet Learning classifier

$\hat{y}(x):=\begin{cases}1&\text{if
}\sum_{i=1}^{n}Y_{i}\mathcal{K}_{t}(X_{i},x)\geq\sum_{i=1}^{n}(1-Y_{i})\mathcal{K}_{t}(X_{i},x)\\\
0&\text{otherwise}\\\ \end{cases},$ (30)

with $\mathcal{K}_{t}(z,x)$ the density-dependent heat kernel defined in (28)
with open, bounded, and connected domain $\mathcal{X}\subset\mathbb{R}^{d}$
and for times $t\in[0,1]$. Let $\\{(X_{i},Y_{i})\\}_{i=1}^{n}$ be samples from
the underlying joint distribution $P(X,Y)$. Assume that the true trend
function $\eta(x)=\mathbb{E}[Y|X=x]$ is uniformly continuous on $\mathcal{X}$
(Assumption 1) and that Assumptions 2 and 3 are satisfied. If $t\rightarrow 0$
and $nt^{\frac{d}{2}}\rightarrow\infty$ as $n\rightarrow\infty$, then for
$\delta\in(0,1)$, the estimated trend function (29) for the Dirichlet Learning
classifier (30) satisfies

$\int_{\mathcal{X}}\left|\hat{\eta}(x)-\eta(x)\right|\rho(x)dx\leq\Xi(\theta)\sqrt{\frac{\ln\left(\frac{1}{\delta}\right)}{n}}$

with probability greater than $1-\delta$, where the constant $\Xi(\theta)$
only depends on the universal constant $\theta<\infty$ as given in Assumption
3.

###### Remark 15

By Devroye et al ((Devroye et al., 1996), Thm. 2.3), the conclusion of Theorem
1 implies that the Dirichlet Learning classifier (30) satisfies universal
consistency, in the sense that

$\mathbb{P}\left(\mathcal{E}(\hat{y})-\mathcal{E}^{\ast}\right)\leq
2e^{-n\epsilon^{2}/\Xi(\theta)^{2}}$

for $\epsilon>0$ where $\mathcal{E}(f)=\mathbb{E}_{\nu|\nu^{n}}[f(X)\neq
Y|X_{1},Y_{1},\ldots,X_{n},Y_{n}]$ is the error probability of a
classification rule, $f$ and $\mathcal{E}^{\ast}$ is the Bayes optimal error
achieved by $y^{\ast}(x)=\mathbbm{1}\\{\eta(x)\geq\frac{1}{2}\\}$.

####  6.4.1 Proof of Theorem 1

We now present the details of the proof of Theorem 1.

Proof  First, note that the expected value of the estimated trend function
(29) satisfies:

$\displaystyle\mathbb{E}_{\nu^{n}}[\hat{\eta}(x)]$
$\displaystyle=\frac{1}{n\mathbb{E}_{\rho}[\mathcal{K}_{t}(Z,x)]}\sum_{i=1}^{n}\mathbb{E}_{\nu}[Y_{i}\mathcal{K}_{t}(X_{i},x)]$
$\displaystyle=\frac{1}{n\mathbb{E}_{\rho}[\mathcal{K}_{t}(Z,x)]}\sum_{i=1}^{n}\int_{\mathcal{X}}\eta(x_{i})\mathcal{K}_{t}(x_{i},x)\rho(x_{i})dx_{i}$
$\displaystyle=\frac{\mathbb{E}_{\rho}[\eta(Z)\mathcal{K}_{t}(Z,x)]}{\mathbb{E}_{\rho}[\mathcal{K}_{t}(Z,x)]}.$

This explicitly highlights the “averaging” nature of this heat kernel
classifier. Further, when $t\rightarrow 0$, we recover
$\mathbb{E}_{\nu^{n}}[\hat{\eta}(x)]\rightarrow\eta(x)$ due to the initial
condition that $\mathcal{K}_{0}(z,x)=\delta_{z}(x)$. This also suggests then
that as we observe more data ($n\rightarrow\infty$), it is natural to consider
$t\rightarrow 0$.

Consider the decomposition

$|\eta(x)-\hat{\eta}(x)|=\mathbb{E}_{\nu^{n}}[|\eta(x)-\hat{\eta}(x)|]+\left(|\eta(x)-\hat{\eta}(x)|-\mathbb{E}_{\nu^{n}}[|\eta(x)-\hat{\eta}(x)|]\right),$

where the expectation is taken over the randomness in the labeled data
$(X_{i},Y_{i})\sim^{iid}\nu$, and we aim to bound the first term by bounding
the terms

$\mathbb{E}_{\nu^{n}}[|\eta(x)-\hat{\eta}(x)|]\leq\underbrace{|\eta(x)-\mathbb{E}_{\nu^{n}}[\hat{\eta}(x)]|}_{A}+\underbrace{\mathbb{E}_{\nu^{n}}|\mathbb{E}_{\nu^{n}}[\hat{\eta}(x)]-\hat{\eta}(x)|}_{B}.$
(31)

##### Part 1:

Given $\tilde{\epsilon}>0$, then by the uniform continuity of $\eta$ there
exists $\xi>0$ such that $|\eta(x)-\eta(z)|\leq\tilde{\epsilon}$ for all $z\in
B(x,\xi)$. Defining
$I(t;x,\xi)=\int_{\mathcal{B}(x,\xi)^{c}}\mathcal{K}_{t}(x,z)\rho(z)dz<\mathbb{E}_{\rho}[\mathcal{K}_{t}(Z,x)]$,
then note that $I(t;x,\xi)$ is continuous with respect to $t\geq 0$ and

$\frac{I(t;x,\xi)}{\mathbb{E}_{\rho}[\mathcal{K}_{t}(Z,x)]}\rightarrow
0,\text{ as }t\rightarrow 0^{+}.$

Thus, with $\xi>0$, we can choose $t^{\ast}(\xi)>0$ such that

$\frac{I(t;x,\xi)}{\mathbb{E}_{\rho}[\mathcal{K}_{t}(Z,x)]}\leq\tilde{\epsilon}.$
(32)

Using the fact that $\eta(x)\in[0,1]$ for all $x$, then for $t\leq
t^{\ast}(\xi)$, we can rewrite the integral of $\mathbb{E}_{\rho}[A]$ as

$\displaystyle\int_{\mathcal{X}}|\eta(x)-\mathbb{E}_{\nu^{n}}[\hat{\eta}(x)]|\rho(x)dx$
$\displaystyle\leq\int_{\mathcal{X}}\int_{\mathcal{X}}|\eta(x)-\eta(z)|\frac{k_{t}(z,x)}{\mathbb{E}_{\rho}[\mathcal{K}_{t}(Z,x)]}\rho(z)dz\rho(x)dx$
$\displaystyle\leq\tilde{\epsilon}\int_{\mathcal{X}}\int_{B(x,\xi)}\mathcal{K}_{t}(z,x)\rho(z)dzdx$
$\displaystyle\qquad+\int_{\mathcal{X}}\overbrace{\int_{B(x,\xi)^{c}}\mathcal{K}_{t}(z,x)\rho(z)dz}^{I(t;x,\xi)}\frac{\rho(x)}{\mathbb{E}_{\rho}[\mathcal{K}_{t}(Z,x)]}dx$
$\displaystyle\leq\tilde{\epsilon}\int_{\mathcal{X}}\rho(x)dx+\int_{\mathcal{X}}\frac{I(t;x,\xi)}{\mathbb{E}_{\rho}[\mathcal{K}_{t}(Z,x)]}\rho(x)dx$
$\displaystyle\leq 2\tilde{\epsilon},$

where we have used (32).

##### Part 2:

For the integral of the term $B$ from (31), we can bound

$\displaystyle\mathbb{E}_{\nu^{n}}[|\mathbb{E}_{\nu^{n}}[\hat{\eta}(x)]-\hat{\eta}(x)|]$
$\displaystyle\leq\sqrt{\mathbb{E}_{\nu^{n}}[|\mathbb{E}_{\nu^{n}}[\hat{\eta}(x)]-\hat{\eta}(x)|^{2}]}$
$\displaystyle=\sqrt{\frac{\mathbb{E}_{\nu^{n}}\left\\{\left(\sum_{i=1}^{n}\mathbb{E}_{\nu}[Y_{i}\mathcal{K}_{t}(X_{i},x)]-\sum_{j=1}^{n}Y_{j}\mathcal{K}_{t}(X_{j},x)\right)^{2}\right\\}}{n^{2}(\mathbb{E}_{\rho}[\mathcal{K}_{t}(Z,x)])^{2}}}$
$\displaystyle=\sqrt{\frac{Var_{\nu}\left(Y\mathcal{K}_{t}(X,x)\right)}{n(\mathbb{E}_{\rho}[\mathcal{K}_{t}(Z,x)])^{2}}}$
$\displaystyle=\frac{1}{\sqrt{n}}\sqrt{\frac{\mathbb{E}_{\nu}\left[Y^{2}k^{2}_{t}(X,x)-2Yk_{t}(X,x)\mathbb{E}_{\rho}[\eta(Z)k_{t}(Z,x)]\right]+(\mathbb{E}_{\rho}[\eta(Z)k_{t}(Z,x)])^{2}}{(\mathbb{E}_{\rho}[\mathcal{K}_{t}(Z,x)])^{2}}}$
$\displaystyle=\frac{1}{\sqrt{n}}\sqrt{\frac{\mathbb{E}_{\rho}[\eta(X)k^{2}_{t}(X,x)]-(\mathbb{E}_{\rho}[\eta(Z)\mathcal{K}_{t}(Z,x)])^{2}}{(\mathbb{E}_{\rho}[\mathcal{K}_{t}(Z,x)])^{2}}}\qquad\qquad(Y^{2}=Y)$
$\displaystyle\leq\frac{1}{\sqrt{n}}\sqrt{\frac{\mathbb{E}_{\rho}[\eta(X)k^{2}_{t}(X,x)]}{(\mathbb{E}_{\rho}[\mathcal{K}_{t}(Z,x)])^{2}}}$
$\displaystyle\leq\frac{1}{\sqrt{n}}\sqrt{\frac{\mathbb{E}_{\rho}[k^{2}_{t}(Z,x)]}{(\mathbb{E}_{\rho}[\mathcal{K}_{t}(Z,x)])^{2}}}$
(33)

where we have used the fact that the data was sampled i.i.d. according to
$\nu$. Assumption 2 with $\|\mathcal{K}_{t}(z,\cdot)\|_{L^{p}(\rho)}$ for some
$p\geq 2$, then by Hölder’s inequality:

$\displaystyle\int_{\mathcal{X}}\mathcal{K}_{t}^{2}(z,x)\rho(z)dz$
$\displaystyle=\int_{\mathcal{X}}\mathcal{K}_{t}^{2}(z,x)\rho^{\frac{2}{p}}(z)\rho^{\frac{p-2}{p}}(z)dz$
$\displaystyle\leq\left(\int_{\mathcal{X}}\mathcal{K}_{t}^{p}(z,x)\rho(z)dz\right)^{\frac{2}{p}}\left(\int_{\mathcal{X}}\rho^{\frac{p-2}{p}\frac{p}{p-2}}(z)\right)^{\frac{p-2}{p}}$
$\displaystyle=\|\mathcal{K}_{t}(\cdot,x)\|_{L^{p}(\rho)}^{2}\left(\int_{\mathcal{X}}\rho(z)dz\right)^{\frac{p-2}{p}}$
$\displaystyle=\|\mathcal{K}_{t}(\cdot,x)\|_{L^{p}(\rho)}^{2}.$

Thus, we can simplify (6.4.1) to

$\mathbb{E}_{\nu^{n}}[|\mathbb{E}_{\nu^{n}}[\hat{\eta}(x)]-\hat{\eta}(x)|]\leq\frac{1}{\sqrt{n}}\frac{\|\mathcal{K}_{t}(\cdot,x)\|_{L^{p}(\rho)}}{\mathbb{E}_{\rho}[\mathcal{K}_{t}(Z,x)]}.$

Noting that as $t\rightarrow 0^{+}$

$\mathbb{E}_{\rho}[\mathcal{K}_{t}(Z,x)]\rightarrow\rho(x)\text{ as
}t\rightarrow 0^{+},$

then we have that there exists $\hat{t}>0$ such that for all $t\leq\hat{t}$

$\frac{\mathbb{E}_{\rho}[\mathcal{K}_{t}(Z,x)]}{\rho(x)}\geq\frac{1}{2}$ (34)

for a sufficiently small $\gamma>0$. The integral of $B$ then becomes

$\displaystyle\int_{\mathcal{X}}\mathbb{E}_{\nu^{n}}[|\mathbb{E}_{\nu^{n}}[\hat{\eta}(x)]-\hat{\eta}(x)|]\rho(x)dx$
$\displaystyle\leq\frac{1}{\sqrt{n}}\int_{\mathcal{X}}\|\mathcal{K}_{t}(z,\cdot)\|_{L^{p}(\rho)}\frac{\rho(x)}{\mathbb{E}_{\rho}[\mathcal{K}_{t}(Z,x)]}dx$
$\displaystyle\leq\frac{2}{\sqrt{nt^{d}}}Vol(\mathcal{X})$
$\displaystyle\rightarrow 0\quad\text{as }n\rightarrow\infty.$ (35)

Thus, as long as $nt^{d}\rightarrow\infty$ as $n\rightarrow\infty$ and
$t\rightarrow 0$, we can take $n$ sufficiently large to bound (31) by
combining the bounds on $\mathbb{E}_{\rho}[A]$ and $\mathbb{E}_{\rho}[B]$ with
(6.4.1) to get

$\mathbb{E}_{\nu^{n}}\left[\int_{\mathcal{X}}|\eta(x)-\hat{\eta}(x)|\rho(x)dx\right]\leq
3\tilde{\epsilon}=\frac{\epsilon}{2},$

where we take $\tilde{\epsilon}=\frac{\epsilon}{6}.$

##### Part 3:

We can then use McDiarmid’s inequality to give a high probability bound on the
difference

$\int_{\mathcal{X}}\left(|\eta(x)-\hat{\eta}(x)|-\mathbb{E}_{\nu^{n}}[|\eta(x)-\hat{\eta}(x)|]\right)dx.$

Let $\mathcal{D}=\\{(x_{i},y_{i})\\}_{i=1}^{n}$ be the fixed training
(labeled) data, and let $\mathcal{D}^{(i)}$ be a copy of $\mathcal{D}$ with
the $i^{th}$ pair replaced with $(\tilde{x}_{i},\tilde{y}_{i})$. Denote by
$\hat{\eta}_{i}(x)$ the corresponding estimated trend function for
$\mathcal{D}^{(i)}$. Then, we can straightforwardly bound the following
difference by using Assumption 3

$\displaystyle\int_{\mathcal{X}}|\eta(x)-\hat{\eta}(x)|\rho(x)dx$
$\displaystyle-\int_{\mathcal{X}}|\eta(x)-\hat{\eta}_{i}(x)|\rho(x)dx\leq\int_{\mathcal{X}}|\hat{\eta}(x)-\hat{\eta}_{i}(x)|\rho(x)dx$
$\displaystyle=\frac{1}{n}\int_{\mathcal{X}}|Y_{i}\mathcal{K}_{t}(X_{i},x)-\tilde{Y_{i}}\mathcal{K}_{t}(\tilde{X}_{i},x)|\frac{\rho(x)}{\mathbb{E}_{\rho}[\mathcal{K}_{t}(Z,x)]}dx$
$\displaystyle\leq\frac{4}{n}\sup_{z}\left\\{\int_{\mathcal{X}}\mathcal{K}_{t}(z,x)dx\right\\}$
$\displaystyle\leq\frac{4\theta}{n}$

where the constant $\theta$ was defined in Assumption 3 and have used (34)
with $n\geq 2$ trivially. McDiarmid’s inequality yields that

$\displaystyle\mathbb{P}$
$\displaystyle\left(\int_{\mathcal{X}}|\eta(x)-\hat{\eta}(x)|\rho(x)dx>\epsilon\right)$
$\displaystyle\leq\mathbb{P}\left(\int_{\mathcal{X}}|\eta(x)-\hat{\eta}(x)|\rho(x)dx-\mathbb{E}_{\nu^{n}}\left[\int_{\mathcal{X}}|\eta(x)-\hat{\eta}(x)|\rho(x)dx\right]>\frac{\epsilon}{2}\right)$
$\displaystyle\leq e^{-\frac{n\epsilon^{2}}{32\theta^{2}}},$

our desired result that, with high probability, the $\ell^{1}$ difference
between $\eta(x)$ and $\hat{\eta}(x)$ is small as $n\rightarrow\infty$ and
$t\rightarrow 0^{+}$ with $nt^{d}\rightarrow\infty$ with $t$ sufficiently
small.

Acknowledgments and Disclosure of Funding

KM acknowledges support from the Peter J. O’Donnell Jr. Postdoctoral
Fellowship and NSF IFML grant 2019844. RM acknowledges partial support from
NSF-DMS 2307971 and the Simons Foundation MP-TSM.

## References

  * Agarwal et al. (2021) Alekh Agarwal, Nan Jiang, Sham M Kakade, and Wen Sun.  Reinforcement Learning: Theory and Algorithms, 2021.  Accessed via https://rltheorybook.github.io/rltheorybook_AJKS.pdf.
  * Ash et al. (2020) Jordan T. Ash, Chicheng Zhang, Akshay Krishnamurthy, John Langford, and Alekh Agarwal.  Deep batch active learning by diverse, uncertain gradient lower bounds.  In _8th International Conference on Learning Representations (ICLR)_ , 2020.
  * Bachem et al. (2017) Olivier Bachem, Mario Lucic, and Andreas Krause.  Practical coreset constructions for machine learning, 2017.  arXiv preprint arXiv:1703.06476 [stat.ML].
  * Balcan et al. (2007) Maria-Florina Balcan, Andrei Broder, and Tong Zhang.  Margin based active learning.  In _International Conference on Computational Learning Theory_ , volume 4539, pages 35–50. Springer Berlin Heidelberg, 2007.  ISBN 978-3-540-72925-9.  doi: 10.1007/978-3-540-72927-3˙5.  URL https://doi.org/10.1007/978-3-540-72927-3_5.
  * Balcan et al. (2009) Maria-Florina Balcan, Alina Beygelzimer, and John Langford.  Agnostic active learning.  _Journal of Computer and System Sciences_ , 75(1):78–89, 2009.  doi: 10.1016/j.jcss.2008.07.003.  URL https://doi.org/10.1016/j.jcss.2008.07.003.
  * Berthelot et al. (2019) David Berthelot, Nicholas Carlini, Ian Goodfellow, Nicolas Papernot, Avital Oliver, and Colin A Raffel.  Mixmatch: A holistic approach to semi-supervised learning.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, _Advances in Neural Information Processing Systems_ , volume 32. Curran Associates, Inc., 2019.  URL https://proceedings.neurips.cc/paper_files/paper/2019/file/1cd138d0499a68f4bb72bee04bbec2d7-Paper.pdf.
  * Bertozzi and Flenner (2016) Andrea L. Bertozzi and Arjuna Flenner.  Diffuse interface models on graphs for classification of high dimensional data.  _SIAM Review_ , 2016.  doi: 10.1137/16M1070426.
  * Bertozzi et al. (2018) Andrea L. Bertozzi, Xiyang Luo, Andrew M. Stuart, and Konstantinos C. Zygalakis.  Uncertainty quantification in graph-based classification of high dimensional data.  _SIAM/ASA Journal on Uncertainty Quantification_ , 6(2):568–595, 2018.
  * Bertozzi et al. (2021) Andrea L Bertozzi, Bamdad Hosseini, Hao Li, Kevin Miller, and Andrew M Stuart.  Posterior consistency of semi-supervised regression on graphs.  _Inverse Problems_ , 37(10):105011, September 2021\.  doi: 10.1088/1361-6420/ac1e80.  URL https://doi.org/10.1088/1361-6420/ac1e80.
  * Cai et al. (2017) Wenbin Cai, Muhan Zhang, and Ya Zhang.  Batch mode active learning for regression with expected model change.  _IEEE Transactions on Neural Networks and Learning Systems_ , 28(7):1668–1681, July 2017.  doi: 10.1109/tnnls.2016.2542184.  URL https://doi.org/10.1109/tnnls.2016.2542184.
  * Calder and Slepčev (2020) Jeff Calder and Dejan Slepčev.  Properly-weighted graph Laplacian for semi-supervised learning.  _Applied Mathematics & Optimization_, 82(3):1111–1159, December 2020.  ISSN 1432-0606.  doi: 10.1007/s00245-019-09637-3.  URL https://doi.org/10.1007/s00245-019-09637-3.
  * Calder et al. (2020) Jeff Calder, Brendan Cook, Matthew Thorpe, and Dejan Slepčev.  Poisson learning: Graph-based semi-supervised learning at very low label rates.  In _Proceedings of the 37th International Conference on Machine Learning_ , pages 1306–1316. Proceedings of Machine Learning Research, November 2020.  URL https://proceedings.mlr.press/v119/calder20a.html.  ISSN: 2640-3498.
  * Calder et al. (2023) Jeff Calder, Dejan Slepčev, and Matthew Thorpe.  Rates of convergence for Laplacian semi-supervised learning with low labeling rates.  _Research in the Mathematical Sciences_ , 10(1):10, February 2023.  ISSN 2197-9847.  doi: 10.1007/s40687-022-00371-x.  URL https://doi.org/10.1007/s40687-022-00371-x.
  * Chapman et al. (2023) James Chapman, Bohan Chen, Zheng Tan, Jeffrey Calder, Kevin Miller, and Andrea Bertozzi.  Novel batch active learning approach and its application on the synthetic aperture radar datasets.  In _Proceedings of Society of Photo-Optical Instrumentation Engineers (SPIE) 2023 Conference on Defense + Commercial Sensing_. SPIE, 2023\.
  * Chen et al. (2023) Yifan Chen, Ethan N. Epperly, Joel A. Tropp, and Robert J. Webber.  Randomly pivoted Cholesky: Practical approximation of a kernel matrix with few entry evaluations, 2023.  arXiv preprint arXiv:2207.06503 [math.NA].
  * Chung (1997) F. R. K. Chung.  _Spectral Graph Theory_.  American Mathematical Society, 1997.
  * Cloninger and Mhaskar (2021) A. Cloninger and H. N. Mhaskar.  Cautious active clustering.  _Applied and Computational Harmonic Analysis_ , 54:44–74, September 2021.  ISSN 1063-5203.  doi: 10.1016/j.acha.2021.02.002.  URL https://www.sciencedirect.com/science/article/pii/S1063520321000154.
  * Cover and Hart (1967) T. Cover and P. Hart.  Nearest neighbor pattern classification.  _IEEE Transactions on Information Theory_ , 13(1):21–27, 1967.  doi: 10.1109/TIT.1967.1053964.
  * Dasarathy et al. (2015) Gautam Dasarathy, Robert Nowak, and Xiaojin Zhu.  S2: An efficient graph based active learning algorithm with application to nonparametric classification.  In Peter Grünwald, Elad Hazan, and Satyen Kale, editors, _Proceedings of The 28th Conference on Learning Theory_ , volume 40 of _Proceedings of Machine Learning Research_ , pages 503–522, Paris, France, July 2015. Proceedings of Machine Learning Research.  URL https://proceedings.mlr.press/v40/Dasarathy15.html.
  * Dasgupta (2006) Sanjoy Dasgupta.  Coarse sample complexity bounds for active learning.  In _Advances in Neural Information Processing Systems_ , volume 18, pages 235–242. MIT Press, 2006.  URL https://papers.nips.cc/paper/2005/hash/6e82873a32b95af115de1c414a1849cb-Abstract.html.
  * Dasgupta (2011) Sanjoy Dasgupta.  Two faces of active learning.  _Theoretical Computer Science_ , 412(19):1767–1781, April 2011.  doi: 10.1016/j.tcs.2010.12.054.  URL https://doi.org/10.1016/j.tcs.2010.12.054.
  * Dasgupta and Hsu (2008) Sanjoy Dasgupta and Daniel Hsu.  Hierarchical sampling for active learning.  In _Proceedings of the 25th International Conference on Machine Learning_ , pages 208–215, Helsinki, Finland, July 2008. Association for Computing Machinery.  ISBN 978-1-60558-205-4.  doi: 10.1145/1390156.1390183.  URL https://doi.org/10.1145/1390156.1390183.
  * Davies (1989) Edward Brian Davies.  _Heat Kernels and Spectral Theory_.  Number 92. Cambridge University Press, 1989.
  * Deshpande and Vempala (2006) Amit Deshpande and Santosh Vempala.  Adaptive sampling and fast low-rank matrix approximation.  In Josep Díaz, Klaus Jansen, José D. P. Rolim, and Uri Zwick, editors, _Approximation, Randomization, and Combinatorial Optimization. Algorithms and Techniques_ , pages 292–303, Berlin, Heidelberg, 2006. Springer Berlin Heidelberg.  ISBN 978-3-540-38045-0.
  * Devroye et al. (1996) Luc Devroye, László Györfi, and Gábor Lugosi.  _A Probabilistic Theory of Pattern Recognition_ , volume 31 of _Stochastic Modelling and Applied Probability_.  Springer, 1996.  ISBN 978-1-4612-0711-5.
  * Gal et al. (2017) Yarin Gal, Riashat Islam, and Zoubin Ghahramani.  Deep Bayesian active learning with image data.  In _Proceedings of the 34th International Conference on Machine Learning_ , pages 1183–1192, Sydney, NSW, Australia, August 2017. Journal of Machine Learning Research.
  * García Trillos et al. (2021) Nicolás García Trillos, Franca Hoffmann, and Bamdad Hosseini.  Geometric structure of graph Laplacian embeddings.  _Journal of Machine Learning Research_ , 22(63):1–55, 2021.
  * Grigoryan (2009) Alexander Grigoryan.  _Heat kernel and analysis on manifolds_ , volume 47.  American Mathematical Soc., 2009.
  * Hanneke (2007) Steve Hanneke.  A bound on the label complexity of agnostic active learning.  In _Proceedings of the 24th International Conference on Machine Learning_ , pages 353–360, New York, NY, USA, June 2007. Association for Computing Machinery.  ISBN 978-1-59593-793-3.  doi: 10.1145/1273496.1273541.
  * Hanneke (2014) Steve Hanneke.  Theory of disagreement-based active learning.  _Foundations and Trends® in Machine Learning_ , 7(2-3):131–309, June 2014.  doi: 10.1561/2200000037.  URL https://doi.org/10.1561/2200000037.
  * Hanneke and Yang (2015) Steve Hanneke and Liu Yang.  Minimax analysis of active learning.  _Journal of Machine Learning Research_ , 16(109):3487–3602, 2015.  ISSN 1533-7928.  URL http://jmlr.org/papers/v16/hanneke15a.html.
  * Huang et al. (2010) Sheng-jun Huang, Rong Jin, and Zhi-Hua Zhou.  Active learning by querying informative and representative examples.  In J. D. Lafferty, C. K. I. Williams, J. Shawe-Taylor, R. S. Zemel, and A. Culotta, editors, _Advances in Neural Information Processing Systems 23_ , pages 892–900. Curran Associates, Inc., 2010.  URL http://papers.nips.cc/paper/4176-active-learning-by-querying-informative-and-representative-examples.pdf.
  * Ji and Han (2012) Ming Ji and Jiawei Han.  A variance minimization criterion to active learning on graphs.  In _Artificial Intelligence and Statistics_ , pages 556–564, March 2012.  URL http://proceedings.mlr.press/v22/ji12.html.
  * Jiang and Gupta (2019) Heinrich Jiang and Maya Gupta.  Minimum-margin active learning.  May 2019.  URL http://arxiv.org/abs/1906.00025.  arXiv preprint arXiv:1906.00025 [cs, stat].
  * Jun and Nowak (2018) Kwang-Sung Jun and Robert Nowak.  Chapter 10 - Bayesian active learning on graphs.  In Petar M. Djurić and Cédric Richard, editors, _Cooperative and Graph Signal Processing_ , pages 283–297. Academic Press, January 2018\.  ISBN 978-0-12-813677-5.  doi: 10.1016/B978-0-12-813677-5.00010-9.  URL https://www.sciencedirect.com/science/article/pii/B9780128136775000109.
  * Kapoor et al. (2007) Ashish Kapoor, Kristen Grauman, Raquel Urtasun, and Trevor Darrell.  Active learning with Gaussian processes for object categorization.  In _2007 IEEE 11th international conference on computer vision_ , pages 1–8. IEEE, 2007.
  * Karzand and Nowak (2020) Mina Karzand and Robert D. Nowak.  Maximin active learning in overparameterized model classes.  _IEEE Journal on Selected Areas in Information Theory_ , 1(1):167–177, May 2020.  ISSN 2641-8770.  doi: 10.1109/JSAIT.2020.2991518.
  * Kirsch et al. (2022) Andreas Kirsch, Sebastian Farquhar, Parmida Atighehchian, Andrew Jesson, Frederic Branchaud-Charron, and Yarin Gal.  Stochastic batch acquisition for deep active learning, 2022.  arXiv preprint arXiv:2106.12059 [cs.LG].
  * Krause and Guestrin (2007) Andreas Krause and Carlos Guestrin.  Nonmyopic active learning of gaussian processes: an exploration-exploitation approach.  In _Proceedings of the 24th international conference on Machine learning_ , pages 449–456, 2007.
  * Kushnir and Venturi (2020) Dan Kushnir and Luca Venturi.  Diffusion-based deep active learning, March 2020.  arXiv preprint arXiv:2003.10339.
  * LeCun and Cortes (2010) Yann LeCun and Corinna Cortes.  MNIST handwritten digit database.  2010\.  URL http://yann.lecun.com/exdb/mnist/.
  * Ma et al. (2013) Yifei Ma, Roman Garnett, and Jeff Schneider.  $\Sigma$-optimality for active learning on Gaussian random fields.  In C. J. C. Burges, L. Bottou, M. Welling, Z. Ghahramani, and K. Q. Weinberger, editors, _Advances in Neural Information Processing Systems 26_ , pages 2751–2759. Curran Associates, Inc., 2013.  URL http://papers.nips.cc/paper/4951–optimality-for-active-learning-on-gaussian-random-fields.pdf.
  * Miller and Bertozzi (2023) Kevin Miller and Andrea L. Bertozzi.  Model-change active learning in graph-based semi-supervised learning.  October 2023.  URL http://arxiv.org/abs/2110.07739.  To appear in Springer Nature Communications on Applied Mathematics and Computation (CAMC).
  * Miller and Calder (2023) Kevin Miller and Jeff Calder.  Poisson reweighted Laplacian uncertainty sampling for graph-based active learning, 2023.  To appear in SIAM Journal on Mathematics of Data Science. arXiv preprint arXiv:2210.15786.
  * Miller et al. (2020) Kevin Miller, Hao Li, and Andrea L. Bertozzi.  Efficient graph-based active learning with probit likelihood via Gaussian approximations.  In _Workshop on Experimental Design and Active Learning_ , July 2020\.  URL http://arxiv.org/abs/2007.11126.  arXiv: 2007.11126.
  * Miller et al. (2022) Kevin Miller, John Mauro, Jason Setiadi, Xoaquin Baca, Zhan Shi, Jeffrey Calder, and Andrea Bertozzi.  Graph-based active learning for semi-supervised classification of sar data.  In _Proceedings of Society of Photo-Optical Instrumentation Engineers (SPIE) 2022 Conference on Defense + Commercial Sensing_. SPIE, 2022\.
  * (47) James Murphy.  Learning by active nonlinear diffusion (LAND) codebase.  URL https://jmurphy.math.tufts.edu/Code/.  https://jmurphy.math.tufts.edu/Code/.
  * Murphy and Maggioni (2019) James M. Murphy and Mauro Maggioni.  Unsupervised clustering and active learning of hyperspectral images with nonlinear diffusion.  _IEEE Transactions on Geoscience and Remote Sensing_ , 57(3):1829–1845, March 2019.  ISSN 1558-0644.  doi: 10.1109/TGRS.2018.2869723.
  * Musco and Woodruff (2017) Cameron Musco and David P. Woodruff.  Sublinear time low-rank approximation of positive semidefinite matrices.  In _2017 IEEE 58th Annual Symposium on Foundations of Computer Science (FOCS)_ , pages 672–683, 2017.  doi: 10.1109/FOCS.2017.68.
  * Nadaraya (1964) E. A. Nadaraya.  On estimating regression.  _Theory of Probability & Its Applications_, 9(1):141–142, 1964.  doi: 10.1137/1109020.  URL https://doi.org/10.1137/1109020.
  * Ouhabaz (2009) El-Maati Ouhabaz.  _Analysis of Heat Equations on Domains.(LMS-31)_.  Princeton University Press, 2009.
  * Parzen (1962) Emanuel Parzen.  On estimation of a probability density function and mode.  _The Annals of Mathematical Statistics_ , 33(3):1065–1076, 1962.  ISSN 00034851.  URL http://www.jstor.org/stable/2237880.
  * Qiao et al. (2019) Yi-Ling Qiao, Chang Xin Shi, Chenjian Wang, Hao Li, Matt Haberland, Xiyang Luo, Andrew M. Stuart, and Andrea L. Bertozzi.  Uncertainty quantification for semi-supervised multi-class classification in image processing and ego-motion analysis of body-worn videos.  _Image Processing: Algorithms and Systems_ , 2019.  doi: 10.2352/issn.2470-1173.2019.11.ipas-264.
  * Rasmussen and Williams (2006) Carl Edward Rasmussen and Christopher K. I. Williams.  _Gaussian Processes for Machine Learning_.  Adaptive computation and machine learning. MIT Press, Cambridge, Mass, 2006.  ISBN 978-0-262-18253-9.
  * Ren et al. (2021) Pengzhen Ren, Yun Xiao, Xiaojun Chang, Po-Yao Huang, Zhihui Li, Brij B. Gupta, Xiaojiang Chen, and Xin Wang.  A survey of deep active learning.  _ACM Comput. Surv._ , 54(9), oct 2021.  ISSN 0360-0300.  doi: 10.1145/3472291.  URL https://doi.org/10.1145/3472291.
  * Riis et al. (2022) Christoffer Riis, Francisco Antunes, Frederik Boe Hüttel, Carlos Lima Azevedo, and Francisco C. Pereira.  Bayesian active learning with fully Bayesian Gaussian processes.  In Alice H. Oh, Alekh Agarwal, Danielle Belgrave, and Kyunghyun Cho, editors, _Advances in Neural Information Processing Systems_ , 2022.  URL https://openreview.net/forum?id=yhZLEvmyHYQ.
  * Rittler and Chaudhuri (2023) Nicholas Rittler and Kamalika Chaudhuri.  A two-stage active learning algorithm for k-nearest neighbors.  In Andreas Krause, Emma Brunskill, Kyunghyun Cho, Barbara Engelhardt, Sivan Sabato, and Jonathan Scarlett, editors, _Proceedings of the 40th International Conference on Machine Learning_ , volume 202 of _Proceedings of Machine Learning Research_ , pages 29103–29129. PMLR, 23–29 Jul 2023.  URL https://proceedings.mlr.press/v202/rittler23a.html.
  * Rosenblatt (1956) Murray Rosenblatt.  Remarks on some nonparametric estimates of a density function.  _The Annals of Mathematical Statistics_ , 27(3):832 – 837, 1956.  doi: 10.1214/aoms/1177728190.  URL https://doi.org/10.1214/aoms/1177728190.
  * Schreiter et al. (2015) Jens Schreiter, Duy Nguyen-Tuong, Mona Eberts, Bastian Bischoff, Heiner Markert, and Marc Toussaint.  Safe exploration for active learning with Gaussian processes.  In _Machine Learning and Knowledge Discovery in Databases: European Conference, ECML PKDD 2015, Porto, Portugal, September 7-11, 2015, Proceedings, Part III 15_ , pages 133–149. Springer, 2015.
  * Sener and Savarese (2018) Ozan Sener and Silvio Savarese.  Active learning for convolutional neural networks: A core-set approach, June 2018.  URL http://arxiv.org/abs/1708.00489.  arXiv preprint arXiv: 1708.00489.
  * Settles (2012) Burr Settles.  _Active Learning_ , volume 6.  Morgan & Claypool Publishers LLC, June 2012.  doi: 10.2200/s00429ed1v01y201207aim018.  URL https://doi.org/10.2200/s00429ed1v01y201207aim018.
  * Sohn et al. (2020) Kihyuk Sohn, David Berthelot, Nicholas Carlini, Zizhao Zhang, Han Zhang, Colin A Raffel, Ekin Dogus Cubuk, Alexey Kurakin, and Chun-Liang Li.  Fixmatch: Simplifying semi-supervised learning with consistency and confidence.  In H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin, editors, _Advances in Neural Information Processing Systems_ , volume 33, pages 596–608. Curran Associates, Inc., 2020.  URL https://proceedings.neurips.cc/paper_files/paper/2020/file/06964dce9addb1c5cb5d6e3d9838f733-Paper.pdf.
  * Sutton and Barto (2018) Richard S. Sutton and Andrew G. Barto.  _Reinforcement Learning: An Introduction_.  The MIT Press, second edition, 2018.  URL http://incompleteideas.net/book/the-book-2nd.html.
  * Tong and Koller (2001) Simon Tong and Daphne Koller.  Support vector machine active learning with applications to text classification.  _Journal of Machine Learning Research_ , 2(Nov):45–66, 2001.  ISSN ISSN 1533-7928.  URL http://www.jmlr.org/papers/v2/tong01a.html.
  * von Luxburg (2007) Ulrike von Luxburg.  A tutorial on spectral clustering.  _Statistics and Computing_ , 17(4):395–416, December 2007.  ISSN 1573-1375.  doi: 10.1007/s11222-007-9033-z.  URL https://doi.org/10.1007/s11222-007-9033-z.
  * Watson (1964) Geoffrey S. Watson.  Smooth regression analysis.  _Sankhyā: The Indian Journal of Statistics, Series A (1961-2002)_ , 26(4):359–372, 1964.  ISSN 0581572X.  URL http://www.jstor.org/stable/25049340.
  * Xiao et al. (2017) Han Xiao, Kashif Rasul, and Roland Vollgraf.  Fashion-mnist: a novel image dataset for benchmarking machine learning algorithms, 2017.  URL http://arxiv.org/abs/1708.07747.
  * Yang et al. (2023) Shuo Yang, Yijun Dong, Rachel Ward, Inderjit S. Dhillon, Sujay Sanghavi, and Qi Lei.  Sample Efficiency of Data Augmentation Consistency Regularization.  In _Proceedings of The 26th International Conference on Artificial Intelligence and Statistics_ , pages 3825–3853. PMLR, April 2023\.  URL https://proceedings.mlr.press/v206/yang23c.html.  ISSN: 2640-3498.
  * Zhan et al. (2022) Xueying Zhan, Yaowei Wang, and Antoni B. Chan.  Asymptotic optimality for active learning processes.  In _The 38th Conference on Uncertainty in Artificial Intelligence_ , 2022.
  * Zhou et al. (2004) Dengyong Zhou, Olivier Bousquet, Thomas Navin Lal, Jason Weston, and Bernhard Schölkopf.  Learning with local and global consistency.  In _Advances in Neural Information Processing Systems 16_ , pages 321–328. MIT Press, 2004.
  * Zhu et al. (2003a) Xiaojin Zhu, Zoubin Ghahramani, and John Lafferty.  Semi-supervised learning using Gaussian fields and harmonic functions.  In _Proceedings of the 20th International Conference on International Conference on Machine Learning_ , pages 912–919, Washington, DC, USA, August 2003a. AAAI Press.  ISBN 978-1-57735-189-4.
  * Zhu et al. (2003b) Xiaojin Zhu, John Lafferty, and Zoubin Ghahramani.  Combining active learning and semi-supervised learning using Gaussian fields and harmonic functions.  In _International Conference on Machine Learning (ICML) 2003 workshop on The Continuum from Labeled to Unlabeled Data in Machine Learning and Data Mining_ , pages 58–65, 2003b.