Files
thesis/content/method/method.tex
T

1020 lines
50 KiB
TeX

\section{Problem Formulation}
\label{sec:problem-formulation}
% Motivation for this section — why accuracy fails under class imbalance, why
% TPR@FPR gives direct operating-point control, and the Neyman-Pearson framing
% that underpins ROLL — is developed in \Cref{chapter:related_work}
% (\Cref{sec:imbalanced-tpr-fpr,sec:neyman-pearson}). This section establishes
% only the notation and objective used throughout the chapter.
Given a dataset $\mathcal{D} = \{(\mathbf{x}_i, y_i)\}_{i=1}^{n}$ with inputs
$\mathbf{x}_i \in \mathcal{X}$ and binary labels $y_i \in \{0, 1\}$, we learn a
scoring function $f_\theta : \mathcal{X} \to \mathbb{R}$ by fitting.
The positive and negative subsets of $\mathcal{D}$ are:
\[
\mathcal{D}_1 = \{(\mathbf{x}_i, y_i) \in \mathcal{D} \mid y_i = 1\},
\qquad
\mathcal{D}_0 = \{(\mathbf{x}_i, y_i) \in \mathcal{D} \mid y_i = 0\}.
\]
Applying a threshold $\tau$ produces predictions $\hat{y} = \mathbf{1}[f_\theta(\mathbf{x}) > \tau]$.
The \emph{true positive rate} (TPR) and \emph{false positive rate} (FPR) at $\tau$ are:
\begin{align}
\operatorname{TPR}(\tau) &= P\!\bigl(f_\theta(\mathbf{x}) > \tau \mid y = 1\bigr), \\
\operatorname{FPR}(\tau) &= P\!\bigl(f_\theta(\mathbf{x}) > \tau \mid y = 0\bigr).
\end{align}
In our scenario, we either wish to learn $f_{\theta}$ in order to maximize $\operatorname{TPR}$
subject to a fixed $\operatorname{FPR}$, or vice versa. Our fixed $\operatorname{TPR}$ or $\operatorname{FPR}$
will be fixed to a value $\alpha \in (0, 1)$.
Therefor, when maximizing our $\operatorname{TPR}$ subject to a fixed $\operatorname{FPR}$, we are solving for:
\begin{equation}
\max_\theta \; \operatorname{TPR}\!\bigl(\tau^*(\theta)\bigr)
\quad \text{s.t.} \quad \operatorname{FPR}\!\bigl(\tau^*(\theta)\bigr) = \alpha,
\label{eq:objective}
\end{equation}
where $\tau^*(\theta)$ is the threshold that achieves exactly $\operatorname{FPR} = \alpha$
for the current model $f_\theta$.
Conversley, when we are minimizing $\operatorname{FPR}$ subject to a fixed $\operatorname{TPR}$ we are solving for:
\begin{equation}
\min_\theta \; \operatorname{FPR}\!\bigl(\tau^*(\theta)\bigr)
\quad \text{s.t.} \quad \operatorname{TPR}\!\bigl(\tau^*(\theta)\bigr) = \alpha,
\label{eq:objective}
\end{equation}
The ROLL framework, introduced in
\Cref{sec:roll-framework}, provides a method of differentiable surrogate for these objectives.
%------------------------------------------------
\section{The ROLL Framework}
\label{sec:roll-framework}
In order to allow gradient-based learning methods to optimize TPR or FPR directly, we introduce
ROLL - Rate Optimized Likelyhood Loss.
The core idea of ROLL is as follows: when training, the user selects that they would like to optimize
TPR given some fixed FPR. Then, during training, on every training step, the actual TPR is calculated for
the given FPR, and the fitting algorithm tries to maximize the TPR directly. This way, while training, we gain direct control over the performance at the threshold that we care about - avoiding needing to pick out a threshold later.
In Neyman-Pearson terms~\cite{tong2016neyman}, our algorithm wishes to directly minimize the Type-1 error given a fixed Type-2 error, or vice versa.
In practice, however, things are not so simple.
\subsection{The Differentiability Problem}
\label{sec:differentiability-problem}
% TPR and FPR are defined via threshold comparisons on model scores, making them
% piecewise-constant functions of the model parameters with zero gradient almost
% everywhere — standard backpropagation cannot move through them. Show this
% concretely with the step-function definition of TPR/FPR and explain why
% common surrogates (AUC proxies, hinge-based losses) only indirectly target
% the operating point a practitioner cares about.
Directly optimizing TPR for FPR, or vice versa, is an indefferentiable problem.
For this sub-section, we will reffer to an example of maximizing TPR given a fixed FPR,
though this can be reversed.
For any loss $\mathcal{L}(\theta)$, in order
to continue fitting our model, we must derive some gradient
$\frac{\partial \mathcal{L}(\theta)}{\partial f_\theta(\mathbf{x}_i)}$ where $f_\theta(\mathbf{x}_i)$ is the score
given to data point $\mathbf{x}_i$ by our model $f_\theta$.
There are two inherent steps built-in to our methods:
\begin{itemize}
\item Based on fixed rate, calculate threshold that gives said rate.
For example, when fixing FPR = $\alpha$, find some threshold $\tau$ s.t. $\mathbb{E}[1_{f_{\theta}\mathbf{x}_i > \tau | y_i = 0}] = \alpha$
\item Given the calculated threshold, calculate the rate we would like to maximize
(in our example, maximize TPR)
\end{itemize}
TPR or FPR calculation, given a threshold, is non-differentiable.
The reason for that is that the TPR and FPR as a function of $f_\theta(\mathbf{x}_i)$ are
non-continuous. If $f_\theta(\mathbf{x}_i) \neq \tau$, that is, the score is not exactly on the threshold, then changing it
by some $\varepsilon$ will not affect either TPR or FPR at all, and the derivative of the function w.r.t.\ $f_\theta(\mathbf{x}_i)$ will be zero. Likewise,
if $f_\theta(\mathbf{x}_i) = \tau$, the derivative of TPR/FPR w.r.t.\ $f_\theta(\mathbf{x}_i)$ becomes infinite. This derivative cannot be used
to continuously fit our model.
This necessitates a different approach, to make these calculations differentiable w.r.t the scores themselves, in order to fit this type of optimization into a learning platform.
\subsection{Score Distribution Fitting as a Solution}
\label{sec:score-distribution-fitting}
% The key insight: instead of thresholding scores directly, fit a differentiable
% probability distribution to the positive-class and negative-class score outputs.
% TPR and FPR then become integrals of that distribution above/below a threshold,
% which are smooth and differentiable in the distribution parameters — and through
% them, in the model weights. Explain intuitively why this preserves the semantics
% of TPR/FPR while making the objective amenable to gradient descent.
The way to overcome this problem in differentiability is not to calculate the TPR and FPR
directly, but to estimate them by fitting probability distributions to the score output of the model at
every step. Contiguous probability distributions are contiguous w.r.t their inputs. Additionally, we can calculate the TPR, FPR given the probability distributions CDF and ICDF functions. If we can gaurantee these functions to be differentiable, which proved to be the main challenge of this work, we can gaurantee this process to be differentiable.
\subsection{General ROLL Formulation and Derivation}
\label{sec:roll-formulation}
% Present the abstract ROLL loss in terms of a generic distribution $\mathcal{F}$
% fitted to model scores, without yet committing to a specific family. Derive the
% expression for TPR@FPR$_0$ as a function of $\mathcal{F}$'s parameters, and show
% that minimizing the loss is equivalent to maximizing TPR subject to the FPR
% constraint. Include the full derivation here; this is the core theoretical
% contribution of the chapter.
Using the dataset $\mathcal{D} = \{(\mathbf{x}_i, y_i)\}_{i=1}^{n}$ as introduced in
\Cref{sec:problem-formulation}, let $\mathcal{B} \subseteq \mathcal{D}$ denote a training batch
over which we compute the loss.
We then subdivide the batch into positive and negative samples:
%
\begin{align*}
\mathcal{B}_1 &= \{(\mathbf{x}_i, y_i) \in \mathcal{B} \mid y_i = 1\} \\
\mathcal{B}_0 &= \{(\mathbf{x}_i, y_i) \in \mathcal{B} \mid y_i = 0\}
\end{align*}
%
and extract the corresponding model score vectors:
%
\begin{equation*}
\mathbf{s}^{(1)} = \bigl(f_\theta(\mathbf{x}_i)\bigr)_{(\mathbf{x}_i,\, y_i)\,\in\,\mathcal{B}_1},
\qquad
\mathbf{s}^{(0)} = \bigl(f_\theta(\mathbf{x}_i)\bigr)_{(\mathbf{x}_i,\, y_i)\,\in\,\mathcal{B}_0}
\end{equation*}
We then fit a distribution family $\mathcal{F} = \{P_\phi : \phi \in \Phi\}$ separately to each
score vector in order to obtain a differentiable model of each class's score distribution:
%
\begin{equation*}
\hat{\phi}_1 = \operatorname{Fit}_{\mathcal{F}}\!\left(\mathbf{s}^{(1)}\right),
\qquad
\hat{\phi}_0 = \operatorname{Fit}_{\mathcal{F}}\!\left(\mathbf{s}^{(0)}\right)
\end{equation*}
%
Let $\hat{F}_1$ and $\hat{F}_0$ denote the CDFs of $P_{\hat{\phi}_1}$ and $P_{\hat{\phi}_0}$,
respectively. Specific choices of $\mathcal{F}$ (Gaussian, Beta, KDE) are discussed in
\Cref{sec:roll-instantiations}.
Using this notation, at a threshold $\alpha$:
\begin{align}
\operatorname{TPR}(\alpha) &= 1 - \hat{F}_1(\alpha), \label{eq:tpr-cdf} \\
\operatorname{FPR}(\alpha) &= \hat{F}_0(\alpha). \label{eq:fpr-cdf}
\end{align}
Let $\hat{F}_1^{-1}$ and $\hat{F}_0^{-1}$ denote the quantile functions of $P_{\hat{\phi}_1}$ and
$P_{\hat{\phi}_0}$, respectively. The threshold achieving a target TPR of $\alpha$ follows directly
from \eqref{eq:tpr-cdf}:
\[
\operatorname{TPR}(\tau) = \alpha \;\Longrightarrow\; 1 - \hat{F}_1(\tau) = \alpha
\;\Longrightarrow\; \tau = \hat{F}_1^{-1}(1 - \alpha)
\]
and similarly for a target FPR of $\alpha$, from \eqref{eq:fpr-cdf}:
\[
\operatorname{FPR}(\tau) = \alpha \;\Longrightarrow\; \hat{F}_0(\tau) = \alpha
\;\Longrightarrow\; \tau = \hat{F}_0^{-1}(\alpha)
\]
This yields two versions of our objective. Substituting the ICDF threshold into the rate
expressions, our estimated TPR at a fixed FPR of $\alpha$ is $1 - \hat{F}_1(\hat{F}_0^{-1}(\alpha))$.
Therefore, to maximise TPR we minimise:
\begin{equation}
\mathcal{L}_{\text{ROLL-TPR@FPR}}(f_{\theta}(\mathcal{B}) ; \alpha) = \hat{F}_1(\hat{F}_0^{-1}(\alpha))
\label{eq:roll-tpr-at-fpr}
\end{equation}
Our estimated FPR at a fixed TPR of $\alpha$ is $\hat{F}_0(\hat{F}_1^{-1}(1 - \alpha))$.
Therefore, to minimise FPR:
\begin{equation}
\mathcal{L}_{\text{ROLL-FPR@TPR}}(f_{\theta}(\mathbf{X}) ; \alpha) = \hat{F}_0(\hat{F}_1^{-1}(1 - \alpha))
\label{eq:roll-fpr-at-tpr}
\end{equation}
Observe that \eqref{eq:roll-tpr-at-fpr} and \eqref{eq:roll-fpr-at-tpr} share the same
functional structure: each evaluates one empirical CDF at the quantile of the other. This
symmetry means \eqref{eq:roll-fpr-at-tpr} can be reduced to \eqref{eq:roll-tpr-at-fpr} by
negating the model scores and exchanging the class labels. Concretely, replacing
$f_\theta(\mathbf{x})$ with $-f_\theta(\mathbf{x})$ and $y$ with $1 - y$ swaps the roles
of the two populations, turning a TPR@FPR objective into a FPR@TPR one. Because the
two objectives are therefore interchangeable at the implementation level, all derivations
and implementation details below are given for \eqref{eq:roll-tpr-at-fpr} only;
\eqref{eq:roll-fpr-at-tpr} follows by the same transformation applied to the inputs.
\Cref{fig:roll-principle} provides a geometric view of \eqref{eq:roll-tpr-at-fpr}.
\begin{figure}[H]
\centering
\includegraphics[width=0.85\textwidth]{content/method/figures/roll_principle.png}
\caption{Geometric picture of the ROLL TPR@FPR objective. \textit{Upper}: score PDFs
for the negative (blue, $y=0$) and positive (orange, $y=1$) classes. The operating
threshold $\tau$ (dashed) is placed at the $(1-\alpha)$-quantile of the negative
class, so that exactly $\alpha$ fraction of negative scores exceed it ($\operatorname{FPR}=\alpha$).
The shaded area under the positive-class PDF to the left of $\tau$ is
$\hat{F}_1(\tau)$ --- the ROLL loss we minimize. Tick marks along the bottom represent
observed positive-class scores. \textit{Lower}: corresponding CDFs $\hat{F}_0$ (blue)
and $\hat{F}_1$ (orange). Reading horizontally at height $1-\alpha$ on $\hat{F}_0$
gives $\tau$; reading $\hat{F}_1(\tau)$ off $\hat{F}_1$ at the same $\tau$ gives the
loss; the distance from $\hat{F}_1(\tau)$ to $1$ is the resulting TPR.}
\label{fig:roll-principle}
\end{figure}
We now derive the gradient of \eqref{eq:roll-tpr-at-fpr} with respect to $f_\theta(\mathbf{x}_i)$.
If $y_i = 1$, then $(\mathbf{x}_i, y_i) \in \mathcal{B}_1$ and $f_\theta(\mathbf{x}_i)$ has no
effect on $\hat{F}_0^{-1}(\alpha)$. Therefore:
\begin{equation}
\left.\frac{\partial \mathcal{L}_{\text{ROLL-TPR@FPR}}(f_\theta(\mathcal{B}) ; \alpha)}
{\partial f_\theta(\mathbf{x}_i)}\right|_{y_i = 1}
= \frac{\partial \hat{F}_1(\hat{F}_0^{-1}(\alpha))}{\partial f_\theta(\mathbf{x}_i)}
\end{equation}
If $y_i = 0$, then $(\mathbf{x}_i, y_i) \in \mathcal{B}_0$ and $f_\theta(\mathbf{x}_i)$ affects
$\hat{F}_0^{-1}(\alpha)$ instead. Applying the chain rule:
\begin{equation}
\left.\frac{\partial \mathcal{L}_{\text{ROLL-TPR@FPR}}(f_\theta(\mathbf{X}) ; \alpha)}
{\partial f_\theta(\mathbf{x}_i)}\right|_{y_i = 0}
= \frac{\partial \hat{F}_1(\hat{F}_0^{-1}(\alpha))}{\partial \hat{F}_0^{-1}(\alpha)}
\cdot \frac{\partial \hat{F}_0^{-1}(\alpha)}{\partial f_\theta(\mathbf{x}_i)}
\end{equation}
Unifying both cases, the general gradient applicable to all ROLL instantiations is:
\begin{equation}
\frac{\partial \mathcal{L}_{\text{ROLL-TPR@FPR}}(f_\theta(\mathbf{X}) ; \alpha)}
{\partial f_\theta(\mathbf{x}_i)} =
\begin{cases}
\dfrac{\partial \hat{F}_1(\hat{F}_0^{-1}(\alpha))}{\partial f_\theta(\mathbf{x}_i)}
& \text{if } y_i = 1 \\[10pt]
\dfrac{\partial \hat{F}_1(\hat{F}_0^{-1}(\alpha))}{\partial \hat{F}_0^{-1}(\alpha)}
\cdot \dfrac{\partial \hat{F}_0^{-1}(\alpha)}{\partial f_\theta(\mathbf{x}_i)}
& \text{if } y_i = 0
\end{cases}
\label{eq:roll-gradient}
\end{equation}
The concrete form of each partial derivative depends on the choice of distribution family
$\mathcal{F}$, and is derived for each instantiation in \Cref{sec:roll-instantiations}.
%------------------------------------------------
\section{Distribution Instantiations}
\label{sec:roll-instantiations}
% Each subsection below instantiates the general ROLL framework with a specific
% distributional assumption. Gaussian and Beta are included as parametric
% special cases; KDE is the primary contribution and receives the fullest treatment.
\subsection{Gaussian ROLL}
\label{sec:roll-gaussian}
% Assume positive- and negative-class scores are Gaussian distributed. The
% TPR/FPR integrals reduce to closed-form expressions in terms of the normal
% CDF, giving a simple and fast-to-compute instantiation of ROLL. Note the
% limitations: the Gaussian assumption can be violated in practice, motivating
% the more flexible instantiations below.
A simple method of implementing the ROLL loss is to approximate the distribution
of the score distributions as gaussian. This method is simple, as there exist
simple closed-form estimations for the parameters of gaussian distributions, and
they are differentiable.
While the distributions of scores cannot be gauranteed to follow any form, there are
some instances where the score distributions might be similar or exactly gaussian.
For example, when finding linear separability between gaussian distributions, the resulting
scores will be distributed normally, and this method will fit exactly.
\subsubsection{Forward}
\label{sec:roll-gaussian-forward}
% Instantiate \eqref{eq:roll-tpr-at-fpr} under Gaussian score distributions.
% Each class's scores are modelled as $\mathcal{N}(\mu_k, \sigma_k^2)$, fitted
% by MLE. The CDF is the standard normal CDF $\Phi$, giving a closed-form loss.
Firstly, we must estimate the parameters to fit the gaussian distributions
to the true and false populations of the score output.
\begin{equation}
\mu_0 = \frac{1}{|\mathcal{B}_0|}\sum_{\mathbf{x}_j \in \mathcal{B}_0}{f_{\theta}(\mathbf{x}_j)}
\qquad
\mu_1 = \frac{1}{|\mathcal{B}_1|}\sum_{\mathbf{x}_j \in \mathcal{B}_1}{f_{\theta}(\mathbf{x}_j)}
\end{equation}
\begin{equation}
\sigma_0^2 = \frac{1}{|\mathcal{B}_0|}\sum_{\mathbf{x}_j \in \mathcal{B}_0}(f_\theta(\mathbf{x}_j) - \mu_0)^2
\qquad
\sigma_1^2 = \frac{1}{|\mathcal{B}_1|}\sum_{\mathbf{x}_j \in \mathcal{B}_1}(f_\theta(\mathbf{x}_j) - \mu_1)^2
\end{equation}
Once these parameters are estimated, we can then define our ROLL objective as required.
The CDF of the normal distribution with mean $\mu$ and standard deviation $\sigma$ is:
\[
F(x) = \frac{1}{2}\left[1 + \text{erf}\!\left(\frac{x - \mu}{\sigma\sqrt{2}}\right)\right]
\]
and its inverse (quantile function) is:
\[
F^{-1}(p) = \mu + \sigma\sqrt{2}\,\text{ierf}(2p - 1)
\]
where $\operatorname{erf}(x) = \frac{2}{\sqrt{\pi}}\int_0^x e^{-t^2}\,\mathrm{d}t$ is the Gauss
error function and $\operatorname{ierf}$ denotes its functional inverse, satisfying
$\operatorname{ierf}(\operatorname{erf}(x)) = x$. Neither has a closed form, but both are
available in standard numerical libraries (e.g.\ \texttt{torch.erf} and
\texttt{torch.erfinv}). Plugging into the general ROLL framework, the threshold
$\tau$ achieving FPR $= \alpha$ and the resulting loss are:
\begin{align}
\tau &= \mu_0 + \sigma_0\sqrt{2}\,\text{ierf}(2\alpha - 1) \nonumber \\
\mathcal{L}_{\text{ROLL-TPR@FPR}}^{\text{GAUSSIAN}}(f_{\theta}(\mathcal{B}) ; \alpha)
&= \frac{1}{2}\left[1 + \text{erf}\!\left(\frac{\tau - \mu_1}{\sigma_1\sqrt{2}}\right)\right]
\label{eq:roll-tpr-at-fpr-gaussian}
\end{align}
\subsubsection{Gradient Derivation}
\label{sec:roll-gaussian-backward}
% Differentiate the Gaussian loss with respect to $f_\theta(\mathbf{x}_i)$.
% Gradients flow through $\mu_k$ and $\sigma_k^2$ (which are differentiable
% functions of the scores), and through $\Phi$ and $\Phi^{-1}$, which have
% simple closed-form derivatives.
In order to compute the gradient derivation for the gaussian estimates we must first compute the gradient to each of the parameters.
Firstly, for $y_i = 1$.
\[
\dfrac{\partial \hat{F}_1(\hat{F}_0^{-1}(\alpha))}{\partial \mu_1} =
\frac{\partial \frac{1}{2}\left[1 + \text{erf}\!\left(\frac{\hat{F}_0^{-1}(\alpha) - \mu_1}{\sigma_1\sqrt{2}}\right)\right]}{\partial \mu_1}
\]
\[
= -\frac{1}{\sigma_1\sqrt{2\pi}}\exp\!\left(-\frac{(\hat{F}_0^{-1}(\alpha) - \mu_1)^2}{2\sigma_1^2}\right)
\]
\[
\frac{\partial \mu_1}{\partial f_{\theta}(\mathbf{x}_i)} = \frac{1}{|\mathcal{B}_1|}
\]
\[
\dfrac{\partial \hat{F}_1(\hat{F}_0^{-1}(\alpha))}{\partial \sigma_1} =
\frac{\partial \frac{1}{2}\left[1 + \text{erf}\!\left(\frac{\hat{F}_0^{-1}(\alpha) - \mu_1}{\sigma_1\sqrt{2}}\right)\right]}{\partial \sigma_1}
\]
\[
= \frac{\mu_1 - \hat{F}_0^{-1}(\alpha)}{\sigma_1^2\sqrt{2\pi}} \exp\!\left( -\frac{(\hat{F}_0^{-1}(\alpha) - \mu_1)^2}{2\sigma_1^2} \right)
\]
\[
\frac{\partial \sigma_1}{\partial f_{\theta}(\mathbf{x}_i)} = \frac{f_\theta(\mathbf{x}_i) - \mu_1}{|\mathcal{B}_1| \sigma_1}
\]
Combining via the chain rule, the gradient for $y_i = 1$ is:
\begin{equation}
\dfrac{\partial \hat{F}_1(\hat{F}_0^{-1}(\alpha))}{\partial f_{\theta}(\mathbf{x}_i)}
= -\frac{1}{|\mathcal{B}_1|\,\sigma_1\sqrt{2\pi}}
\exp\!\left(-\frac{(\tau - \mu_1)^2}{2\sigma_1^2}\right)
\left(1 + \frac{(\tau - \mu_1)(f_\theta(\mathbf{x}_i) - \mu_1)}{\sigma_1^2}\right)
\label{eq:gauss-grad-y1}
\end{equation}
where $\tau = \hat{F}_0^{-1}(\alpha) = \mu_0 + \sigma_0\sqrt{2}\,\operatorname{ierf}(2\alpha - 1)$.
Given $y_i = 0$, the score $f_\theta(\mathbf{x}_i)$ affects the loss only through the threshold
$\tau = \hat{F}_0^{-1}(\alpha)$, which depends on $\mu_0$ and $\sigma_0$.
Applying the chain rule:
\[
\frac{\partial \hat{F}_1(\tau)}{\partial f_\theta(\mathbf{x}_i)}
= \underbrace{\frac{\partial \hat{F}_1(\tau)}{\partial \tau}}_{\text{PDF of class 1 at }\tau}
\cdot \frac{\partial \tau}{\partial f_\theta(\mathbf{x}_i)}
\]
The first factor is the Gaussian PDF evaluated at $\tau$:
\[
\frac{\partial \hat{F}_1(\tau)}{\partial \tau}
= \frac{1}{\sigma_1\sqrt{2\pi}}\exp\!\left(-\frac{(\tau - \mu_1)^2}{2\sigma_1^2}\right)
\]
For the second factor, $\tau = \mu_0 + \sigma_0\sqrt{2}\,\operatorname{ierf}(2\alpha - 1)$, so
$\frac{\partial \tau}{\partial \mu_0} = 1$ and
$\frac{\partial \tau}{\partial \sigma_0} = \sqrt{2}\,\operatorname{ierf}(2\alpha - 1) = \frac{\tau - \mu_0}{\sigma_0}$.
Combined with $\frac{\partial \mu_0}{\partial f_\theta(\mathbf{x}_i)} = \frac{1}{|\mathcal{B}_0|}$ and
$\frac{\partial \sigma_0}{\partial f_\theta(\mathbf{x}_i)} = \frac{f_\theta(\mathbf{x}_i) - \mu_0}{|\mathcal{B}_0|\sigma_0}$:
\[
\frac{\partial \tau}{\partial f_\theta(\mathbf{x}_i)}
= \frac{1}{|\mathcal{B}_0|}\left(1 + \frac{(\tau - \mu_0)(f_\theta(\mathbf{x}_i) - \mu_0)}{\sigma_0^2}\right)
\]
Therefore, the gradient for $y_i = 0$ is:
\begin{equation}
\dfrac{\partial \hat{F}_1(\hat{F}_0^{-1}(\alpha))}{\partial f_\theta(\mathbf{x}_i)}
= \frac{1}{|\mathcal{B}_0|\,\sigma_1\sqrt{2\pi}}
\exp\!\left(-\frac{(\tau - \mu_1)^2}{2\sigma_1^2}\right)
\left(1 + \frac{(\tau - \mu_0)(f_\theta(\mathbf{x}_i) - \mu_0)}{\sigma_0^2}\right)
\label{eq:gauss-grad-y0}
\end{equation}
\paragraph{Combined gradient.}
Substituting \Cref{eq:gauss-grad-y1} and \Cref{eq:gauss-grad-y0} into \Cref{eq:roll-gradient},
and letting $\tau = \hat{F}_0^{-1}(\alpha) = \mu_0 + \sigma_0\sqrt{2}\,\operatorname{ierf}(2\alpha-1)$:
\begin{equation}
\frac{\partial \mathcal{L}_{\text{ROLL-TPR@FPR}}^{\text{GAUSSIAN}}}{\partial f_\theta(\mathbf{x}_i)} =
\begin{cases}
\dfrac{-1}{|\mathcal{B}_1|\,\sigma_1\sqrt{2\pi}}
\exp\!\left(-\dfrac{(\tau-\mu_1)^2}{2\sigma_1^2}\right)
\!\left(1 + \dfrac{(\tau-\mu_1)(f_\theta(\mathbf{x}_i)-\mu_1)}{\sigma_1^2}\right)
& \text{if } y_i = 1 \\[18pt]
\dfrac{1}{|\mathcal{B}_0|\,\sigma_1\sqrt{2\pi}}
\exp\!\left(-\dfrac{(\tau-\mu_1)^2}{2\sigma_1^2}\right)
\!\left(1 + \dfrac{(\tau-\mu_0)(f_\theta(\mathbf{x}_i)-\mu_0)}{\sigma_0^2}\right)
& \text{if } y_i = 0
\end{cases}
\label{eq:gauss-grad-combined}
\end{equation}
\subsection{Beta ROLL}
\label{sec:roll-beta}
Beta ROLL applies when model scores are constrained to the open interval $(0,1)$,
as is the case when $f_\theta$ uses a sigmoid output activation. The Beta distribution
is supported on $(0,1)$ and its two shape parameters allow it to capture a wide range
of unimodal and skewed score distributions that a Gaussian cannot represent.
\subsubsection{Forward}
\label{sec:roll-beta-forward}
Given score vectors $\mathbf{s}^{(k)} \in (0,1)^{|\mathcal{B}_k|}$ for class $k$, let
$\mathbf{t}^{(k)} = \mathbf{1} - \mathbf{s}^{(k)}$ denote the complementary scores.
For any vector $\mathbf{v}$, write $\bar{v} = \frac{1}{|\mathbf{v}|}\sum_i v_i$.
We estimate the shape parameters using the closed-form log-moment estimator
of~\cite{chen2022beta}: define for each class $k$,
\begin{equation}
D_k = \overline{s^{(k)}\!\ln s^{(k)}} - \bar{s}^{(k)}\cdot\overline{\ln s^{(k)}}
+ \overline{t^{(k)}\!\ln t^{(k)}} - \bar{t}^{(k)}\cdot\overline{\ln t^{(k)}}
\label{eq:beta-denom}
\end{equation}
\begin{equation}
\hat{a}_k = \frac{\bar{s}^{(k)}}{D_k}, \qquad \hat{b}_k = \frac{\bar{t}^{(k)}}{D_k}
\label{eq:beta-params}
\end{equation}
These estimators follow from two simultaneous moment conditions:
$\mathbb{E}[X] = a/(a+b)$ and
$\operatorname{Cov}(X,\ln X) + \operatorname{Cov}(Y,\ln Y) = 1/(a+b)$,
where $Y=1-X\sim\operatorname{Beta}(b,a)$.
Replacing expectations by sample means and solving in closed form yields~\eqref{eq:beta-params}.
The CDF of $\operatorname{Beta}(\hat{a}_k,\hat{b}_k)$ at $x\in(0,1)$ is the
regularised incomplete beta function:
\begin{equation}
\hat{F}_k(x) = I_x(\hat{a}_k,\hat{b}_k)
= \frac{B(x;\,\hat{a}_k,\hat{b}_k)}{B(\hat{a}_k,\hat{b}_k)},
\qquad B(x;\,a,b) = \int_0^x u^{a-1}(1-u)^{b-1}\,\mathrm{d}u
\label{eq:beta-cdf}
\end{equation}
Substituting into the general ROLL framework \eqref{eq:roll-tpr-at-fpr}:
\begin{equation}
\mathcal{L}_{\text{ROLL-TPR@FPR}}^{\text{BETA}}(f_\theta(\mathcal{B});\,\alpha)
= I_{\hat\tau}(\hat{a}_1,\hat{b}_1),
\qquad
\hat\tau = I^{-1}_\alpha(\hat{a}_0,\hat{b}_0)
\label{eq:roll-tpr-at-fpr-beta}
\end{equation}
where $I^{-1}_\alpha(a,b)$ denotes the $\alpha$-quantile of $\operatorname{Beta}(a,b)$.
\subsubsection{Gradient Derivation}
\label{sec:roll-beta-backward}
Unlike the Gaussian case, we do not derive closed-form gradients for Beta ROLL.
The regularised incomplete beta function $I_x(a,b)$ and its inverse are smooth
functions of all their arguments and are available as differentiable operations in
PyTorch (\texttt{torch.special.betainc} and \texttt{torch.special.betaincinv}).
Gradients of \eqref{eq:roll-tpr-at-fpr-beta} with respect to $f_\theta(\mathbf{x}_i)$
are therefore propagated via automatic differentiation through the estimator
\eqref{eq:beta-params}, the CDF evaluation, and the quantile inversion.
Note that the gradient balance property (\Cref{prop:gradient-balance}) does not hold
for Beta ROLL: the sigmoid activation that maps raw scores into $(0,1)$ is nonlinear,
so a uniform additive shift of the raw scores produces a non-uniform shift of the
transformed inputs, violating the translation-invariance assumption required by the
proof.
\subsection{KDE ROLL}
\label{sec:roll-kde}
% KDE background: see \Cref{sec:kde-background} in \Cref{chapter:related_work}.
% The primary instantiation: non-parametric, no assumption on score shape.
% Discuss the trade-off between flexibility and computational cost, and how
% the bandwidth (see \Cref{sec:kde-bandwidth}) controls bias-variance.
Unfortunately, the score output of trained models remains unpredictable. A solution for this is to use a Kernel Density Estimation (KDE) for probability
estimation. This benefits from being far more adaptable to real-world distributions of outputs of our trained models. The tradeoff of course is computational
complexity - computing the gradient derivation in such a case is non-trivial.
\subsubsection{Forward}
\label{sec:roll-kde-forward}
% Instantiate \eqref{eq:roll-tpr-at-fpr} under KDE. The CDF becomes a sum of
% kernel CDFs; the ICDF is computed numerically. Derive the loss expression
% explicitly and discuss computational cost relative to the parametric cases.
As before, to calculate the loss, we calculate:
\[
\mathcal{L}_{\text{ROLL-TPR@FPR}}^{\text{KDE}}(f_{\theta}(\mathcal{B}) ; \alpha) = \hat{F}_1(\hat{F}_0^{-1}(\alpha))
\]
To calculate this, we first define the kernel function as $\sigma '$, and the CDF of the kernel function as $\sigma$. The approach used here is to define $\sigma$ as a sigmoid function,.
\begin{equation}
\sigma(x ; v) = \left( 1 + \exp(-vx)\right)^{-1} \qquad \sigma'(x ; v) = \frac{v \cdot \exp(-v \cdot x)}{(1 + \exp(-vx))^2}
\end{equation}
Where $v$ is the bandwidth parameter. This affects the calculation. We calculate $v_0$ as the bandwidth for the false population and $v_1$ for the true population. Bandwidth calculation
will be discussed later.
We then denote $\sigma_1(x) = \sigma(x ; v_1), \sigma_0(x) = \sigma(x ; v_0)$ to differentiate between the kernel functions for the true and false population, which only differentiate by bandwidth.
The CDF of the KDE function then becomes:
\[
\hat{F}_{KDE}(\tau ; \mathbf{X}) = \frac{1}{|X|}\sum_i \sigma(\tau - x_i)
\]
While the CDF has a nice formula, the inverse CDF has no known closed formula. It instead must be calculated numerically. Thus, the calculation of the decision threshold $\tau$ must be achevied
using a numerical calculation algorithm. In this work, the Newton--Raphson method~\cite{press2007numerical} is chosen.
Firstly we must calculate $\hat{F}_0^{-1}(\alpha)$ = \tau. We start with our initial guess, $\tau_0$, and for each step, calculate:
\[
\tau_{n+1} = \tau_n - \frac{\hat{F}_0(\tau_n ; \mathcal{X}_0) - \alpha}{\frac{\partial \hat{F}_0(\tau_n ; \mathcal{X}_0)}{\partial \tau}}
\]
This step is taken until $\hat{F}_0(\tau_n ; \mathcal{B}_0) \approx \alpha$ within some acceptable error (in our case $1e-4$).
Our derivative w.r.t. $\tau$ is
\[ %
\frac{\partial \hat{F}_0(\tau_n ; \mathcal{X}_0)}{\partial \tau} = \frac{1}{|\mathcal{X}_0|}\sum_i \sigma'(\tau - x_i)
\]
At this point, we have succesfully calculated $\tau$ for which $\tau = \hat{F}_0^{-1}(\alpha ; \mathbf{X}_0)$.
We can now relatively easily calculate $\hat{F}_1(\tau ; \mathbf{X}_1)$ using the formulat for $\hat{F}_{KDE}(\tau ; \mathbf{X})$ above.
Importantly, we must save the calculation of $\tau$ for the backward derivation, as well as the calculated bandwidths $v_1, v_0$.
\subsubsection{Gradient Derivation}
\label{sec:roll-kde-backward}
% The primary technical contribution. Standard autograd cannot differentiate
% through the numerical ICDF inversion; describe the custom backward pass.
% Gradients flow through the kernel evaluations back to the model scores.
% This is the section that connects to \Cref{sec:roll-backward} motivation.
We derive the concrete form of the two terms from \Cref{eq:roll-gradient} for the KDE instantiation,
where $\tau = \hat{F}_0^{-1}(\alpha ; \mathcal{B}_0)$.
\paragraph{Case $y_i = 1$ (positive class).}
When $\mathbf{x}_i \in \mathcal{B}_1$, the score $f_\theta(\mathbf{x}_i)$ appears directly in $\hat{F}_1$
and has no effect on $\tau$ (which depends only on $\mathcal{B}_0$). Therefore:
\[
\dfrac{\partial \hat{F}_1(\hat{F}_0^{-1}(\alpha))}{\partial f_{\theta}(\mathbf{x}_i)} =
\frac{\partial \frac{1}{|\mathcal{B}_1|}\sum_{\mathbf{x}_j \in \mathcal{B}_1} \sigma_1(\tau - f_{\theta}(\mathbf{x}_j))}{\partial f_{\theta}(\mathbf{x}_i)}
\]
\[
= \frac{\partial \left(\sum_{\mathbf{x}_j \in \mathcal{B}_1 ; j \neq i} \frac{1}{|\mathcal{B}_1|}\sigma_1(\tau - f_{\theta}(\mathbf{x}_j)) + \frac{1}{|\mathcal{B}_1|}\sigma_1(\tau - f_{\theta}(\mathbf{x}_i))\right)}{\partial f_{\theta}(\mathbf{x}_i)}
\]
Applying the chain rule to the only $f_\theta(\mathbf{x}_i)$-dependent term,
$\frac{\partial}{\partial f_\theta(\mathbf{x}_i)}(\tau - f_\theta(\mathbf{x}_i)) = -1$:
\begin{equation}
\frac{\partial \hat{F}_1(\hat{F}_0^{-1}(\alpha))}{\partial f_{\theta}(\mathbf{x}_i)} = -\frac{1}{|\mathcal{B}_1|}\sigma_1'(\tau - f_{\theta}(\mathbf{x}_i))
\label{eq:kde-grad-y1}
\end{equation}
\paragraph{Case $y_i = 0$ (negative class).}
When $\mathbf{x}_i \in \mathcal{B}_0$, the score $f_\theta(\mathbf{x}_i)$ does not appear directly in
$\hat{F}_1$, but influences it through the threshold $\tau$. Per \Cref{eq:roll-gradient}, we compute
each factor of the chain rule separately.
\medskip\noindent\textit{Derivative of $\hat{F}_1$ with respect to $\tau$.}
\[
\frac{\partial \hat{F}_1(\hat{F}_0^{-1}(\alpha))}{\partial \hat{F}_0^{-1}(\alpha)} = \frac{\partial \frac{1}{|\mathcal{B}_1|}\sum_{\mathbf{x}_j \in \mathcal{B}_1} \sigma_1(\tau - f_{\theta}(\mathbf{x}_j))}{\partial \tau}
\]
Applying the chain rule, $\frac{\partial}{\partial\tau}(\tau - f_\theta(\mathbf{x}_j)) = +1$:
\begin{equation}
= +\frac{1}{|\mathcal{B}_1|}\sum_{\mathbf{x}_j \in \mathcal{B}_1}\sigma_1'(\tau - f_{\theta}(\mathbf{x}_j))
\label{eq:kde-dF1-dtau}
\end{equation}
\medskip\noindent\textit{Derivative of $\tau$ with respect to $f_\theta(\mathbf{x}_i)$.}
This requires differentiating a quantity computed numerically --- $\tau = \hat{F}_0^{-1}(\alpha ; \mathcal{B}_0)$
is found via Newton's method rather than in closed form.
We apply the inverse derivative rule:
\[
f'(x) = \frac{1}{f^{-1}'(f(x))}
\]
This allows us to use the derivative of the inverse function instead. We need the inverse of
$\hat{F}_0^{-1}$ viewed as a function of $f_\theta(\mathbf{x}_i)$ --- that is, an expression for
$f_\theta(\mathbf{x}_i)$ in terms of $\tau$, $\alpha$, and $\mathcal{B}_0 \setminus \{\mathbf{x}_i\}$.
Solving $|\mathcal{B}_0|\cdot\alpha = \sum_{\mathbf{x}_j \in \mathcal{B}_0}\sigma_0(\tau - f_\theta(\mathbf{x}_j))$
for $f_\theta(\mathbf{x}_i)$ gives:
\[
f_\theta(\mathbf{x}_i) = \tau - \sigma_0^{-1}\left( |\mathcal{B}_0|\cdot\alpha - \sum_{\mathbf{x}_j \in \mathcal{B}_0 ; j \neq i}\sigma_0(\tau - f_{\theta}(\mathbf{x}_j))\right)
\]
Applying the inverse derivative rule to this inverse:
\[
\frac{\partial \hat{F}_0^{-1}(\alpha ; \mathcal{B}_0)}{\partial f_\theta(\mathbf{x}_i)} = \frac{1}{\frac{\partial \left(\tau - \sigma_0^{-1}\left( |\mathcal{B}_0|\cdot \alpha - \sum_{\mathbf{x}_j \in \mathcal{B}_0 ; j \neq i}\sigma_0(\tau - f_{\theta}(\mathbf{x}_j))\right)\right)}{\partial \tau}}
\]
Applying the chain rule:
\[
= \frac{1}{\sigma_0^{-1}{}'\left( |\mathcal{B}_0|\cdot \alpha - \sum_{\mathbf{x}_j \in \mathcal{B}_0 ; j \neq i}\sigma_0(\tau - f_{\theta}(\mathbf{x}_j))\right)\cdot\left( \sum_{\mathbf{x}_j \in \mathcal{B}_0 ; j \neq i}\sigma_0'(\tau - f_{\theta}(\mathbf{x}_j))\right) + 1}
\]
Since $|\mathcal{B}_0|\cdot\alpha = \sum_{\mathbf{x}_j \in \mathcal{B}_0}\sigma_0(\tau - f_\theta(\mathbf{x}_j))$,
the argument of $\sigma_0^{-1}{}'$ simplifies:
\[
|\mathcal{B}_0|\cdot \alpha - \sum_{\mathbf{x}_j \in \mathcal{B}_0 ; j \neq i}\sigma_0(\tau - f_{\theta}(\mathbf{x}_j)) = \sigma_0(\tau - f_{\theta}(\mathbf{x}_i))
\]
Therefore:
\[
\frac{\partial \hat{F}_0^{-1}(\alpha ; \mathcal{B}_0)}{\partial f_\theta(\mathbf{x}_i)} = \frac{1}{\sigma_0^{-1}{}'\!\left( \sigma_0(\tau - f_{\theta}(\mathbf{x}_i)) \right)\cdot\left( \sum_{\mathbf{x}_j \in \mathcal{B}_0 ; j \neq i}\sigma_0'(\tau - f_{\theta}(\mathbf{x}_j))\right) + 1 }
\]
Substituting $\sigma_0^{-1}{}'(\sigma_0(u)) = \frac{1}{\sigma_0'(u)}$ and multiplying numerator and denominator
by $\sigma_0'(\tau - f_{\theta}(\mathbf{x}_i))$:
\begin{equation}
\frac{\partial \hat{F}_0^{-1}(\alpha ; \mathcal{B}_0)}{\partial f_\theta(\mathbf{x}_i)} =
\frac{\sigma_0'(\tau - f_{\theta}(\mathbf{x}_i))}{\displaystyle\sum_{\mathbf{x}_j \in \mathcal{B}_0}\sigma_0'(\tau - f_{\theta}(\mathbf{x}_j))}
\label{eq:kde-grad-tau}
\end{equation}
\medskip\noindent\textit{Efficient computation via the sigmoid identity.}
Direct evaluation of \eqref{eq:kde-grad-tau} requires computing $\sigma_0'$ for every
point, then dividing by the per-point value — potentially unstable when a point is far
from $\tau$. For the sigmoid kernel, we can avoid this by exploiting the identity
\begin{equation}
\sigma(u;\,v)\,\bigl(1 - \sigma(u;\,v)\bigr) = \frac{\sigma'(u;\,v)}{v}
\label{eq:sigmoid-identity}
\end{equation}
which holds for all $u \in \mathbb{R}$ (proved by direct substitution: both sides equal
$\exp(-v|u|)/(1+\exp(-v|u|))^2$). Applying this with $u = \tau - f_\theta(\mathbf{x}_j)$:
\[
\sigma_0'(\tau - f_\theta(\mathbf{x}_j)) = v_0\,\sigma_0(\tau - f_\theta(\mathbf{x}_j))
\,\bigl(1-\sigma_0(\tau - f_\theta(\mathbf{x}_j))\bigr)
\]
Since $v_0$ cancels between numerator and denominator of \eqref{eq:kde-grad-tau},
the ratio can be written entirely in terms of the sigmoid values already computed
during the forward pass:
\begin{equation}
\frac{\partial \hat{F}_0^{-1}(\alpha ; \mathcal{B}_0)}{\partial f_\theta(\mathbf{x}_i)}
=
\frac{\sigma_0(\tau-f_\theta(\mathbf{x}_i))\,\bigl(1-\sigma_0(\tau-f_\theta(\mathbf{x}_i))\bigr)}
{\displaystyle\sum_{\mathbf{x}_j \in \mathcal{B}_0}
\sigma_0(\tau-f_\theta(\mathbf{x}_j))\,\bigl(1-\sigma_0(\tau-f_\theta(\mathbf{x}_j))\bigr)}
\label{eq:kde-grad-tau-efficient}
\end{equation}
This avoids recomputing $\sigma_0'$ (which requires an extra exponential), instead
reusing the sigmoid activations cached from the CDF computation.
\paragraph{Combined gradient.}
Substituting \Cref{eq:kde-grad-y1}, \Cref{eq:kde-dF1-dtau}, and \Cref{eq:kde-grad-tau}
into \Cref{eq:roll-gradient}:
\begin{equation}
\frac{\partial \mathcal{L}_{\text{ROLL-TPR@FPR}}}{\partial f_\theta(\mathbf{x}_i)} =
\begin{cases}
-\dfrac{1}{|\mathcal{B}_1|}\,\sigma_1'(\tau - f_{\theta}(\mathbf{x}_i))
& \text{if } y_i = 1 \\[14pt]
+\dfrac{1}{|\mathcal{B}_1|}\!\displaystyle\sum_{\mathbf{x}_j \in \mathcal{B}_1}\!\sigma_1'(\tau - f_{\theta}(\mathbf{x}_j))
\;\cdot\;
\dfrac{\sigma_0'(\tau - f_{\theta}(\mathbf{x}_i))}{\displaystyle\sum_{\mathbf{x}_j \in \mathcal{B}_0}\sigma_0'(\tau - f_{\theta}(\mathbf{x}_j))}
& \text{if } y_i = 0
\end{cases}
\label{eq:kde-grad-combined}
\end{equation}
where $\tau = \hat{F}_0^{-1}(\alpha ; \mathcal{B}_0)$.
\subsubsection{Sigmoid Kernel: Explicit Equations}
\label{sec:roll-kde-sigmoid}
The derivations above hold for any differentiable CDF kernel $\sigma$.
In our implementation we choose the logistic sigmoid:
%
\begin{equation}
\sigma(x ; v) = \frac{1}{1 + \exp(-vx)},
\qquad
\sigma'(x ; v) = \frac{v\,\exp(-v|x|)}{\bigl(1 + \exp(-v|x|)\bigr)^{2}}
\label{eq:sigmoid-kernel}
\end{equation}
%
where $v > 0$ is the bandwidth parameter. The derivative formula evaluates at
$|x|$ rather than $x$: since $\sigma'$ is an even function
($\sigma'(x) = \sigma'(-x)$ for any symmetric logistic), both forms are
mathematically identical, but $\exp(-v|x|) \to 0$ as $|x| \to \infty$ while
$\exp(-vx) \to \infty$ for $x \to -\infty$, so the $|x|$ form avoids
floating-point overflow when scores fall well below the threshold.
\paragraph{Forward pass.}
Substituting \eqref{eq:sigmoid-kernel} into the KDE CDF:
%
\begin{equation}
\hat{F}_{k}(\tau ; \mathbf{s}^{(k)}) =
\frac{1}{|\mathcal{B}_k|}
\sum_{\mathbf{x}_j \in \mathcal{B}_k}
\frac{1}{1 + \exp\!\bigl(-v_k\bigl(\tau - f_\theta(\mathbf{x}_j)\bigr)\bigr)}
\label{eq:kde-sigmoid-cdf}
\end{equation}
%
and its derivative with respect to $\tau$ — used both in the Newton--Raphson
inversion and in the gradient computation below — is:
%
\begin{equation}
\frac{\partial \hat{F}_{k}(\tau)}{\partial \tau}
=
\frac{1}{|\mathcal{B}_k|}
\sum_{\mathbf{x}_j \in \mathcal{B}_k}
\frac{v_k\,\exp\!\bigl(-v_k\bigl|\tau - f_\theta(\mathbf{x}_j)\bigr|\bigr)}
{\Bigl(1 + \exp\!\bigl(-v_k\bigl|\tau - f_\theta(\mathbf{x}_j)\bigr|\bigr)\Bigr)^{2}}
\label{eq:kde-sigmoid-pdf}
\end{equation}
\paragraph{Gradient.}
Define the sigmoid kernel PDF shorthand:
\begin{equation}
\varphi_k(u) \;=\;
\frac{v_k\,\exp(-v_k|u|)}{\bigl(1+\exp(-v_k|u|)\bigr)^{2}}
\label{eq:sigmoid-phi}
\end{equation}
This is simply $\sigma'(u ; v_k)$ from \eqref{eq:sigmoid-kernel}, a bell-shaped function
centred at zero with width controlled by bandwidth $v_k$.
Substituting \eqref{eq:sigmoid-kernel} into \eqref{eq:kde-grad-combined} and using this
notation, the explicit sigmoid-kernel gradient is:
\begin{equation}
\frac{\partial \mathcal{L}_{\text{ROLL-TPR@FPR}}}{\partial f_\theta(\mathbf{x}_i)} =
\begin{cases}
-\dfrac{\varphi_1\!\left(\tau - f_\theta(\mathbf{x}_i)\right)}{|\mathcal{B}_1|}
& \text{if } y_i = 1\\[14pt]
+\,\dfrac{\displaystyle\sum_{\mathbf{x}_j \in \mathcal{B}_1}
\varphi_1\!\left(\tau - f_\theta(\mathbf{x}_j)\right)}{|\mathcal{B}_1|}
\;\cdot\;
\dfrac{\varphi_0\!\left(\tau-f_\theta(\mathbf{x}_i)\right)}
{\displaystyle\sum_{\mathbf{x}_j \in \mathcal{B}_0}
\varphi_0\!\left(\tau-f_\theta(\mathbf{x}_j)\right)}
& \text{if } y_i = 0
\end{cases}
\label{eq:kde-sigmoid-grad}
\end{equation}
where $\tau = \hat{F}_0^{-1}(\alpha ; \mathcal{B}_0)$.
The signs reflect the direction of each class's effect on the loss: raising a
positive-class score pushes $\hat{F}_1(\tau)$ down (gradient negative), while
raising a negative-class score shifts $\tau$ right, also lowering $\hat{F}_1(\tau)$
(gradient positive). Positive-class samples are weighted by their individual kernel
PDF at $\tau$; negative-class samples are weighted by the mean positive-class kernel
PDF scaled by each sample's share of the negative-class kernel mass.
%------------------------------------------------
\subsection{Properties of the ROLL Formulation}
\label{sec:roll-properties}
% State and prove (or sketch) key properties: consistency of the estimator as
% the number of samples grows, monotonicity of the loss with respect to the
% TPR/FPR trade-off, and behaviour at the constraint boundary. Discuss what
% guarantees ROLL provides and under what assumptions they hold.
% NOTE - these properties are auto-generated, consider suggestions.
A structural consequence of the ROLL formulation is that gradient mass is automatically
balanced between the two classes, regardless of their relative sizes in the batch.
This stands in contrast to standard losses such as cross-entropy, where the total gradient
contribution from each class is proportional to its count, making heavy class imbalance
a direct source of biased updates. Therefor, instead of needing to compensate for this
with loss-weighting per-class, fitting updates are made based equally off the results
of both classes.
\begin{proposition}[Gradient balance]
\label{prop:gradient-balance}
Let $\mathcal{F}$ be any distribution family whose CDF estimator is translation-invariant,
meaning $\hat{F}_k(t\,;\,\mathbf{s}+\delta) = \hat{F}_k(t-\delta\,;\,\mathbf{s})$
for all $\delta \in \mathbb{R}$. Then for any batch $\mathcal{B} = \mathcal{B}_0 \cup \mathcal{B}_1$:
\begin{equation}
\sum_{\mathbf{x}_i \in \mathcal{B}_1}
\frac{\partial \mathcal{L}_{\emph{\text{ROLL-TPR@FPR}}}}{\partial f_\theta(\mathbf{x}_i)}
\;=\;
-\!\sum_{\mathbf{x}_i \in \mathcal{B}_0}
\frac{\partial \mathcal{L}_{\emph{\text{ROLL-TPR@FPR}}}}{\partial f_\theta(\mathbf{x}_i)}
\label{eq:grad-balance}
\end{equation}
with both sides equal in magnitude to $\dfrac{\partial \hat{F}_1(\tau)}{\partial \tau}$,
independently of $|\mathcal{B}_0|$ and $|\mathcal{B}_1|$.
\end{proposition}
\begin{proof}
We treat each class separately via a uniform-shift argument.
\medskip\noindent\textit{Sum over $\mathcal{B}_0$.}
From \Cref{eq:roll-gradient}, for $y_i = 0$:
\[
\frac{\partial \mathcal{L}}{\partial f_\theta(\mathbf{x}_i)}
= \frac{\partial \hat{F}_1(\tau)}{\partial \tau}
\cdot \frac{\partial \tau}{\partial f_\theta(\mathbf{x}_i)}
\]
We claim $\sum_{\mathbf{x}_i \in \mathcal{B}_0}\frac{\partial \tau}{\partial f_\theta(\mathbf{x}_i)} = 1$.
Shift all negative scores uniformly: $f_\theta(\mathbf{x}_i) \to f_\theta(\mathbf{x}_i) + \delta$ for
$\mathbf{x}_i \in \mathcal{B}_0$. By translation invariance, $\hat{F}_0$ simply translates, so
$\tau = \hat{F}_0^{-1}(\alpha)$ shifts by exactly $\delta$. Differentiating with respect to $\delta$ at
$\delta = 0$:
\[
1 = \frac{\mathrm{d}\tau}{\mathrm{d}\delta}\bigg|_{\delta=0}
= \sum_{\mathbf{x}_i \in \mathcal{B}_0} \frac{\partial \tau}{\partial f_\theta(\mathbf{x}_i)}
\]
Summing the $y_i = 0$ gradients therefore gives:
\begin{equation}
\sum_{\mathbf{x}_i \in \mathcal{B}_0}
\frac{\partial \mathcal{L}}{\partial f_\theta(\mathbf{x}_i)}
= \frac{\partial \hat{F}_1(\tau)}{\partial \tau} \cdot 1
= \frac{\partial \hat{F}_1(\tau)}{\partial \tau}
\label{eq:grad-balance-b0}
\end{equation}
\medskip\noindent\textit{Sum over $\mathcal{B}_1$.}
For $y_i = 1$, $\frac{\partial \mathcal{L}}{\partial f_\theta(\mathbf{x}_i)} = \frac{\partial \hat{F}_1(\tau)}{\partial f_\theta(\mathbf{x}_i)}$.
Shifting all positive scores by $\delta$ and applying translation invariance of $\hat{F}_1$:
$\hat{F}_1(\tau\,;\,\mathbf{s}+\delta\mathbf{1}) = \hat{F}_1(\tau-\delta\,;\,\mathbf{s})$.
Differentiating at $\delta = 0$:
\[
\sum_{\mathbf{x}_i \in \mathcal{B}_1}
\frac{\partial \hat{F}_1(\tau)}{\partial f_\theta(\mathbf{x}_i)}
= -\frac{\partial \hat{F}_1(\tau)}{\partial \tau}
\]
Combining with \eqref{eq:grad-balance-b0} gives \eqref{eq:grad-balance}.
\end{proof}
The translation-invariance assumption holds for all distribution families considered in this
work: for Gaussian MLE, a uniform shift in scores moves $\mu_k$ by the same amount while
leaving $\sigma_k$ unchanged, which translates $\hat{F}_k$; for KDE, each kernel term is
evaluated at $\tau - f_\theta(\mathbf{x}_j)$, so a uniform shift in the $f_\theta(\mathbf{x}_j)$
is equivalent to an equal shift in $\tau$. The condition is necessary: without it, a uniform
score shift need not produce a rigid translation of $\hat{F}_k$, so the shift argument in the
proof breaks down and the gradient sums need not be equal. A natural example where this fails
is a family that first maps scores through a nonlinear link function before fitting, such as
passing scores through a sigmoid into $[0,1]$ before fitting a Beta distribution. A uniform
additive shift of the raw scores then produces a non-uniform shift of the transformed inputs,
changing the shape of the fitted distribution rather than merely translating it, and the
balance property no longer holds.
\paragraph{Gradient locality.}
\label{para:gradient-locality}
A second structural consequence of the ROLL formulation is that gradient mass is
\emph{concentrated near the operating threshold $\tau$}. Samples whose scores lie
close to $\tau$ receive disproportionately large gradient updates, while samples
already far on the correct side contribute very little.
This follows directly from \Cref{eq:kde-grad-combined}. For the positive class ($y_i = 1$),
the gradient magnitude is
\[
\left|\frac{\partial \mathcal{L}}{\partial f_\theta(\mathbf{x}_i)}\right|
= \frac{1}{|\mathcal{B}_1|}\,\sigma_1'(f_\theta(\mathbf{x}_i) - \tau),
\]
which is a bell-shaped kernel PDF that peaks when $f_\theta(\mathbf{x}_i) = \tau$ and
decays to zero as the score moves away from the threshold in either direction.
For the negative class ($y_i = 0$), the gradient is proportional to
$\sigma_0'(\tau - f_\theta(\mathbf{x}_i))$ (normalized over $\mathcal{B}_0$), which is
again bell-shaped and centred on $\tau$.
The effect is that each gradient step focuses the model's attention on the region that
directly determines the operating-point performance: samples that could plausibly be
re-ranked relative to the threshold at the current step. Samples comfortably classified
on the correct side contribute negligible signal, not because they are deemed unimportant,
but because their contribution to the CDF estimate at $\tau$ is already saturated.
This stands in contrast to standard losses such as cross-entropy, where the per-sample
gradient is determined by the predicted probability, with no direct dependence on an
operating threshold. Under heavy class imbalance this can cause the dominant class to
``push'' the threshold away from the intended operating point, a problem that gradient
locality avoids by construction.
\begin{observation}[Gradient locality]
\label{obs:gradient-locality}
For KDE-ROLL, the gradient magnitude of sample $\mathbf{x}_i$ is proportional to
the kernel PDF evaluated at the distance between $f_\theta(\mathbf{x}_i)$ and the
operating threshold $\tau$. Samples nearest to $\tau$ receive the largest updates;
samples far from $\tau$ in either direction contribute negligible gradient.
\end{observation}
\Cref{fig:threshold-locality} illustrates this for a simulated two-class setting. The
upper panel shows the score distributions and the operating threshold; the lower panel
shows the gradient magnitude of each sample. Both classes exhibit a clear bell-shaped
concentration of gradient mass around $\tau$, with smooth decay on either side.
\begin{figure}[H]
\centering
\includegraphics[width=0.82\textwidth]{content/method/figures/threshold_weighting.pdf}
\caption{Gradient locality in KDE-ROLL. \textit{Upper}: score distributions for the
negative (blue, $y=0$) and positive (orange, $y=1$) classes; dashed line marks the
operating threshold $\tau$ at FPR $= \alpha$. \textit{Lower}: gradient magnitude
$|\partial\mathcal{L}/\partial f_\theta(\mathbf{x})|$ for each sample. Both classes
show a bell-shaped concentration of gradient mass near $\tau$, with samples far from
the threshold contributing negligible updates. Smooth curves are the kernel PDFs scaled
for reference.}
\label{fig:threshold-locality}
\end{figure}
\subsection{Gradient Computation and the Custom Backward Pass}
\label{sec:roll-backward}
% Because the distribution parameters are themselves functions of model scores
% (not fixed), the backward pass must differentiate through the distribution
% fitting step. Describe the custom autograd implementation: how gradients flow
% from the ROLL loss back through the distribution parameters to the model
% weights. Highlight where standard autograd would fail and what the manual
% backward pass does differently.
%------------------------------------------------
\section{Implementation Considerations}
\label{sec:roll-implementation}
\subsection{Numerical Stability and the Scaling Trick}
\label{sec:roll-numerical-stability}
% The KDE loss involves sums of kernel evaluations that can span many orders of
% magnitude, leading to floating-point underflow or overflow during training.
% Describe the scaling trick that normalizes these sums to a numerically stable
% range without changing the gradient direction, and show that it is equivalent
% to computing in log-space with a max-subtraction stabilizer (analogous to the
% log-sum-exp trick).
\subsection{Bandwidth Estimation for KDE}
\label{sec:kde-bandwidth}
% The KDE bandwidth $h$ controls how smoothly the score distribution is estimated:
% too small and the loss becomes noisy; too large and the gradient signal is
% washed out. Describe the bandwidth estimation strategy used in ROLL (e.g.
% Silverman's rule, a learned or scheduled $h$), justify the choice, and discuss
% sensitivity of results to this hyperparameter.
% NOTE - not silverman's rule, we do use scheduling, but we also use different
% estimator - see impl.
When computing KDE, kernel selection and bandwidth estimation affect the final result greatly. Large bandwidths
don't accurately track the probability distribution, while small kernels result in a very jagged CDF, which for us means very low gradient values unless unless a point is very close to the decision threshold.
Use of the Improved Sheather-Jones (ISJ) method~\cite{botev2010kde} of bandwidth selection worked reliably well at creating very close-tracking, yet smooth CDFs.
\subsubsection{Bandwidth Scaling for KDE}
\label{sec:kde-bandwidth-scaling}
While the ISJ method worked very well for selecting an acceptable bandwidth, this wasn't necessarily wanted during training. Near the start of the training session, scores tended to be all low. This naturally made the selected bandwidth to be too low. This caused very low training performance near the start of a training session, which could result in non-convergence of the model even after significant training epochs elapsed.
The solution was to use a bandwidth scheduler. Starting with large bandwidths resulted in gradients that spread out more evenly across all scores, and less non-convergence runs occurred. The bandwidth scheduler artificially scales up the bandwidth at the beginning to speed up initial training, and then tapers back to restore the intended bandwidth in order to maximize resuslts.
% TODO: algorithm box — full training loop (forward → ROLL loss → custom backward → weight update)
%%% Local Variables:
%%% mode: latex
%%% TeX-master: "../../main"
%%% TeX-engine: xetex
%%% End: