396 lines
15 KiB
TeX
396 lines
15 KiB
TeX
This appendix contains the full PyTorch implementation of the ROLL loss
|
|
functions used in the experiments. The code is reproduced from
|
|
\texttt{src/roll.py} in the accompanying implementation repository,
|
|
reformatted for print width and with unrelated losses (baselines,
|
|
Beta-ROLL) omitted; the logic is unchanged.
|
|
|
|
\section{Gaussian ROLL}
|
|
|
|
The Gaussian variant fits a normal distribution to each class's scores and
|
|
evaluates TPR at the requested FPR analytically via the Gaussian CDF/ICDF.
|
|
|
|
\begin{lstlisting}[language=Python]
|
|
import torch
|
|
from functools import partial
|
|
import numpy as np
|
|
|
|
|
|
# -- helpers --
|
|
|
|
def _torch_mean_std(yh):
|
|
# MLE variance (n denominator), not torch.std's Bessel-corrected default.
|
|
yhmean = torch.mean(yh)
|
|
yhstd = torch.std(yh, correction=0)
|
|
return yhmean, yhstd
|
|
|
|
def _torch_normal_fit(yh):
|
|
yhmean, yhstd = _torch_mean_std(yh)
|
|
return torch.distributions.normal.Normal(yhmean, yhstd)
|
|
|
|
|
|
# -- batch splitting --
|
|
|
|
def split_true_false_indeces(yh, y):
|
|
true_indeces = torch.argwhere(y)[:, 0]
|
|
false_indeces = torch.argwhere(torch.logical_not(y))[:, 0]
|
|
return true_indeces, false_indeces
|
|
|
|
def split_true_false(yh, y):
|
|
true_indeces, false_indeces = split_true_false_indeces(yh, y)
|
|
return yh[true_indeces], yh[false_indeces]
|
|
|
|
|
|
# -- Gaussian ROLL --
|
|
|
|
def _roll_loss_from_fpr_internal(yh, y, fpr):
|
|
true_yh, false_yh = split_true_false(yh, y)
|
|
true_normal = _torch_normal_fit(true_yh)
|
|
false_normal = _torch_normal_fit(false_yh)
|
|
return true_normal.cdf(false_normal.icdf(torch.tensor(1 - fpr)))
|
|
|
|
def roll_loss_from_fpr(fpr):
|
|
return partial(_roll_loss_from_fpr_internal, fpr=fpr)
|
|
\end{lstlisting}
|
|
|
|
\section{Kernelized (KDE) ROLL}
|
|
|
|
The KDE variant replaces the parametric Gaussian with a sigmoid-kernel KDE,
|
|
enabling a custom autograd backward pass that directly implements the
|
|
analytic gradients derived in \Cref{appendix:kde-grad}.
|
|
|
|
\begin{lstlisting}[language=Python]
|
|
import torch
|
|
import logging
|
|
from functools import partial
|
|
from torch.autograd import Function
|
|
import numpy as np
|
|
from KDEpy.bw_selection import silvermans_rule, improved_sheather_jones
|
|
|
|
|
|
# -- bandwidth selection --
|
|
|
|
def _bandwidth(scores_1d, max_v=10000.0):
|
|
"""Return 1/h (kernel precision v), clipped to max_v.
|
|
|
|
Tries ISJ; falls back to Silverman on convergence failure
|
|
(common on small or near-constant batches early in training).
|
|
"""
|
|
data = scores_1d[:, np.newaxis]
|
|
try:
|
|
h = float(improved_sheather_jones(data))
|
|
except Exception:
|
|
h = float(silvermans_rule(data))
|
|
return min(1.0 / h, max_v)
|
|
|
|
|
|
# -- FPR grid for AOC objective --
|
|
# 0.01, 0.02, ..., 0.99 (99 values); sum / 100 approximates
|
|
# integral_0^1 (1 - TPR(alpha)) d alpha = AOC = 1 - AUC.
|
|
|
|
_AOC_FPRS = np.arange(0.01, 1.0, 0.01).tolist()
|
|
_AOC_ALPHAS = [1.0 - fpr for fpr in _AOC_FPRS]
|
|
_AOC_DIVISOR = 100
|
|
|
|
|
|
class KernelizedROLLoss(Function):
|
|
"""KDE-ROLL loss with custom backward pass.
|
|
|
|
Accepts a list of alpha values (CDF targets, alpha_k = 1 - FPR_k) and a
|
|
divisor, so one class covers both objectives:
|
|
point-FPR : alphas=[1-fpr], divisor=1
|
|
AOC : alphas=_AOC_ALPHAS, divisor=_AOC_DIVISOR
|
|
|
|
Notation:
|
|
v0, v1 -- KDE bandwidths for B_0 (negative) and B_1 (positive)
|
|
taus -- per-alpha thresholds from inverting F_0_hat
|
|
scores_pos -- positive-class model scores in the current batch
|
|
scores_neg -- negative-class model scores in the current batch
|
|
"""
|
|
|
|
# -- sigmoid kernel primitives --
|
|
|
|
@staticmethod
|
|
def _sigma(v, x):
|
|
"""sigma(x; v) = 1 / (1 + exp(-v*x)) [kernel CDF]."""
|
|
return 1 / (1 + torch.exp(-v * x))
|
|
|
|
@staticmethod
|
|
def _sigma_prime(v, x):
|
|
"""sigma'(x; v) = v*exp(-v|x|) / (1+exp(-v|x|))^2 [kernel PDF].
|
|
|
|
Uses |x| to avoid overflow: exp(-v|x|)->0 for large |x|,
|
|
whereas exp(-v*x)->inf for large negative x.
|
|
sigma' is even so the two forms are equivalent.
|
|
"""
|
|
x = torch.abs(x)
|
|
return (v * torch.exp(-v * x)) / ((1 + torch.exp(-v * x)) ** 2)
|
|
|
|
@staticmethod
|
|
def _log_sigma_prime(v, x):
|
|
"""log sigma'(x; v) = -v|x| - 2*softplus(-v|x|) [stable log-derivative].
|
|
|
|
Used wherever sigma' only ever feeds a normalised ratio (softmax):
|
|
computing sigma*(1-sigma) directly saturates to exactly 0 once
|
|
v|x| exceeds ~17 in fp32 (sigma itself rounds to exactly 1 there),
|
|
while the true derivative stays representable to ~88. This form
|
|
never materialises the raw derivative, so it has no such failure
|
|
point.
|
|
"""
|
|
vx = v * torch.abs(x)
|
|
return -vx - 2 * torch.nn.functional.softplus(-vx)
|
|
|
|
# -- KDE CDF and ICDF --
|
|
|
|
@staticmethod
|
|
def _cdf(v, tau, scores):
|
|
"""F_hat(tau; scores) = (1/n) sum_j sigma(tau - s_j; v)."""
|
|
return torch.mean(KernelizedROLLoss._sigma(v, tau - scores))
|
|
|
|
@staticmethod
|
|
def _icdf(v0, scores_neg, alpha, v1, scores_pos,
|
|
loss_tol=1e-4, max_iters=60):
|
|
"""Find tau s.t. F_hat_0(tau; scores_neg) = alpha via bisection.
|
|
|
|
Bracket: lo/hi = min/max(scores_neg) + h0*log(alpha/(1-alpha)),
|
|
h0 = 1/v0. Gives exact containment F_hat_0(lo) < alpha <
|
|
F_hat_0(hi) for every alpha in (0,1), not just an asymptotic
|
|
range like a fixed +-10/v0 margin.
|
|
|
|
Stopping criterion certifies the LOSS, not just tau: F_hat_1 is
|
|
itself a CDF (monotonic), and the bracket invariant keeps the
|
|
true root tau* in [lo, hi] throughout, so F_hat_1(tau*) lies in
|
|
[F_hat_1(lo), F_hat_1(hi)] regardless of where tau* falls in the
|
|
bracket. Stop once that interval is narrower than loss_tol.
|
|
"""
|
|
if v0 > 10000.:
|
|
logging.warning(
|
|
f'Bandwidth v0={v0:.1f} exceeds maximum 10000, capping')
|
|
v0 = 10000.
|
|
|
|
h0 = 1.0 / v0
|
|
offset = (h0 * torch.log(alpha / (1 - alpha))
|
|
if isinstance(alpha, torch.Tensor)
|
|
else h0 * np.log(alpha / (1 - alpha)))
|
|
lo = scores_neg.min() + offset
|
|
hi = scores_neg.max() + offset
|
|
|
|
for _ in range(max_iters):
|
|
mid = (lo + hi) / 2
|
|
F0_mid = KernelizedROLLoss._cdf(v0, mid, scores_neg)
|
|
if F0_mid < alpha:
|
|
lo = mid
|
|
else:
|
|
hi = mid
|
|
|
|
F1_lo = KernelizedROLLoss._cdf(v1, lo, scores_pos)
|
|
F1_hi = KernelizedROLLoss._cdf(v1, hi, scores_pos)
|
|
if F1_hi - F1_lo < loss_tol:
|
|
break
|
|
|
|
return (lo + hi) / 2
|
|
|
|
@staticmethod
|
|
def _icdf_batch(v0, scores_neg, alphas, v1, scores_pos,
|
|
loss_tol=1e-4, max_iters=60):
|
|
"""Vectorized bisection: find tau_k s.t. F_hat_0(tau_k) = alpha_k
|
|
for all K alphas at once.
|
|
|
|
Each bracket uses the exact per-alpha offset
|
|
h0*log(alpha_k/(1-alpha_k)) (same containment guarantee as
|
|
_icdf) rather than a shared fixed margin. Stopping criterion is
|
|
the same certified-loss bound as _icdf, evaluated for all K
|
|
brackets at once.
|
|
"""
|
|
if float(v0) > 10000.:
|
|
logging.warning(
|
|
f'Bandwidth v0={float(v0):.1f} exceeds maximum 10000, capping')
|
|
v0 = (v0.clamp(max=10000.)
|
|
if isinstance(v0, torch.Tensor) else 10000.)
|
|
|
|
K = len(alphas)
|
|
alphas_t = torch.tensor(alphas,
|
|
dtype=scores_neg.dtype, device=scores_neg.device)
|
|
h0 = 1.0 / float(v0)
|
|
# log(alpha/(1-alpha)) computed in float64 first: casting alpha
|
|
# to the working dtype before subtracting from 1 loses all
|
|
# precision once alpha is within ~1e-7 of 1, collapsing
|
|
# 1-alpha to exactly 0 and the offset to inf.
|
|
alphas_f64 = np.asarray(alphas, dtype=np.float64)
|
|
offsets_np = h0 * np.log(alphas_f64 / (1 - alphas_f64))
|
|
offsets = torch.tensor(offsets_np, dtype=scores_neg.dtype,
|
|
device=scores_neg.device)
|
|
lo = float(scores_neg.min()) + offsets
|
|
hi = float(scores_neg.max()) + offsets
|
|
|
|
for _ in range(max_iters):
|
|
mid = (lo + hi) / 2
|
|
diff0 = mid.unsqueeze(1) - scores_neg.unsqueeze(0) # (K, n0)
|
|
F0_mid = torch.mean(
|
|
KernelizedROLLoss._sigma(v0, diff0), dim=1) # (K,)
|
|
too_low = F0_mid < alphas_t
|
|
lo = torch.where(too_low, mid, lo)
|
|
hi = torch.where(too_low, hi, mid)
|
|
|
|
diff1_lo = lo.unsqueeze(1) - scores_pos.unsqueeze(0) # (K, n1)
|
|
diff1_hi = hi.unsqueeze(1) - scores_pos.unsqueeze(0) # (K, n1)
|
|
F1_lo = torch.mean(
|
|
KernelizedROLLoss._sigma(v1, diff1_lo), dim=1) # (K,)
|
|
F1_hi = torch.mean(
|
|
KernelizedROLLoss._sigma(v1, diff1_hi), dim=1) # (K,)
|
|
if torch.all(F1_hi - F1_lo < loss_tol):
|
|
break
|
|
|
|
return (lo + hi) / 2 # (K,)
|
|
|
|
# -- forward pass --
|
|
|
|
@staticmethod
|
|
def forward(ctx, alphas, divisor, yh, y, gamma=1.0, normalize=True):
|
|
# Normalise scores by std for numerical stability.
|
|
# Chain-rule correction is applied in backward().
|
|
# clamp avoids inf when all batch scores are identical (zero variance).
|
|
if normalize:
|
|
score_scale = 1 / torch.sqrt(
|
|
torch.var(yh, correction=1).clamp(min=1e-8))
|
|
else:
|
|
score_scale = torch.tensor(1.0, device=yh.device)
|
|
yh_norm = yh * score_scale
|
|
|
|
pos_idx, neg_idx = split_true_false_indeces(yh, y)
|
|
scores_pos = yh_norm[pos_idx] # s^(1)
|
|
scores_neg = yh_norm[neg_idx] # s^(0)
|
|
|
|
# Bandwidths via ISJ; gamma > 1 widens kernels early in training.
|
|
v0 = (torch.tensor(
|
|
_bandwidth(scores_neg.cpu().numpy().astype(np.float64)))
|
|
/ gamma).to(yh.device)
|
|
v1 = (torch.tensor(
|
|
_bandwidth(scores_pos.cpu().numpy().astype(np.float64)))
|
|
/ gamma).to(yh.device)
|
|
|
|
# tau_k = F_0_hat^{-1}(alpha_k):
|
|
# threshold where the negative-class CDF equals alpha_k.
|
|
if len(alphas) == 1:
|
|
taus = KernelizedROLLoss._icdf(
|
|
v0, scores_neg, alphas[0], v1, scores_pos).unsqueeze(0)
|
|
else:
|
|
taus = KernelizedROLLoss._icdf_batch(
|
|
v0, scores_neg, alphas, v1, scores_pos) # (K,)
|
|
|
|
# loss_k = F_1_hat(tau_k) = 1 - TPR@FPR_k
|
|
diff_pos = taus.unsqueeze(1) - scores_pos.unsqueeze(0) # (K, n1)
|
|
losses = torch.mean(
|
|
KernelizedROLLoss._sigma(v1, diff_pos), dim=1) # (K,)
|
|
|
|
total_loss = losses.sum() / divisor
|
|
|
|
logging.debug(
|
|
f'KDE-ROLL K={len(alphas)} gamma={gamma:.4f}'
|
|
f' v0={v0:.4f} v1={v1:.4f} loss={total_loss:.4f}')
|
|
ctx.save_for_backward(yh, scores_pos, scores_neg,
|
|
pos_idx, neg_idx, taus, v1, v0,
|
|
score_scale)
|
|
ctx.divisor = divisor
|
|
return total_loss
|
|
|
|
# -- backward pass --
|
|
|
|
@staticmethod
|
|
def backward(ctx, grad_output):
|
|
yh, scores_pos, scores_neg, pos_idx, neg_idx, \
|
|
taus, v1, v0, score_scale = ctx.saved_tensors
|
|
divisor = ctx.divisor
|
|
n1 = scores_pos.shape[0]
|
|
|
|
diff_pos = taus.unsqueeze(1) - scores_pos.unsqueeze(0) # (K, n1)
|
|
sp_pos = KernelizedROLLoss._sigma_prime(v1, diff_pos) # (K, n1)
|
|
|
|
# d(sum_k loss_k)/d(s_i^(1)) = -(1/n1) sum_k sigma'(tau_k - s_i)
|
|
grad_pos = -(1.0 / n1) * sp_pos.sum(dim=0) / divisor # (n1,)
|
|
|
|
# dF_1_hat(tau_k)/d(tau_k) = (1/n1) sum_j sigma'(tau_k - s_j^(1))
|
|
dF1_dtau = sp_pos.mean(dim=1) # (K,)
|
|
|
|
# d(tau_k)/d(s_i^(0)), computed in log-space for numerical
|
|
# stability -- replaces sigma*(1-sigma), which saturates to
|
|
# exactly 0 in fp32 once v0*diff is large, even though the true
|
|
# derivative stays representable much further out.
|
|
diff_neg = taus.unsqueeze(1) - scores_neg.unsqueeze(0) # (K, n0)
|
|
log_kernel = KernelizedROLLoss._log_sigma_prime(
|
|
v0, diff_neg) # (K, n0)
|
|
grad_tau_neg = torch.softmax(log_kernel, dim=1) # (K, n0)
|
|
|
|
# d(sum_k loss_k)/d(s_i^(0))
|
|
grad_neg = (dF1_dtau.unsqueeze(1)
|
|
* grad_tau_neg).sum(dim=0) / divisor # (n0,)
|
|
|
|
grad_input = torch.zeros(yh.shape[0], device=yh.device)
|
|
grad_input[pos_idx] = grad_pos * grad_output
|
|
grad_input[neg_idx] = grad_neg * grad_output
|
|
|
|
if torch.any(torch.isnan(grad_input)):
|
|
logging.error('NaN gradient detected in KDE-ROLL backward')
|
|
|
|
# Chain-rule for score normalisation: yh_norm = yh * score_scale
|
|
return None, None, grad_input * score_scale, None, None, None
|
|
|
|
|
|
# -- public constructors --
|
|
|
|
def _kernelized_roll_fpr_internal(yh, y, cdf_target, gamma=1.0, normalize=True,
|
|
bce_weight=0.0, bce_pos_weight=1.0):
|
|
roll = KernelizedROLLoss.apply([cdf_target], 1, yh, y, gamma, normalize)
|
|
if bce_weight == 0.0:
|
|
return roll
|
|
pos_weight = torch.tensor([bce_pos_weight],
|
|
dtype=yh.dtype, device=yh.device)
|
|
bce = torch.nn.functional.binary_cross_entropy_with_logits(
|
|
yh, y.float(), pos_weight=pos_weight)
|
|
return roll + bce_weight * bce
|
|
|
|
def kernelized_roll_fpr(fpr, normalize=True,
|
|
bce_weight=0.0, bce_pos_weight=1.0):
|
|
# Pass 1-fpr as CDF target: F_0_hat(tau) = 1-fpr iff P(neg > tau) = fpr
|
|
return partial(_kernelized_roll_fpr_internal, cdf_target=1 - fpr,
|
|
normalize=normalize, bce_weight=bce_weight,
|
|
bce_pos_weight=bce_pos_weight)
|
|
|
|
|
|
def _kernelized_roll_tpr_internal(yh, y, tpr, gamma=1.0,
|
|
bce_weight=0.0, bce_pos_weight=1.0):
|
|
# Reduce TPR@FPR to FPR@TPR by negating scores and swapping labels.
|
|
# `1 - yh` (not `-yh`) is used here; safe only because it differs
|
|
# from `-yh` by a constant shift, and KernelizedROLLoss is
|
|
# translation-invariant in its scores.
|
|
roll = KernelizedROLLoss.apply([tpr], 1, 1 - yh, 1 - y, gamma)
|
|
if bce_weight == 0.0:
|
|
return roll
|
|
pos_weight = torch.tensor([bce_pos_weight],
|
|
dtype=yh.dtype, device=yh.device)
|
|
bce = torch.nn.functional.binary_cross_entropy_with_logits(
|
|
yh, y.float(), pos_weight=pos_weight)
|
|
return roll + bce_weight * bce
|
|
|
|
def kernelized_roll_tpr(tpr, bce_weight=0.0, bce_pos_weight=1.0):
|
|
return partial(_kernelized_roll_tpr_internal, tpr=tpr,
|
|
bce_weight=bce_weight, bce_pos_weight=bce_pos_weight)
|
|
|
|
|
|
def _kernelized_roll_aoc_internal(yh, y, gamma=1.0, normalize=True,
|
|
bce_weight=0.0, bce_pos_weight=1.0):
|
|
roll = KernelizedROLLoss.apply(
|
|
_AOC_ALPHAS, _AOC_DIVISOR, yh, y, gamma, normalize)
|
|
if bce_weight == 0.0:
|
|
return roll
|
|
pos_weight = torch.tensor([bce_pos_weight],
|
|
dtype=yh.dtype, device=yh.device)
|
|
bce = torch.nn.functional.binary_cross_entropy_with_logits(
|
|
yh, y.float(), pos_weight=pos_weight)
|
|
return roll + bce_weight * bce
|
|
|
|
def kernelized_roll_aoc(normalize=True,
|
|
bce_weight=0.0, bce_pos_weight=1.0):
|
|
return partial(_kernelized_roll_aoc_internal, normalize=normalize,
|
|
bce_weight=bce_weight, bce_pos_weight=bce_pos_weight)
|
|
\end{lstlisting}
|