Fixes fixes fixes

content/appendixs/appendix_kde_bisection.tex

@@ -2,15 +2,51 @@ We describe the bisection procedure used to compute $\tau = \hat{F}_0^{-1}(1-\al
 in the KDE-ROLL forward pass (\Cref{sec:roll-kde-forward}).
 
 \paragraph{Bracket initialisation.}
-Because $\sigma(u;\,v) \to 0$ exponentially as $u \to -\infty$ and $\sigma(u;\,v) \to 1$
-as $u \to +\infty$, with tails decaying within roughly $10/v_0$ of the data range, the
-interval
+A fixed margin of $10/v_0 = 10h_0$ beyond the data range,
 \[
-  \mathrm{lo} = \min(\mathbf{s}^{(0)}) - \tfrac{10}{v_0},
+  \mathrm{lo} = \min(\mathbf{s}^{(0)}) - 10h_0,
   \qquad
-  \mathrm{hi} = \max(\mathbf{s}^{(0)}) + \tfrac{10}{v_0}
+  \mathrm{hi} = \max(\mathbf{s}^{(0)}) + 10h_0,
 \]
-is guaranteed to satisfy $\hat{F}_0(\mathrm{lo}) < 1-\alpha < \hat{F}_0(\mathrm{hi})$.
+only guarantees $\hat{F}_0(\mathrm{lo}) < 1-\alpha < \hat{F}_0(\mathrm{hi})$ for
+$1-\alpha \in (\sigma(-10;1),\,\sigma(10;1)) \approx (4.5\times10^{-5},\,1-4.5\times10^{-5})$,
+since that is the range the sigmoid tails have decayed into at distance $10h_0$ — not for
+every $\alpha \in (0,1)$ as a fixed margin might suggest.
+
+Instead we use an exact, target-dependent bracket. For any $p \in (0,1)$, define
+\[
+  \mathrm{lo}(p) = \min(\mathbf{s}^{(0)}) + h_0\log\!\frac{p}{1-p},
+  \qquad
+  \mathrm{hi}(p) = \max(\mathbf{s}^{(0)}) + h_0\log\!\frac{p}{1-p}.
+\]
+
+\begin{proposition}[Exact bracket containment]
+\label{prop:bisection-bracket}
+If $\min(\mathbf{s}^{(0)}) < \max(\mathbf{s}^{(0)})$, then for every $p \in (0,1)$,
+\[
+  \hat{F}_0(\mathrm{lo}(p)) < p < \hat{F}_0(\mathrm{hi}(p)).
+\]
+\end{proposition}
+
+\begin{proof}
+$\sigma(u;v_0)$ is strictly increasing in $u$ (its derivative $v_0\,\sigma(1-\sigma) > 0$).
+Direct substitution, using $v_0 h_0 = 1$, gives
+$\sigma\bigl(\mathrm{lo}(p) - \min(\mathbf{s}^{(0)});\,v_0\bigr)
+  = \sigma\bigl(\log\tfrac{p}{1-p};\,1\bigr)
+  = \bigl(1+\tfrac{1-p}{p}\bigr)^{-1} = p$
+exactly. For any other score $s_j > \min(\mathbf{s}^{(0)})$,
+$\mathrm{lo}(p) - s_j < \mathrm{lo}(p) - \min(\mathbf{s}^{(0)})$, so monotonicity gives
+$\sigma(\mathrm{lo}(p)-s_j;v_0) < p$. Averaging over $\mathcal{B}_0$, with at least one
+score exceeding the minimum, $\hat{F}_0(\mathrm{lo}(p)) < p$ strictly. The symmetric
+argument at $\max(\mathbf{s}^{(0)})$ gives $\hat{F}_0(\mathrm{hi}(p)) > p$.
+\end{proof}
+
+Setting $p = 1-\alpha$ gives a bracket valid for every $\alpha \in (0,1)$, unlike the
+fixed-margin version above. The degenerate case $\min(\mathbf{s}^{(0)}) = \max(\mathbf{s}^{(0)})$
+— a batch with constant negative-class scores — is excluded from
+\Cref{prop:bisection-bracket} above; it also breaks the Silverman bandwidth estimate
+itself ($h_0 = 0$), and is addressed as part of the degenerate-batch handling still
+outstanding for \Cref{sec:kde-bandwidth}.
 
 \paragraph{Bisection iterations.}
 At each step we evaluate the midpoint,
@@ -24,7 +60,40 @@ and update the bracket:
     \mathrm{hi} \leftarrow \tau_{\mathrm{mid}} & \text{otherwise.}
   \end{cases}
 \]
-Iteration continues until $\lvert\hat{F}_0(\tau_{\mathrm{mid}};\,\mathcal{B}_0)-(1-\alpha)\rvert < \varepsilon$
-(with $\varepsilon = 10^{-4}$), or until a maximum iteration count is reached.
+\paragraph{Stopping criterion.}
+A residual check on $\hat{F}_0$ alone, $\lvert\hat{F}_0(\tau_{\mathrm{mid}})-(1-\alpha)\rvert < \varepsilon$,
+certifies that the \emph{threshold} $\tau$ has been located accurately, but not that the
+\emph{loss} $\hat{F}_1(\tau)$ has: if $\hat{F}_0$ is comparatively flat near the target
+(e.g.\ $v_0 \ll v_1$, or the negative-class kernel density is low at $\tau$), a wide range
+of $\tau$ values can satisfy the residual bound while $\hat{F}_1(\tau)$ still varies
+substantially across that range.
+
+Instead, we exploit that $\hat{F}_1$ is itself a CDF — monotonically increasing — and
+that the bracket invariant $\hat{F}_0(\mathrm{lo}) < 1-\alpha < \hat{F}_0(\mathrm{hi})$
+keeps the true root $\tau^\star$ inside $[\mathrm{lo},\mathrm{hi}]$ throughout. Hence
+\[
+  \hat{F}_1(\mathrm{lo}) \;\leq\; \hat{F}_1(\tau^\star) \;\leq\; \hat{F}_1(\mathrm{hi})
+\]
+regardless of exactly where $\tau^\star$ falls in the bracket. Iteration continues until
+\begin{equation}
+  \hat{F}_1(\mathrm{hi};\,\mathcal{B}_1) - \hat{F}_1(\mathrm{lo};\,\mathcal{B}_1) < \varepsilon_L
+  \label{eq:bisection-loss-cert}
+\end{equation}
+(with $\varepsilon_L = 10^{-4}$), or until a maximum iteration count is reached. Unlike a
+Lipschitz-based bound (e.g.\ using $\max_u \sigma'(u;v_1) = v_1/4$ to bound $\hat{F}_1$'s
+slope and stopping purely on bracket width), \eqref{eq:bisection-loss-cert} is exact
+rather than a worst-case estimate, at the cost of two extra $\hat{F}_1$ evaluations per
+iteration.
+
+\paragraph{Why the residual-only rule fails.} Take $\mathcal{B}_0 = \{-30,-25,5,45,60\}$
+with $v_0 = 0.01$ (wide, smooth kernel) and $\alpha=\tfrac12$; the residual-only rule with
+$\varepsilon=10^{-4}$ stops after 8 iterations at $\tau \approx 10.918$, versus the
+converged root $\tau^\star \approx 10.932$, a gap of only $0.014$ in $\tau$. But with a
+sharp positive-class kernel ($v_1 = 1000$) and a cluster of $\mathcal{B}_1$ sitting inside
+that gap, $\hat{F}_1(\tau)$ differs by $0.40$ between the two: the reported loss is wrong
+by nearly half its range, despite the threshold residual being $8\times$ tighter than
+$\varepsilon$. Certifying via \eqref{eq:bisection-loss-cert} instead catches this
+directly, converging to $\hat{F}_1$ accurate to $3\times 10^{-5}$.
+
 The method converges unconditionally and requires no derivative evaluation, in contrast
 to the Newton--Raphson alternative documented in Appendix~\ref{appendix:kde-nr}.

content/appendixs/appendix_kde_grad.tex

@@ -140,3 +140,38 @@ entirely in terms of sigmoid values already cached from the forward pass:
 
 Substituting \eqref{eq:kde-grad-y1}, \eqref{eq:kde-dF1-dtau}, and \eqref{eq:kde-grad-tau}
 into \eqref{eq:roll-gradient} yields \eqref{eq:kde-grad-combined}.
+
+\paragraph{Numerical stability of the efficient form.}
+\label{sec:kde-grad-stability}
+
+\eqref{eq:kde-grad-tau-efficient} is vulnerable in floating point: $\sigma_0(\tau-f_i)$
+rounds to exactly $1$ once its argument $v_0(\tau-f_i)$ exceeds $\approx 17$ in single
+precision, so $1-\sigma_0(\tau-f_i)$ evaluates to exactly $0$.
+$\sigma_0'$ itself stays representable until the argument exceeds $\approx 88$, meaning we are giving
+up orders of magnitude of valid numerical representation. In extreme cases, such as those
+where a class is well-separated from the threshold $\tau$, and multiple scores fall
+over $5$ deviations away from said threshold, it is possible for many multiple scorewise
+gradients to fail, returning a value of $0$ instead of returning a low magnitude, yet usable, gradient.
+
+The ratio being computed is unchanged; only the arithmetic path is. Writing
+$z_i = v_0(\tau-f_\theta(\mathbf{x}_i))$ and using
+$\sigma'(z)=\sigma(z)(1-\sigma(z))$ together with the stable log-sigmoid identity
+$\log\sigma(z) = -\log(1+e^{-z})$:
+\begin{equation}
+  \log\sigma_0'(z_i) = -|z_i| - 2\log\bigl(1+e^{-|z_i|}\bigr) =: \ell_i
+  \label{eq:kde-grad-log-sigma-prime}
+\end{equation}
+(valid for either sign of $z_i$, since $\sigma'$ is even), which is safe to evaluate for
+any $|z_i|$: the only large-magnitude term is the linear $-|z_i|$, while
+$\log(1+e^{-|z_i|}) \in [0,\log 2]$ regardless of $|z_i|$. The normalised weight is then
+recovered by exponentiating exactly once, inside a $\operatorname{softmax}$-style
+normalisation that is itself overflow-safe:
+\begin{equation}
+  \frac{\sigma_0'(z_i)}{\sum_j \sigma_0'(z_j)} = \exp\Bigl(\ell_i - \operatorname*{logsumexp}_j \ell_j\Bigr).
+  \label{eq:kde-grad-softmax}
+\end{equation}
+This has no failure point short of the representable exponent range. The denominator
+can never collapse as it can under \eqref{eq:kde-grad-tau-efficient}.
+\eqref{eq:kde-grad-tau-efficient} is retained above as the clearest statement of the
+derivation; \eqref{eq:kde-grad-log-sigma-prime}, \eqref{eq:kde-grad-softmax} is what is
+implemented.

content/appendixs/appendix_score_normalization.tex

@@ -0,0 +1,143 @@
+This appendix proves \Cref{prop:normalization-invariance} (\Cref{sec:score-normalization-proof}),
+shows numerically that ISJ does not satisfy its scale-equivariance hypothesis
+(\Cref{sec:score-normalization-isj}), and derives the chain-rule term that a frozen
+normalisation constant $\hat\sigma_\mathcal{B}$ omits, together with a concrete instance
+where omitting it is necessary to preserve \Cref{obs:gradient-sign-correctness}
+(\Cref{sec:score-normalization-missing-term}).
+
+\section{Proof of Proposition~\ref{prop:normalization-invariance}}
+\label{sec:score-normalization-proof}
+
+Write $s$ for the raw batch scores, $c := \hat\sigma_\mathcal{B}(s) > 0$, and
+$\tilde{s} = s/c$. Let $v_k^A := 1/\hat{h}(s^{(k)})$ (raw path) and
+$v_k^B := 1/\hat{h}(\tilde{s}^{(k)})$ (normalised path), and let $\tau^A, \tau^B$ solve
+$\hat{F}_0(\tau; s^{(0)}, v_0^A) = 1-\alpha$ and
+$\hat{F}_0(\tau; \tilde{s}^{(0)}, v_0^B) = 1-\alpha$ respectively.
+
+By hypothesis $\hat{h}(\tilde{s}^{(k)}) = \hat{h}(s^{(k)}/c) = \hat{h}(s^{(k)})/c$, so
+\begin{equation}
+  v_k^B = c \cdot v_k^A.
+  \label{eq:sn-v-scale}
+\end{equation}
+Substituting $\tau = \tau^A/c$ into $\hat{F}_0(\tau; \tilde{s}^{(0)}, v_0^B)$:
+\[
+  v_0^B\bigl(\tau^A/c - \tilde{s}_j^{(0)}\bigr)
+  = c v_0^A\Bigl(\frac{\tau^A - s_j^{(0)}}{c}\Bigr)
+  = v_0^A\bigl(\tau^A - s_j^{(0)}\bigr)
+\]
+for every $j$, so $\hat{F}_0(\tau^A/c; \tilde{s}^{(0)}, v_0^B) = \hat{F}_0(\tau^A; s^{(0)}, v_0^A) = 1-\alpha$
+— i.e.\ $\tau^A/c$ solves the normalised-path equation, and by uniqueness of the
+bisection root, $\tau^B = \tau^A/c$ exactly. The same substitution gives, for every $j$,
+\begin{equation}
+  v_k^B\bigl(\tau^B - \tilde{s}_j^{(k)}\bigr) = v_k^A\bigl(\tau^A - s_j^{(k)}\bigr),
+  \label{eq:sn-argument-match}
+\end{equation}
+i.e.\ every sigmoid-kernel argument is identical between the two paths. Consequently
+$\sigma(\cdot)$ evaluates identically in both paths (so the loss value matches:
+$\hat{F}_1(\tau^B;\tilde{s}^{(1)},v_1^B) = \hat{F}_1(\tau^A;s^{(1)},v_1^A)$), while
+$\sigma_k'(u;v) = v\,\sigma(u;v)(1-\sigma(u;v))$ picks up exactly a factor of
+$v_k^B/v_k^A = c$ by \eqref{eq:sn-v-scale}:
+\begin{equation}
+  \sigma_k'\bigl(\tau^B - \tilde{s}_j^{(k)}; v_k^B\bigr)
+  = c \cdot \sigma_k'\bigl(\tau^A - s_j^{(k)}; v_k^A\bigr).
+  \label{eq:sn-derivative-scale}
+\end{equation}
+
+\paragraph{Case $y_i=1$.} By \eqref{eq:kde-grad-y1} and the frozen-$c$ chain rule
+$\partial/\partial s_i = (1/c)\,\partial/\partial\tilde{s}_i$:
+\[
+  \frac{\partial L^B}{\partial s_i}
+  = \frac{1}{c}\left(-\frac{1}{|\mathcal{B}_1|}\sigma_1'\bigl(\tau^B-\tilde{s}_i;v_1^B\bigr)\right)
+  \overset{\eqref{eq:sn-derivative-scale}}{=}
+  \frac{1}{c}\left(-\frac{c}{|\mathcal{B}_1|}\sigma_1'\bigl(\tau^A-s_i;v_1^A\bigr)\right)
+  = \frac{\partial L^A}{\partial s_i}.
+\]
+
+\paragraph{Case $y_i=0$.} By \eqref{eq:kde-grad-combined}, the gradient is
+$\frac{1}{c}\cdot D_1^B \cdot
+\dfrac{\sigma_0'(\tau^B-\tilde{s}_i;v_0^B)}{\sum_j \sigma_0'(\tau^B-\tilde{s}_j;v_0^B)}$,
+where $D_1^B := \frac{1}{|\mathcal{B}_1|}\sum_j\sigma_1'(\tau^B-\tilde{s}_j;v_1^B) = c\cdot D_1^A$
+by \eqref{eq:sn-derivative-scale}. The ratio's numerator and denominator both pick up the
+same factor $c$ (again by \eqref{eq:sn-derivative-scale}), so the ratio itself is
+unchanged:
+\[
+  \frac{\partial L^B}{\partial s_i}
+  = \frac{1}{c}\cdot c D_1^A \cdot
+    \frac{\sigma_0'(\tau^A-s_i;v_0^A)}{\sum_j\sigma_0'(\tau^A-s_j;v_0^A)}
+  = D_1^A \cdot \frac{\sigma_0'(\tau^A-s_i;v_0^A)}{\sum_j\sigma_0'(\tau^A-s_j;v_0^A)}
+  = \frac{\partial L^A}{\partial s_i}.
+\]
+
+In both cases the outer $1/c$ from the chain rule is exactly cancelled by the $c$
+introduced into $v_k^B$, confirming \Cref{prop:normalization-invariance}. (Verified
+numerically to float32 precision — $\max_i|\partial L^A/\partial s_i - \partial L^B/\partial s_i|
+\approx 5.6\times10^{-9}$ on a random batch, using Silverman's rule for $\hat{h}$.)
+
+\section{ISJ Does Not Satisfy the Hypothesis}
+\label{sec:score-normalization-isj}
+
+Silverman's rule satisfies $\hat{h}(a\cdot x)=a\cdot\hat{h}(x)$ exactly, since
+$\hat{h}\propto\hat\sigma(x)$ and $\hat\sigma(ax)=a\hat\sigma(x)$. ISJ, computed via
+\texttt{KDEpy}'s \texttt{improved\_sheather\_jones} (the primary estimator,
+\Cref{sec:kde-bandwidth}), does not: rescaling a fixed sample by $a$ and re-running ISJ
+gives $\hat{h}(ax)/a$ ratios that vary with $a$ rather than staying constant:
+
+\begin{center}
+\begin{tabular}{rrr}
+\toprule
+scale $a$ & $\hat{h}_{\text{Silverman}}(ax)/a$ & $\hat{h}_{\text{ISJ}}(ax)/a$ \\
+\midrule
+$10^{-6}$ & $0.353732$ & $0.001157$ \\
+$10^{-3}$ & $0.353732$ & $0.001155$ \\
+$10^{0}$  & $0.353732$ & $0.116810$ \\
+$10^{3}$  & $0.353732$ & $0.217671$ \\
+$10^{6}$  & $0.353732$ & $0.217671$ \\
+\bottomrule
+\end{tabular}
+\end{center}
+
+Consequently, for a batch where ISJ (rather than the Silverman fallback) is used,
+neither the loss value nor the gradient is invariant to normalisation: a directly
+measured example (random batch, $n{=}10$) gives loss $0.734$ on raw scores versus
+$0.767$ on normalised scores, and a maximum per-sample gradient difference of $0.014$
+— both far from the exact invariance \Cref{prop:normalization-invariance} guarantees
+for Silverman. This is the genuine numerical-stability role of normalisation for the
+method as actually used: not ``preventing overflow via an unchanged, merely-rescaled
+argument'' (\Cref{prop:normalization-invariance} shows that mechanism is a no-op under
+exact scale-equivariance), but pinning ISJ's input to the fixed scale regime — roughly
+unit variance — it behaves consistently at, rather than leaving it to whatever scale the
+model's raw outputs happen to have at a given point in training.
+
+\section{The Missing Chain-Rule Term for a Differentiated $\hat\sigma_\mathcal{B}$}
+\label{sec:score-normalization-missing-term}
+
+\Cref{sec:roll-numerical-stability} treats $c=\hat\sigma_\mathcal{B}$ as frozen:
+gradients with respect to normalised scores are multiplied by $1/c$ alone. Since $c$ is
+itself a function of every score in the batch, $\tilde{s}_j = s_j/c(s_1,\ldots,s_n)$,
+the fully differentiated chain rule is
+\begin{equation}
+  \frac{\partial L}{\partial s_i}
+  = \sum_j \frac{\partial L}{\partial \tilde{s}_j}\frac{\partial \tilde{s}_j}{\partial s_i}
+  = \frac{1}{c}\frac{\partial L}{\partial \tilde{s}_i}
+    \;-\; \frac{1}{c^2}\frac{\partial c}{\partial s_i}\sum_j s_j\, q_j,
+  \qquad q_j := \frac{\partial L}{\partial \tilde{s}_j},
+  \label{eq:sn-missing-term}
+\end{equation}
+using $\partial \tilde{s}_j/\partial s_i = \delta_{ij}/c - (s_j/c^2)\,\partial c/\partial s_i$.
+For the sample standard deviation used in \Cref{sec:roll-numerical-stability},
+$\partial c/\partial s_i = (s_i-\bar{s})/((n-1)c)$. The frozen-$c$ convention sets
+$\partial c/\partial s_i := 0$ by declaration, which is exactly why the implemented
+gradient (\Cref{sec:roll-numerical-stability}) has no such term.
+
+If instead $\hat\sigma_\mathcal{B}$ is differentiated through, the term in
+\eqref{eq:sn-missing-term} does not generally vanish. Unlike the bandwidth case, there
+is no exact-cancellation argument available here: the symmetry that would zero it out
+via Euler's identity ($\sum_j \tilde{s}_j\, q_j = 0$, from $L$ being homogeneous of
+degree $0$ in $\tilde{s}$) requires the bandwidth to be differentiated through the same
+rescaling too — which \Cref{sec:roll-kde-backward} deliberately does not do. Numerically
+(random batch, $n{=}10$, one $y_i{=}1$ sample): the frozen-$c$ gradient for that sample
+is $-0.00181$ — correctly satisfying \Cref{obs:gradient-sign-correctness} ($\leq 0$
+required) — while the fully-differentiated-$c$ gradient for the same sample is
+$+0.00020$, violating it. Freezing $c$ is therefore adopted for the same reason as
+freezing $v_0,v_1$ (Appendix~\ref{appendix:bandwidth-freeze}): the alternative is not
+merely more complex, it can produce gradients of the wrong sign.