###### Abstract

We study the minimal extended seesaw mechanism with one sterile neutrino in a
3+1 framework using modular $S_{3}$ symmetry. The active-sterile neutrino
models are classified based on the assignments of $S_{3}$ representations and
modular weights of the left-handed lepton doublets, triplet right-handed
neutrino, and sterile neutrino. No scalar flavons are considered, and the
flavor symmetry is broken by the vacuum expectation value (vev) of the modulus
$\tau$. For a particular set of representations of the Leptons and Higgs
field, we obtain eleven (11) different models based on different modular
weights of charged lepton $(k_{L})$ and right-handed neutrino ($k_{N}$). Out
of these, we consider two models, which are discriminated by carrying out the
numerical analysis so that the parameter space in each model can fit the
latest neutrino oscillation data at 3$\sigma$. The Planck cosmological bound
on the upper limit of the sum of the active neutrino masses $\sum
m_{i}<0.12$eV is also considered. Finally, the best-fit parameters of the
neutrino observables and model predictions are evaluated using the minimum
$\chi^{2}$ analysis.

Modular $S_{3}$ symmetry, Sterile neutrino, Beyond Standard Model , MES
mechanism.

##  I Introduction

Among various problems in the Standard Model(SM) of particle physics, the
origin of neutrino masses and flavour mixing, mass hierarchy, Dirac CP-
violating phase, octant of $\theta_{23}$, etc. are some of the crucial areas
that scientists are trying to explain both experimentally and theoretically.
In addition, certain discrepancies in the results of LSND and MiniBooNE
compelled the need to extend the three neutrino theory. These extra states of
neutrinos do not have weak interaction; therefore, they are referred to as
$sterile$ neutrinos. However, sterile neutrinos can mix with the active
neutrinos. Such neutrinos are incorporated with the three neutrino theory in
various literatures in a 3+1 scheme, 3+1+1, 3+2, etc. One of the most
straightforward extensions is the 3+1 scheme, where a singlet sterile neutrino
is added to the three active neutrinos, and the sterile neutrino gets mass
through a minimal extended seesaw mechanism(MES)[1].

Neutrino mass models with non-Abelian discrete symmetries have been explored
to address the problems of flavour mixing. However, such approaches require
ad-hoc scalar fields called $flavons$, which are singlet under SM but
transform non-trivially as irreducible representations of the symmetry. After
the electroweak symmetry breaking, these flavons acquire vacuum expectation
values along specific vacuum alignments, which shape the structure of fermion
mass matrices. The main areas for improvement of such models are the presence
of additional new flavons, difficulty in generating the right vacuum
alignments of the scalars, and requirements of heavy fine-tuning of the
parameters. Very recently, an interesting method was proposed in which modular
flavour symmetry is used with the discrete symmetries as its subgroups[2]. In
such frameworks, the Yukawa couplings transform non-trivially as modular
forms. These modular forms depend on a single parameter $\tau$ called the
modulus. Minimal or no scalar flavons are required to break the symmetry as
the fermion masses are generated from the symmetry breaking by the vev of the
modulus $\tau$. The modular symmetry defined by a finite modular group
$\Gamma_{N}$ is isomorphic to the conventional discrete symmetry groups
$S_{3},S_{4},A_{4},A_{5},etc.$. For instance, the simplest group for $N=2,\
\Gamma_{2}$ is isomorphic to $S_{3}$. Detailed analysis of modular groups and
their application in neutrino model building are studied in Refs. [2, 3].
There are numerous studies based on modular groups such as $S_{3}$[4], $A_{4}$
[5, 6] $A_{5}$ [7, 8], $S_{4}$[9, 10], etc.

In modular $S_{3}$ symmetry, the Yukawa coupling transforms as an $S_{3}$
doublet represented by $Y_{2}^{(k)}$. Here, $k$ denotes the weight of the
modular form, and the subscript $2$ denotes the doublet representation in
$S_{3}$. The doublet modular form of lowest weight 2,
$Y_{2}^{(2)}=(y_{1},y_{2})$ can be expressed in terms of the Dedekind-eta
function $\eta(\tau)$ and its derivative $\eta^{\prime}(\tau)$ as

$\displaystyle y_{1}(\tau)=$
$\displaystyle\frac{i}{4\pi}\left(\frac{\eta^{\prime}(\tau/2)}{\eta(\tau/2)}+\frac{\eta^{\prime}((\tau+1)/2)}{\eta((\tau+1)/2)}-\frac{8\eta^{\prime}(2\tau)}{\eta(2\tau)}\right)$
$\displaystyle y_{2}(\tau)=$
$\displaystyle\frac{\sqrt{3}i}{4\pi}\left(\frac{\eta^{\prime}(\tau/2)}{\eta(\tau/2)}-\frac{\eta^{\prime}((\tau+1)/2)}{\eta((\tau+1)/2)}\right)$
(1)

where $\eta(\tau)$ is defined as

$\eta(\tau)=q^{1/24}\prod_{n=1}^{\infty}(1-q^{n}),\ \ \ \ q\equiv
e^{i2\pi\tau}.$ (2)

Multiplication rules of $S_{3}$ can construct other higher-weighted modular
forms. The modular forms of higher weights $k=4,6$ are given as

$\displaystyle Y_{1}^{(4)}=y_{1}^{2}+y_{2}^{2},\ \ \ \ \ \
Y_{1}^{(6)}=3y_{1}^{2}y_{2}-y_{2}^{3},\ \ \ \ \ \ \ \
Y_{1^{\prime}}^{(6)}=y_{1}^{3}-3y_{1}y_{2}^{2},$ $\displaystyle
Y_{2}^{(4)}=\left(\begin{array}[]{ccc}2y_{1}y_{2}\\\
y_{1}^{2}-y_{2}^{2}\end{array}\right),\ \ \ \ \ \
Y_{2}^{(6)}=\left(\begin{array}[]{ccc}y_{1}^{3}+y_{1}y_{2}^{2}\\\
y_{2}^{3}+y_{1}^{2}y_{2}\end{array}\right)$ (7)

It is important to note that a doublet modular form must have weight $k\geq 2$
while the singlet modular form $Y_{1}$ cannot have weight $2$. The lowest
weight of an $S_{3}$ singlet modular form starts from $4$. The number of
independent modular forms of weight k increases with k/2+1. For minimality, we
consider modular forms with a maximum weight of 6.

The present work will analyze active-sterile neutrino mass models based on
modular $S_{3}$ symmetry in a minimal extended seesaw (MES) framework. The SM
particle contents are extended by three right-handed neutrino fields
$N_{i}(i=1,2),N_{3}$, and a singlet sterile neutrino field $S$. We shall
explore different models corresponding to different assignments of the fields’
modular weight $`k$’. We develop the neutrino mass models without using any
hypothetical scalar flavons in the charged lepton as well as the neutrino
sector, which is an advantage of using modular flavor symmetry.

##  II Building strategy of the models

The $S_{3}$ group has a doublet $(2)$ and two singlets $(1,1^{\prime})$ as its
irreducible representations. We extend the SM particle contents by three
right-handed neutrinos where two components $N_{i}$ $(i=1,2)$ and a third
component $N_{3}$ transform as a doublet and a singlet under $S_{3}$. A right-
handed SM singlet sterile neutrino field $S$ with modular weight $k_{S}$ which
transforms as a singlet under $S_{3}$, is also considered. The three
generations of the SM lepton doublet $L$ are also grouped into two irreps. of
$S_{3}$ as doublet $L_{i}$ $(i=1,2)$ and a singlet $L_{3}$. Similarly, the
right-handed charged leptons are represented as singlets $e_{1},e_{2}$ and
$e_{3}$ of $S_{3}$, respectively. The Higgs represented as $H_{u}$ and $H_{d}$
having hypercharges $+1/2$ and $-1/2$ are taken as singlets of $S_{3}$. The
Higgs develop VEVs after electroweak symmetry breaking along $\langle
H_{u}\rangle=\left(0,v_{u}/\sqrt{2}\right)^{T}$ and $\langle
H_{d}\rangle=\left(v_{d}/\sqrt{2},0\right)^{T}$ inducing SM lepton mass terms.

To simplify our analysis and to avoid trivial interaction terms in the
superpotential, we consider that the weights of SM lepton doublets
$L_{i},L_{3}$ are denoted by $k_{L}$, charged lepton singlets
$e_{1},e_{2},e_{3}$ are denoted by $k_{e_{i}},$ ($i=1,2,3$) respectively. The
heavy right-handed neutrino fields $N_{i},N_{3}$ have modular weight $k_{N}$.
It should be noted that $N_{3}$ and $S$ have the exact transformation charges
under the $S_{3}$ symmetry. All the interactions allowed for $N_{3}$ will also
be possible for $S.$ To avoid weak interactions of the sterile neutrino field,
we also assume that the weights of sterile neutrino and the heavy right-handed
neutrino $N_{3}$ are different, i.e., $k_{N}\neq k_{S}.$ We have fixed three
different weights $k_{S}=(1,4,5)$ for $S$ field in our model111We have taken
this choice to avoid Majorana mass terms of the sterile neutrinos of the form
$\sim\Lambda YSS.$. The complete particle contents of the model and their
group charges are given in Table 1.

Table 1: Particle contents of the model and their group charges. $\frac{Fields}{Charges}$ |  $L_{i}$ |  $L_{3}$ |  $e_{1}^{c}$ |  $e_{2}^{c}$ |  $e_{3}^{c}$ |  $H_{u,d}$ |  $N_{i}^{c}$ |  $N_{3}^{c}$ |  $S^{c}$
---|---|---|---|---|---|---|---|---|---
$SU(2)_{L}$ |  2 |  1 |  1 |  1 |  1 |  1 |  1 |  1 |  1
$U(1)_{Y}$ |  $-1/2$ |  $-1/2$ |  1 |  1 |  1 |  $\pm$ 1/2 |  0 |  0 |  0
$S_{3}$ |  2 |  1′ |  1 |  1 |  1′ |  1 |  2 |  1 |  1
$-k_{i}$ |  $k_{L}$ |  $k_{L}$ |  $k_{e_{1}}$ |  $k_{e_{2}}$ |  $k_{e_{3}}$ |  $k_{H}$ |  $k_{N}$ |  $k_{N}$ |  $k_{S}$

The most general superpotential of the charged-lepton sector invariant under
$S_{3}$ symmetry without specific assignments of modular weights is given by

$\displaystyle W_{e}$ $\displaystyle=$ $\displaystyle\
\alpha(Y_{2}^{(k)}H_{u}\bar{L_{i}})e_{1}+\beta(Y_{2}^{(k)}H_{u}\bar{L_{i}})e_{2}+\gamma(Y_{1}^{(k)}H_{u}\bar{L_{3}})e_{3}$
(8) $\displaystyle+$
$\displaystyle\alpha^{\prime}(Y_{1}^{(k)}H_{u}\bar{L_{3}})e_{1}+\beta^{\prime}(Y_{1^{\prime}}^{(k)}H_{u}\bar{L_{3}})e_{2}+\gamma^{\prime}(Y_{2}^{(k)}H_{u}\bar{L_{i}})_{1^{\prime}}e_{3}$

Here, the superscript $(k)$ in the Yukawa couplings denotes the weights of the
modular forms while the subscripts $1,1^{\prime},2$ denote the representation
of the couplings under $S_{3}$ symmetry. The parameters
$\alpha,\beta,\gamma,\alpha^{\prime},\beta^{\prime},$ and $\gamma^{\prime}$
are the dimensionless coefficients which are the free parameters of the
charged lepton sector. Terms in the superpotential are invariant under a
modular symmetry when the overall weight of the superfields in each term is
zero, and the term is invariant under the discrete symmetry. Based on the
assignments of modular weights $k_{i}$ for each field, certain interaction
terms are allowed or disallowed in the superpotential. For a minimal structure
in the charged lepton sector, we adopt the following constraints regarding the
weights of the right-handed charged lepton singlets

$\displaystyle k_{e_{1}}$ $\displaystyle=$ $\displaystyle 4-k_{L}$ (9)
$\displaystyle k_{e_{2}}$ $\displaystyle=$ $\displaystyle 2-k_{L}$ (10)
$\displaystyle k_{e_{3}}$ $\displaystyle=$ $\displaystyle-k_{L}$ (11)

We assume that $k_{L}\geq 0$ and the Higgs field have zero modular weight. As
a consequence of the above constraints, the coefficients $\beta^{\prime}=0$
since there are no singlet modular forms of weight 2. Since doublet modular
forms cannot have weight zero, the coefficient $\gamma^{\prime}=0.$ Finally,
the term $\alpha^{\prime}(Y_{1^{\prime}}^{(k)}H_{u}\bar{L_{3}})e_{1}$ cannot
exist because $Y_{1^{\prime}}^{(4)}=0.$ Thus, only three terms remain in the
charged lepton superpotential. After the electroweak symmetry breaking, the
Higgs gets a vev, and the resulting charged lepton mass matrix is given by

$\displaystyle m_{L}=\frac{v_{u}}{\sqrt{2}}\left(\begin{array}[]{ccc}\alpha
2y_{1}y_{2}&\beta y_{1}&0\\\ \alpha(y_{1}^{2}-y_{2}^{2})&\beta y_{2}&0\\\
0&0&\gamma\end{array}\right)$ (15)

Similarly, in the case of the neutrino sector, the general superpotential is
given by

$\displaystyle W_{\nu}$ $\displaystyle=$ $\displaystyle\
g_{1}Y_{1}^{(k)}H_{d}\bar{L_{i}}N_{i}+g_{2}Y_{2}^{(k)}H_{d}\bar{L_{i}}N_{3}+g_{3}Y_{2}^{(k)}H_{d}\bar{L_{3}}N_{i}+g_{4}Y_{2}^{(k)}H_{d}\bar{L_{i}}N_{i}$
(16) $\displaystyle+$ $\displaystyle
g_{5}Y_{1^{\prime}}^{(k)}H_{d}\bar{L_{3}}N_{3}+\lambda_{1}\left(Y_{1}^{(k)}\bar{N}^{c}_{i}N_{i}+Y_{2}^{(k)}\bar{N_{i}}^{c}N_{i}\right)$
$\displaystyle+$
$\displaystyle\lambda_{2}\left(Y_{2}^{(k)}\bar{N_{i}}^{c}N_{3}+Y_{2}^{(k)}\bar{N_{3}}^{c}N_{i}\right)+\lambda_{3}Y_{1}^{(k)}\bar{N_{3}}^{c}N_{3}$
$\displaystyle+$
$\displaystyle\lambda_{s1}Y_{1}^{(k)}\bar{S}^{c}N_{3}+\lambda_{s2}Y_{2}^{(k)}\bar{S}^{c}N_{i},$

where the coefficients $g_{i}$ ($i=1,2,3,4,5$) are the dimensionless
coefficients. Other parameters $\lambda_{i}\ (i=1,2,3)$, $\lambda_{sj}\
(j=1,2)$ are the constants having dimension of mass, and they define the mass
scales of the heavy right-handed Majorana neutrino and sterile neutrino
respectively. We have divided our models based on the choice of $k_{L}$ as

  1. 1.

For $k_{L}=0,$ there are two possible models : $A1$ ($k_{N}=0,k_{S}=4$) and
$A2$ ($k_{N}=2,k_{S}=4$). In model $A1$, the vanishing coefficients in eq.(16)
are $g_{2}=g_{3}=g_{4}=0$ and $\lambda_{2}=0,$ whereas, in model $A2,$ the
coefficients $g_{1}=g_{5}=0$.

  2. 2.

For $k_{L}=1,$ the two models are : $B1\ (k_{N}=1,k_{S}=5)$ and $B2\
(k_{N}=3,k_{S}=1)$. In model $B1$, the vanishing coefficients in eq.(16) are
$g_{1}=g_{5}=0$ and $\lambda_{3}=0$ while in model $B2,$ the vanishing
coefficient is $g_{5}=0$.

  3. 3.

For $k_{L}=2,$ two possible models are : $C1\ (k_{N}=0,k_{S}=4)$ and $C2\
(k_{N}=2,k_{S}=4).$ In model $C4$, the coefficients $g_{1}=g_{5}=0$ and
$\lambda_{2}=0.$ On the other hand, in model $C2,$ all the coefficients in
eq.(16) are non-zero.

  4. 4.

For $k_{L}=3,$ two models $D1\ (k_{N}=1,k_{S}=5)$ and $D2\ (k_{N}=3,k_{S}=1)$
are possible. Only the coefficient $g_{5}=0$ in eq.(16) while all the terms
are allowed for model $D2.$

  5. 5.

For $k_{L}=4,$ two models $E1\ (k_{N}=0,k_{S}=4)$ and $E2\ (k_{N}=2,k_{S}=4)$
are possible. All the coefficients in eq.(16) are allowed for these two
models.

  6. 6.

Finally, for $k_{L}=5,$ only one model is possible which is $F\
(k_{N}=1,k_{S}=5)$. All the terms in the superpotential $W_{\nu}$ are allowed
in this model.

We have summarized these models in Table 3.

In the present work, we shall discuss the numerical results and consistency
concerning neutrino oscillation data for only two models, A2 and B1. It is
important to note that in the models A1, C1, and E1, the weight of heavy
right-handed neutrino $k_{N}=0$. This condition gives $\lambda_{2}=0$ in
eq.(16). The resulting heavy Majorana mass matrices $M_{R}$ in each model are
singular; hence, these models are impossible.

Table 2: Different models for different choices of $k_{L}$ keeping the maximum weight of modular forms at 6. $k_{L}$ | $k_{N}$ | $k_{S}$ | Models | Coefficients($g_{i}$)
---|---|---|---|---
0 | 0 | 4 | A1 | $g_{2}=g_{3}=g_{4}=0$
2 | 4 | A2 | $g_{1}=g_{5}=0$
1 | 1 | 5 | B1 | $g_{1}=g_{5}=0$
3 | 1 | B2 | $g_{5}=0$
2 | 0 | 4 | C1 | $g_{1}=g_{5}=0$
2 | 4 | C2 | Nil
3 | 1 | 5 | D1 | $g_{5}=0$
3 | 1 | D2 | Nil
4 | 0 | 4 | E1 | Nil
2 | 4 | E2 | Nil
5 | 1 | 5 | F | Nil

Table 3: Ranges of input parameters used in numerical analysis. Free Parameters | Range
---|---
$g_{o2}$ | 1 - 10
$g_{o3}$ | 1 - 3
$|\lambda|$ | 0.01 - 10
$\phi_{\lambda}$ | [$-\pi,\pi]$
$\lambda_{2o}$ | 0.01 - 10
$|\lambda_{3o}|$ | 0.01 - 10
$\phi_{\lambda_{3o}}$ | [$-\pi,\pi]$
$m_{s}$ | (0.8 - 50) eV

###  II.1 Model A2 ($k_{L}=0,k_{N}=2,k_{S}=4$)

When we take $k_{L}=0,$ the possible choices of $k_{N}$ are $0,2,4,6.$
However, $k_{N}=\\{4,6\\}$ will be possible222For instance, if $k_{N}=6,$ the
Majorana mass terms $\lambda Y^{(k)}N_{i}^{c}N_{i}$ requires weight $k=12$.
only when the weight of the Yukawa, $k>6$. For model A2, we consider
$k_{L}=0,k_{N}=2$ and $k_{S}=4$. The coefficients $g_{1}=g_{5}=0$ and eq.(16)
becomes

$\displaystyle W_{\nu}$ $\displaystyle=$ $\displaystyle
g_{2}Y_{2}^{(2)}H_{d}\bar{L_{i}}N_{3}+g_{3}Y_{2}^{(2)}H_{d}\bar{L_{3}}N_{i}+g_{4}Y_{2}^{(2)}H_{d}\bar{L_{i}}N_{i}$
(17) $\displaystyle+$
$\displaystyle\lambda_{1}\left(Y_{1}^{(4)}\bar{N}^{c}_{i}N_{i}+Y_{2}^{(4)}\bar{N_{i}}^{c}N_{i}\right)+\lambda_{2}\left(Y_{2}^{(4)}\bar{N_{i}}^{c}N_{3}+Y_{2}^{(4)}\bar{N_{3}}^{c}N_{i}\right)$
$\displaystyle+$
$\displaystyle\lambda_{3}Y_{1}^{(4)}\bar{N_{3}}^{c}N_{3}+\lambda_{s1}Y_{1}^{(6)}\bar{S}^{c}N_{3}+\lambda_{s2}Y_{2}^{(6)}\bar{S}^{c}N_{i}.$

After EWSB, the resulting Dirac, Majorana, and sterile neutrino mass matrices
are

$\displaystyle M_{D}$ $\displaystyle=$
$\displaystyle\frac{vg_{4}}{\sqrt{2}}\left(\begin{array}[]{ccc}y_{2}&y_{1}&g_{o2}y_{1}\\\
y_{1}&-y_{2}&g_{o2}\ y_{2}\\\ -g_{o3}\ y_{2}&g_{o3}\ y_{1}&0\\\
\end{array}\right),\
M_{S}=\lambda_{s1}\left(\begin{array}[]{ccc}y_{1}^{3}+y_{1}y_{2}^{2}&y_{2}^{3}+y_{1}^{2}y_{2}&\lambda(3y_{1}^{2}y_{2}-y_{2}^{3})\end{array}\right)$
(22) $\displaystyle M_{R}$ $\displaystyle=$
$\displaystyle\lambda_{1}\left(\begin{array}[]{ccc}2y_{1}^{2}&2y_{2}y_{1}&\lambda_{2o}2y_{1}y_{2}\\\
2y_{2}y_{1}&2y_{2}^{2}&\lambda_{2o}(y_{1}^{2}-y_{2}^{2})\\\
\lambda_{2o}2y_{1}y_{2}&\lambda_{2o}(y_{1}^{2}-y_{2}^{2})&\lambda_{3o}(y_{1}^{2}+y_{2}^{2})\\\
\end{array}\right).$ (26)

where the coefficients
$g_{o2}=g_{2}/g_{4},g_{o3}=g_{3}/g_{4},\lambda=\lambda_{s1}/\lambda_{s2},\lambda_{2o}=\lambda_{2}/\lambda_{1}$
and $\lambda_{3o}=\lambda_{3}/\lambda_{1}$.

The active and sterile neutrino masses in the MES mechanism are given by

$m_{\nu}\simeq
M_{D}M_{R}^{-1}M_{S}^{T}\left(M_{S}M_{R}^{-1}M_{S}^{T}\right)^{-1}M_{S}\left(M_{R}^{-1}\right)^{T}M_{D}^{T}-M_{D}M_{R}^{-1}M_{D}^{T};$
(27) $m_{s}\simeq-M_{S}M_{R}^{-1}M_{S}^{T}.$ (28)

Using these,the resulting sterile neutrino mass is given by

$m_{s}=-\frac{\lambda_{s1}^{2}}{2\lambda_{1}}(y_{1}^{2}+y_{2}^{2})^{2}.$ (29)

###  II.2 Model B1 $(k_{L}=1,k_{N}=1,k_{S}=5)$

In this case, the weights of the fermion fields will be $k_{L}=1,k_{N}=1$, and
$k_{S}=5.$ These assignments of the modular weights imply that the
coefficients $g_{1}=g_{5}=0$ and $\lambda_{3}=0$ since all the singlet modular
forms of weight 2 will vanish. The modular invariant terms of $W_{\nu}$ in
eq.(16) are

$\displaystyle W_{\nu}$ $\displaystyle=$ $\displaystyle\
g_{2}Y_{2}^{(2)}H_{d}\bar{L_{i}}N_{3}+g_{3}Y_{2}^{(2)}H_{d}\bar{L_{3}}N_{i}+g_{4}Y_{2}^{(2)}H_{d}\bar{L_{i}}N_{i}$
(30) $\displaystyle+$
$\displaystyle\lambda_{1}Y_{2}^{(2)}\bar{N_{i}}^{c}N_{i}+\lambda_{2}(Y_{2}^{(2)}\bar{N_{i}}^{c}N_{3}+Y_{2}^{(2)}\bar{N_{3}}^{c}N_{i})+\lambda_{s1}Y_{1}^{(6)}\bar{S}^{c}N_{3}$
$\displaystyle+$ $\displaystyle\lambda_{s2}Y_{2}^{(6)}\bar{S}^{c}N_{i}.$

After EWSB, the Dirac, Majorana, and sterile neutrino mass matrices are
obtained as

$\displaystyle M_{D}$ $\displaystyle=$
$\displaystyle\frac{vg_{4}}{\sqrt{2}}\left(\begin{array}[]{ccc}y_{2}&y_{1}&g_{o2}y_{1}\\\
y_{1}&-y_{2}&g_{o2}\ y_{2}\\\ -g_{o3}\ y_{2}&g_{o3}\ y_{1}&0\\\
\end{array}\right),\
M_{S}=\lambda_{s1}\left(\begin{array}[]{ccc}y_{1}^{3}+y_{1}y_{2}^{2}&y_{2}^{3}+y_{1}^{2}y_{2}&\lambda(3y_{1}^{2}y_{2}-y_{2}^{3})\end{array}\right)$
(35) $\displaystyle M_{R}$ $\displaystyle=$
$\displaystyle\lambda_{1}\left(\begin{array}[]{ccc}y_{2}&y_{1}&\lambda_{2o}y_{1}\\\
y_{1}&-y_{2}&\lambda_{2o}y_{2}\\\ \lambda_{2o}y_{1}&\lambda_{2o}y_{2}&0\\\
\end{array}\right).$ (39)

Here, $M_{D}$ and $M_{S}$ in this case are same as model A1 but $M_{R}$ is
different. The resulting sterile neutrino mass matrix becomes

$m_{s}=\frac{\lambda_{s1}^{2}}{\lambda_{1}}\frac{\lambda\
y_{2}}{\lambda_{2o}^{2}}(y_{2}^{2}-3y_{1}^{2})(y_{1}^{2}+y_{2}^{2})(2\lambda_{2o}-\lambda).$
(40)

##  III Numerical analysis

The steps of numerical analysis are conducted as follows. The modulus
parameter $\tau$ is scanned in the fundamental domain as given below

$Re[\tau]=[-0.5,0.5],\ Im[\tau]=[0.6,2].$ (41)

For the charged lepton sector, the unknown coefficients $\alpha,\beta,\gamma$
are solved using the identities

$\displaystyle Tr[M_{L}]$ $\displaystyle=m_{e}^{2}+m_{\mu}^{2}+m_{\tau}^{2},$
(42) $\displaystyle Det[M_{L}]$ $\displaystyle=m_{e}^{2}\times
m_{\mu}^{2}\times m_{\tau}^{2},$ (43) $\displaystyle
Tr[M_{L}]^{2}/2-Tr[M_{L}^{2}]/2$
$\displaystyle=m_{e}m_{\mu}+m_{\mu}m_{\tau}+m_{\tau}m_{e}.$ (44)

where $m_{e},m_{\mu}$ and $m_{\tau}$ are the charged lepton masses and
$M_{L}=m_{L}m_{L}^{\dagger}$ is a hermitian matrix. The charged lepton mass
matrix can be diagonalized using the relation
$M_{L}^{diag}=U_{L}M_{L}U_{L}^{\dagger}.$

In the neutrino sector, the unknown free parameters of the model are randomly
scanned in the ranges given in Table 3 along with the values of $\tau$, which
reproduce the charged lepton masses. Now, the active neutrino mass matrix can
be numerically diagonalised using the relation
$M_{\nu}^{diag}=U_{\nu}M_{\nu}U_{\nu}^{\dagger},$ where
$M_{\nu}=m_{\nu}m_{\nu}^{\dagger}$ and
$M_{\nu}^{diag}=diag(m_{1}^{2},m_{2}^{2},m_{3}^{2}).$ Then, the $3\times 3$
active neutrino mixing matrix will be given by

$\mathcal{U}_{PMNS}=U_{L}^{\dagger}U_{\nu}$ (45)

Due to active-sterile mixing, the $U_{PMNS}$ is no longer unitary. The full
$(4\times 4)$ mixing matrix in the 3+1 scheme is given by

$U\simeq\left(\begin{array}[]{ccc}(1-\frac{1}{2}RR^{\dagger})\mathcal{U}_{PMNS}&R\\\
-R^{\dagger}\mathcal{U}_{PMNS}&1-\frac{1}{2}R^{\dagger}R\end{array}\right),$
(46)

where $R$ gives the strength of active-sterile mixing given by

$R=M_{D}M_{R}^{-1}M_{S}^{T}(M_{S}M_{R}^{-1}M_{S}^{T})^{-1}.$ (47)

The three active neutrino mixing angles are determined using the following
relations

$\displaystyle\sin^{2}\theta_{12}\ =\
\frac{|U_{e2}|^{2}}{1-|U_{e4}|^{2}-|U_{e3}|^{2}},$ $\displaystyle\ \
\sin^{2}\theta_{13}\ =\ \frac{|U_{e3}|^{2}}{1-|U_{e4}|^{2}},$ (48)
$\displaystyle\sin^{2}\theta_{23}\ =\
\frac{|U_{e3}|^{2}(1-|U_{e4}|^{2})-|U_{e4}|^{2}|U_{\mu
4}|^{2}}{1-|U_{e4}|^{2}-|U_{\mu 4}|^{2}}$ $\displaystyle+\frac{|U_{e1}U_{\mu
1}+U_{e2}U_{\mu
2}|^{2}(1-|U_{e4}|^{2})}{1-|U_{e4}|^{2}-|U_{e3}|^{2})(1-|U_{e4}|^{2}-|U_{\mu
4}|^{2}}.$

The free parameters of the model are filtered by constraining the three mixing
angles and the ratio of two mass squared differences $r=\sqrt{\frac{\Delta
m_{21}^{2}}{|\Delta m_{3l}^{2}|}}=m_{2}/m_{3}$ using the latest neutrino
oscillation data from Nufit 5.20. The absolute scale of the active neutrino
mass is fixed by the overall factor $g_{4}^{2}v^{2}/2\lambda_{1}$. We have
also imposed the Planck cosmological upper limit on the sum of active neutrino
masses, $\sum m_{i}<0.12$ eV.

Figure 1: For model A2 : (a). Plot between coefficients $\alpha,\beta$ and
$\gamma$ . (b). Variation between the Yukawa components $|y_{1}|$ and
$|y_{2}|$.

Figure 2: (a). Plot between Jarlskog invariant $J$ with active neutrino mixing
angles $\sin^{2}\theta_{23}$ and $\sin^{2}\theta_{12}$ predicted from the
model A2. (b). Variation between $m_{2}$, $m_{3}$ ($m_{1}=0$) with $\sum
m_{i}$.

Figure 3: Results for model B1 : (a). Plot between coefficients $\alpha,\beta$
and $\gamma$ . (b). Variation between the Yukawa components $|y_{1}|$ and
$|y_{2}|$.

Figure 4: (a). Plot between Jarlskog invariant $J$ with active neutrino mixing
angles $\sin^{2}\theta_{23}$ and $\sin^{2}\theta_{12}$ predicted from the
model B1. (b). Variation between $m_{2}$, $m_{3}$ ($m_{1}=0$) with $\sum
m_{i}$.

##  IV Results

In this section, we discuss the results and predictions of neutrino
observables of the model for the normal hierarchy(NH) of neutrino masses only.
For model A2, the values of unknown parameters evaluated for the charged
lepton sector are shown as variation plots between them in Fig.(1)(a), whereas
the components of the Yukawa couplings are plotted in Fig.(1)(b). The
variation of Jarlskog invariant $J$ with the active neutrino mixing angles is
shown in Fig.(2)(a). The absolute scale of the active neutrino masses is fixed
by choosing the overall factor $g_{4}^{2}v^{2}/2\lambda_{1}=0.01$ eV.
Variation of active neutrino masses $m_{2},m_{3}$ (where, $m_{1}=0$ for NH in
MES mechanism) with sum of active neutrino mass $\sum m_{i}$ is shown in
Fig.(2)(b). Similarly, the numerical results of model B1 are shown in Fig.(3)
and Fig.(4). We have also calculated the best-fit values of the free
parameters and predictions of neutrino observables based on the $\chi^{2}$
function defined as

$\chi^{2}(x_{i})=\sum_{j}\left(\frac{y_{j}(x_{i})-y_{j}^{bf}}{\sigma_{j}}\right)^{2},$
(49)

where $x_{i}$ are the free parameters in the model and $j$ is summed over the
observables
$\\{\sin^{2}\theta_{12},\sin^{2}\theta_{13},\sin^{2}\theta_{23},r\\}$. Here,
$y_{j}(x_{i})$ denotes the model predictions for the observables, and
$y_{j}^{bf}$ are their best-fit values obtained from the global analysis.
$\sigma_{j}$ denotes the corresponding uncertainties obtained by symmetrizing
$1\sigma$ range of the neutrino observables taken from Nufit 5.2 data. The
best-fit values correspond to the condition when $\chi^{2}(x_{i})$ is minimum.
We have summarised these values for models A2 and B1 in Table 4.

Table 4: Best-fit values of the model parameters and the corresponding
predictions of neutrino observables for $\chi^{2}_{min}=0.74$ for model $A2$
and $\chi^{2}_{min}=6.89$ for model $B1$. Here, the two Majorana phases are
$\phi_{1}$ and $\phi_{2}$.

Observables | Model A2 | Model B1
---|---|---
$\sin^{2}\theta_{23}$ | 0.573 | 0.568
$\sin^{2}\theta_{13}$ | 0.021 | 0.022
$\sin^{2}\theta_{12}$ | 0.295 | 0.278
$r$ | 0.173 | 0.169
$\delta_{CP}/^{o}$ | 315.48 | 344.06
$\phi_{1}/^{o}$ | 339.55 | 351.91
$\phi_{2}/^{o}$ | 347.08 | 354.73
$m_{2}$ | 0.0035 eV | 0.00039 eV
$m_{3}$ | 0.0202 eV | 0.0023 eV
$\sum$ $m_{i}$ | 0.0237 eV | 0.00269 eV

Model parameters | Model A2 | Model B1
---|---|---
Re$[\tau]$ | 0.438 | 0.364
Im$[\tau]$ | 1.14 | 0.739
$\alpha$ | 0.0289 | 0.15 $\times 10^{-4}$
$\beta$ | 0.2 $\times 10^{-4}$ | 0.0102
$\gamma$ | 0.0102 | 0.0123
$g_{o2}$ | 5.996 | 3.32
$g_{o3}$ | 1.608 | 1.98
$\lambda$ | 0.156 | 1.34
$\lambda_{2o}$ | 3.441 | 7.89
$\lambda_{3o}$ | 1.232 | -

##  V Summary and Discussion

We have successfully developed a model-building strategy of leptons using
modular $S_{3}$ symmetry without using any hypothetical flavons in the 3+1
framework of active-sterile neutrino mixing. We obtain a non-diagonal charged
lepton mass matrix containing three unknown parameters. Based on assignments
of modular weights of SM lepton doublet, we have classified eleven different
models with the condition that the maximum allowed weight of the modular form
Yukawa coupling is 6. Out of these, only two models, A2 and B1, are discussed
in the present work through numerical analysis constrained by the latest
neutrino oscillation data and Planck upper limit on the sum of active neutrino
masses $\sum m_{i}<0.12$eV. In this work, only the normal hierarchy is
discussed. The best fit-values of the model predictions of both models suggest
$\sin^{2}\theta_{23}>0.5$, implying the higher octant of atmospheric mixing
angle.

In the numerical analysis, we performed a random scan of about $10^{6}$
iterations of the free parameters and fitted the results of neutrino
observables with experimental data. It is important to note that the number of
allowed data points obtained from both the models are different: significantly
less number of data points are obtained for model B1 compared to A2. Further
analysis and phenomenological implications of all the models in both the
normal and inverted hierarchies will be carried out in the future.

## Acknowledgements

One of the authors(MKS) would like to thank DST-INSPIRE, Govt. of India for
providing the fellowship for the research under INSPIRE fellowship(ID
IF180349).

## References

  * [1] He Zhang.  Light Sterile Neutrino in the Minimal Extended Seesaw.  Phys. Lett. B, 714:262–266, 2012.
  * [2] Ferruccio Feruglio.  Are neutrino masses modular forms?  In From My Vast Repertoire… Guido Altarelli’s Legacy, pages 227–266. World Scientific, 2019.
  * [3] Tatsuo Kobayashi, Naoya Omoto, Yusuke Shimizu, Kenta Takagi, Morimitsu Tanimoto, and Takuya H Tatsuishi.  Modular a4 invariance and neutrino mixing.  Journal of High Energy Physics, 2018(11), 2018.
  * [4] Davide Meloni and Matteo Parriciatu.  A simplest modular $s\\_3$ model for leptons.  arXiv preprint arXiv:2306.09028, 2023.
  * [5] Di Zhang.  A modular a4 symmetry realization of two-zero textures of the majorana neutrino mass matrix.  Nuclear Physics B, 952:114935, 2020.
  * [6] Hiroshi Okada and Yutaro Shoji.  A radiative seesaw model with three higgs doublets in modular a4 symmetry.  Nuclear Physics B, 961:115216, 2020.
  * [7] Gui-Jun Ding, Stephen F King, and Xiang-Gan Liu.  Neutrino mass and mixing with a 5 modular symmetry.  Physical Review D, 100(11):115005, 2019.
  * [8] Xin Wang, Bingrong Yu, and Shun Zhou.  Double covering of the modular a 5 group and lepton flavor mixing in the minimal seesaw model.  Physical Review D, 103(7):076005, 2021.
  * [9] JT Penedo and ST Petcov.  Lepton masses and mixing from modular s4 symmetry.  Nuclear Physics B, 939:292–307, 2019.
  * [10] Tatsuo Kobayashi, Yusuke Shimizu, Kenta Takagi, Morimitsu Tanimoto, and Takuya H. Tatsuishi.  ${A}_{4}$ lepton flavor model and modulus stabilization from ${S}_{4}$ modular symmetry.  Phys. Rev. D, 100:115045, Dec 2019.
  * [11] We have taken this choice to avoid Majorana mass terms of the sterile neutrinos of the form $\sim\Lambda YSS.$.
  * [12] For instance, if $k_{N}=6,$ the Majorana mass terms $\lambda Y^{(k)}N_{i}^{c}N_{i}$ requires weight $k=12$.


task [3.8em] 2.3em

# Deep Discriminative to Kernel Density Networks for Calibrated Inference

Jayanta Dey  1,∗ Will LeVine  1,2 Haoyin Xu  1,† Ashwin De Silva  1,† Tyler M.
Tomita  1 Ali Geisa  1 Tiffany Chu  1 Jacob Desman  1 and Joshua T.
Vogelstein1

###### Abstract

Deep discriminative approaches like random forests and deep neural networks
have recently found applications in many important real-world scenarios.
However, deploying these learning algorithms in safety-critical applications
raises concerns, particularly when it comes to ensuring confidence calibration
for both in-distribution and out-of-distribution data points. Many popular
methods for in-distribution (ID) calibration, such as isotonic regression and
Platt’s sigmoidal regression, exhibit excellent ID calibration performance but
often at the cost of classification accuracy. Moreover, these methods are not
calibrated for the entire feature space, leading to overconfidence in the case
of out-of-distribution (OOD) samples. In this paper, we leveraged the fact
that deep models, including both random forests and deep-nets, learn internal
representations which are unions of polytopes with affine activation functions
to conceptualize them both as partitioning rules of the feature space. We
replace the affine function in each polytope populated by the training data
with a Gaussian kernel. We propose sufficient conditions for our proposed
methods to be consistent estimators of the corresponding class conditional
densities. Moreover, our experiments on both tabular and vision benchmarks
show that the proposed approaches obtain well-calibrated posteriors while
mostly preserving or improving the classification accuracy of the original
algorithm for in-distribution region, and extrapolates beyond the training
data to handle out-of-distribution inputs appropriately.

††1Johns Hopkins University (JHU), 2Scale AI †† ∗ corresponding author:
<EMAIL_ADDRESS>
##  1 Introduction

Machine learning methods, specially deep neural networks and random forests
have shown excellent performance in many real-world tasks, including drug
discovery, autonomous driving and clinical surgery. However, calibrating
confidence over the whole feature space for these approaches remains a key
challenge in the field. Although these learning algorithms can achieve near
optimal performance at inferring on the samples lying in the high density
regions of the training data [1, 2, 3], they yield highly confident
predictions for the samples lying far away from the training data [4].
Calibrated confidence within the training or in-distribution (ID) region as
well as in the out-of-distribution (OOD) region is crucial for safety critical
applications like autonomous driving and computer-assisted surgery, where any
aberrant reading should be detected and taken care of immediately [4, 5]. A
well-calibrated model capable of quantifying the uncertainty associated with
inference for any points from the training distribution as well as detecting
OOD data can be a life-saver in these cases.

The approaches to calibrate OOD confidence for learning algorithms described
in the literature can be roughly divided into two groups: discriminative and
generative. Discriminative approaches try to scale the posteriors based on OOD
detection or modify the learning loss function. Intuitively, the easiest
solution for OOD confidence calibration is to learn a function that gives
higher scores for in-distribution samples and lower scores for OOD samples,
and thereby re-scale the posterior or confidence score from the original model
accordingly [6]. There are a number of approaches in the literature which try
to either modify the loss function [7, 8, 9] or adversarially train the
network to be less confident on OOD samples [10, 4]. Recently, as shown by
Hein et al. [4], the ReLU networks produce arbitrarily high confidence as the
inference point moves far away from the training data. Therefore, calibrating
ReLU networks for the whole OOD region is not possible without fundamentally
changing the network architecture. As a result, all of the aforementioned
algorithms are unable to provide any guarantee about the performance of the
network through out the whole feature space. On the other end of the spectrum,
the generative group tries to learn generative models for both the in-
distribution as well as the out-of-distribution samples. The general idea is
to conduct likelihood ratio tests for a particular sample out of the
generative models [11], or threshold the ID likelihoods to determine OOD
samples. However, it is not obvious how to control likelihoods far away from
the training data for powerful generative models like variational autoencoders
(VAEs) [12] and generative adversarial networks (GAN) [13]. Moreover,
Nalisnick et al. [14] and Hendrycks et al. [10] showed VAEs and GANs can also
yield overconfident likelihoods far away from the training data.

The algorithms described so far are concerned with OOD confidence calibration
for deep-nets only. However, in this paper, we show that other approaches
which partition the feature space, for example random forest, can also suffer
from poor confidence calibration both in the ID and the OOD regions. Moreover,
the algorithms described above are concerned about the confidence of the
algorithms in the OOD region only and they do not address the confidence
calibration within the training distribution at all. This issue is addressed
separately in a different group of literature [15, 16, 17, 18, 19, 20]. In
this paper, we consider both calibration problem jointly and propose an
approach that achieves good calibration throughout the whole feature space.

In this paper, we conceptualize both random forest and ReLU networks as
partitioning rules with an affine activation over each polytope. We consider
replacing the affine functions learned over the polytopes with Gaussian
kernels. We propose two novel kernel density estimation techniques named
Kernel Density Forest (KDF) and Kernel Density Network (KDN). We theoretically
show that they asymptotically converge to the true training distribution under
certain conditions. At the same time, the estimated likelihood from the kernel
density models decreases for samples far away from the training samples. By
adding a suitable bias to the kernel density estimate, we can achieve
calibrated posterior over the classes in the OOD region. It completely
excludes the need for providing OOD training examples to the model. We conduct
several simulation and real data studies that show both KDF and KDN are robust
against OOD samples while they maintain good performance in the in-
distribution region.

##  2 Related Works and Our Contributions

There are a number of approaches in the literature which attempt to learn a
generative model and control the likelihoods far away from the training data.
For example, Ren et al. [11] employed likelihood ratio test for detecting OOD
samples. Wan et al. [8] modified the training loss so that the downstream
projected features follow a Gaussian distribution. However, there is no
guarantee of performance for OOD detection for the above methods. To the best
of our knowledge, only Meinke et al. [5] has proposed an approach to guarantee
asymptotic performance for OOD detection. They model the training and the OOD
distribution using Gaussian mixture models which enable them to control the
class conditional posteriors far away. Compared to the aforementioned methods,
our approach differs in several ways:

  * •

We address the confidence calibration problems for both ReLU-nets and random
forests from a common ground.

  * •

We address in-distribution (ID) and out-of-distribution (OOD) calibration
problem as a continuum rather than two separate problems.

  * •

We provide guarantees for asymptotic convergence of our proposed approach
under certain conditions for both ID and OOD regions.

  * •

We propose an unsupervised OOD calibration approach, i.e., we do not need to
train exhaustively on different OOD samples.

##  3 Methods

###  3.1 Setting

Consider a supervised learning problem with independent and identically
distributed training samples $\\{(\mathbf{x}_{i},y_{i})\\}_{i=1}^{n}$ such
that $(X,Y)\sim P_{X,Y}$, where $X\sim P_{X}$ is a
$\mathcal{X}\subseteq\mathbb{R}^{d}$ valued input and $Y\sim P_{Y}$ is a
$\mathcal{Y}=\\{1,\cdots,K\\}$ valued class label. We define in-distribution
region as the high density region of $P_{X,Y}$ and denote it by $\mathcal{S}$.
Here the goal is to learn a confidence score,
$\mathbf{g}:\mathbb{R}^{d}\rightarrow[0,1]^{K}$,
$\mathbf{g}(\mathbf{x})=[g_{1}(\mathbf{x}),g_{2}(\mathbf{x}),\dots,g_{K}(\mathbf{x})]$
such that,

(1)
$g_{y}(\mathbf{x})=\begin{cases}P_{Y|X}(y|\mathbf{x}),&\leavevmode\nobreak\
\text{if}\leavevmode\nobreak\ \mathbf{x}\in\mathcal{S}\\\
P_{Y}(y),&\leavevmode\nobreak\ \text{if}\leavevmode\nobreak\
\mathbf{x}\notin\mathcal{S}\end{cases},\quad\forall y\in\mathcal{Y}$

where $P_{Y|X}(y|\mathbf{x})$ is the posterior probability for class $y$ given
by the Bayes formula:

(2)
$P_{Y|X}(y|\mathbf{x})=\frac{P_{X|Y}(\mathbf{x}|y)P_{Y}(y)}{\sum_{k=1}^{K}P_{X|Y}(\mathbf{x}|k)P_{Y}(k)},\quad\forall
y\in\mathcal{Y}.$

Here $P_{X|Y}(\mathbf{x}|y)$ is the class conditional density for the training
data which we will refer as $f_{y}(\mathbf{x})$ hereafter for brevity.

###  3.2 Background and Main Idea

Deep discriminative networks partition the feature space $\mathbb{R}^{d}$ into
a union of $p$ affine polytopes $Q_{r}$ such that
$\bigcup_{r=1}^{p}Q_{r}=\mathbb{R}^{d}$, and learn an affine function over
each polytope [4, 21]. Mathematically, the class-conditional density for the
label $y$ estimated by these deep discriminative models at a particular point
$\mathbf{x}$ can be expressed as:

(3)
$\hat{f}_{y}(\mathbf{x})=\sum_{r=1}^{p}(\mathbf{a}_{r}^{\top}\mathbf{x}+b_{r})\mathbbm{1}(\mathbf{x}\in
Q_{r}).$

For example, in the case of a decision tree, $\mathbf{a}_{r}=\mathbf{0}$,
i.e., decision tree assumes uniform distribution for the class-conditional
densities over the leaf nodes. Among these polytopes, the ones that lie on the
boundary of the training data extend to the whole feature space and hence
encompass all the OOD samples. Since the posterior probability for a class is
determined by the affine activation over each of these polytopes, the
algorithms tend to be overconfident when making predictions on the OOD inputs.
Moreover, there exist some polytopes that are not populated with training
data. These unpopulated polytopes serve to interpolate between the training
sample points. If we replace the affine activation function of the populated
polytopes with Gaussian kernel $\mathcal{G}$ learned using maximum likelihood
approach on the training points within the corresponding polytope and prune
the unpopulated ones, the tail of the kernel will help interpolate between the
training sample points while assigning lower likelihood to the low density or
unpopulated polytope regions of the feature space. This may result in better
confidence calibration for the proposed modified approach.

###  3.3 Proposed Model

Consider the collection of polytope indices $\mathcal{P}$ which contains the
indices of total $\tilde{p}$ polytopes populated by the training data. We
replace the affine functions over the populated polytopes with Gaussian
kernels $\mathcal{G}(\cdot;\hat{\mu}_{r},\hat{\Sigma}_{r})$. For a particular
inference point $\mathbf{x}$, we consider the Gaussian kernel with the minimum
distance from the center of the kernel to the corresponding point:

(4)
$r^{*}_{\mathbf{x}}=\operatornamewithlimits{argmin}_{r}\|\mu_{r}-\mathbf{x}\|,$

where $\|\cdot\|$ denotes a distance. As we will show later, the type of
distance metric considered in Equation 4 highly impacts the performance of the
proposed model. In short, we modify Equation 3 from the parent ReLU-net or
random forest to estimate the class-conditional density as:

(5)
$\tilde{f}_{y}(\mathbf{x})=\frac{1}{n_{y}}\sum_{r\in\mathcal{P}}n_{ry}\mathcal{G}(\mathbf{x};\mu_{r},\Sigma_{r})\mathbbm{1}(r=r^{*}_{\mathbf{x}}),$

where $n_{y}$ is the total number of samples with label $y$ and $n_{ry}$ is
the number of samples from class $y$ that end up in polytope $Q_{r}$. We add a
bias to the class conditional density $\tilde{f}_{y}$:

(6) $\hat{f}_{y}(\mathbf{x})=\tilde{f}_{y}(\mathbf{x})+\frac{b}{\log(n)}.$

Note that in Equation 6, $\frac{b}{\log(n)}\rightarrow 0$ as the total
training points, $n\rightarrow\infty$. The confidence score
$\hat{g}_{y}(\mathbf{x})$ for class $y$ given a test point $\mathbf{x}$ is
estimated using the Bayes rule as:

(7)
$\hat{g}_{y}(\mathbf{x})=\frac{\hat{f}_{y}(\mathbf{x})\hat{P}_{Y}(y)}{\sum_{k=1}^{K}\hat{f}_{k}(\mathbf{x})\hat{P}_{Y}(k)},$

where $\hat{P}_{Y}(y)$ is the empirical prior probability of class $y$
estimated from the training data. We estimate the class for a particular
inference point $\mathbf{x}$ as:

(8)
$\hat{y}=\operatornamewithlimits{argmax}_{y\in\mathcal{Y}}\hat{g}_{y}(\mathbf{x}).$

###  3.4 Desiderata

We desire our proposed model to estimate confidence score $\hat{g}_{y}$ to
satisfy the following two desiderata:

  1. 1.

Asymptotic Performance: We want point-wise convergence for our estimated
confidence as $n\to\infty$, i.e.,

$\max_{y\in\mathcal{Y}}\sup_{\mathbf{x}\in\mathbb{R}^{d}}|g_{y}(\mathbf{x})-\hat{g}_{y}(\mathbf{x})|\to
0.$

  2. 2.

Finite Sample Performance: We want better posterior calibration for
$\hat{g}_{y}(\mathbf{x})$ both in ID and OOD regions compared to that of its
parent model.

We theoretically derive the conditions under which we achieve Desiderata 1 in
Section 4. However, we run extensive experiments on various simulation and
benchmark datasets in Section 6 to empirically verify that our proposed
approach achieves Desiderata 2.

##  4 Theoretical Results

######  Theorem 4.1 (Asymptotic Convergence to the True Distribution).

Consider a partition rule that partitions $\mathbb{R}^{d}$ into an infinite
number of hypercubes of the same size $h_{n}>0$. Formally, let
$\mathcal{P}_{n}=\\{Q_{1},Q_{2},\cdots\\}$ be a partition of $\mathbb{R}^{d}$,
that is, it partitions $\mathbb{R}^{d}$ into sets of the type
$\Pi_{i=1}^{d}[\psi_{i}h_{n},(\psi_{i}+1)h_{n})$, where $\psi_{i}$’s are
integers. Let $n$ be the total number of samples and $n_{r}$ be the number of
data points within the polytope $Q_{r}$. Consider the true probability density
$f$ which is estimated using the samples populating the polytopes using
Equation 5, denoted as $\hat{f}$. The conditions for choosing the Gaussian
kernel parameters are:

  1. 1.

The center of the kernel can be any point $z_{r}$ within the polytope $Q_{r}$
as $n\to\infty$,

  2. 2.

The kernel bandwidth along any dimension $\sigma_{r}$ is any positive number
always bounded by the polytope bandwidth $h_{n}$ as $n\to\infty$, i.e.,
$\sigma_{r}=C_{r}h_{n}$, where $0<C_{r}<1$.

Consider the following assumptions as well:

  1. 1.

The polytope bandwidth $h_{n}\to 0$ as $n\to\infty$.

  2. 2.

$n$ grows faster than the shrinkage of $h_{n}$, i.e., $nh_{n}\to\infty$ in
probability as $h_{n}\to 0$.

Given these assumptions, we have that as $n\to\infty$:

(9)
$\displaystyle\sup_{\mathbf{x}\in\mathbb{R}^{d}}|f(\mathbf{x})-\hat{f}(\mathbf{x})|\to
0,$

where $|\cdot|$ denotes absolute value of the scalar it operates on.

######  Proof 4.2.

Please see Appendix A for the proof.

######  Theorem 4.3 (Asymptotic OOD Convergence).

Given $n$ independent and identically distributed training samples
$\\{(\mathbf{x}_{i},y_{i})\\}_{i=1}^{n}$, we define the distance of an
inference point $\mathbf{x}$ from the training points as:
$d_{\mathbf{x}}=\min_{i=1,\cdots,n}\|\mathbf{x}-\mathbf{x}_{i}\|$. Here
$\|\cdot\|$ denotes a distance as mentioned in Equation 4. Given non-zero and
bounded bandwidth of the Gaussians, then we have almost sure convergence for
$\hat{g}_{y}$ as:

$\lim_{d_{\mathbf{x}}\to\infty}\hat{g}_{y}(\mathbf{x})\overset{as}{\to}\hat{P_{Y}}(y).$

######  Proof 4.4.

Please see Appendix A for the proof.

######  Corollary 4.5.

Given the conditions in Theorem 4.1 and 4.3, we have:

$\displaystyle\max_{y\in\mathcal{Y}}\sup_{\mathbf{x}\in\mathbb{R}^{d}}|g_{y}(\mathbf{x})-\hat{g}_{y}(\mathbf{x})|\to
0.$

######  Proof 4.6.

Using the law of large numbers, we have
$\hat{P}_{Y}(y)=\frac{n_{y}}{n}\overset{as}{\to}P_{Y}(y)$ as $n\to\infty$. The
rest of the proof follows from Theorem 4.1 and 4.3.

######  Corollary 4.7.

There exists a class of random forest algorithms which partitions the feature
space satisfying the conditions in Theorem 4.1.

######  Proof 4.8.

Please see Appendix A for the proof.

######  Corollary 4.9.

There exists a class of deep-net algorithms which partitions the feature space
satisfying the conditions in Theorem 4.1.

######  Proof 4.10.

Please see Appendix A for the proof.

##  5 Model Parameter Estimation

###  5.1 Gaussian Kernel Parameter Estimation

Theorem 4.1 implies that the Gaussian kernel parameters for each polytope need
to maintain two key properties for ensuring asymptotic optimal performance of
our proposed approach. We fit Gaussian kernel parameters to the samples that
end up in the $r$-th polytope using maximum likelihood approach. To satisfy
the first condition, we set the kernel center as:

(10)
$\hat{\mu}_{r}=\frac{1}{n_{r}}\sum_{i=1}^{n}\mathbf{x}_{i}\mathbbm{1}(\mathbf{x}_{i}\in
Q_{r}).$

To satisfy the second condition, we set the kernel covariance matrix as:

(11) $\hat{\Sigma}_{r}=\frac{\sum_{i=1}^{n}\mathbbm{1}(\mathbf{x}_{i}\in
Q_{r})(\mathbf{x}_{i}-\hat{\mu}_{r})(\mathbf{x}_{i}-\hat{\mu}_{r})^{\top}+\lambda
I_{d}}{\sum_{i=1}^{n}\mathbbm{1}(\mathbf{x}_{i}\in Q_{r})},$

where we regularize the Frobenius norm of precision matrix $\Sigma^{-1}$ so
that $\Sigma$ does not become singular and $\lambda$ is the regularization
parameter. Note that $\hat{\mu}_{r}$ in Equation 10 resides always within the
corresponding polytope $Q_{r}$ and the estimated Gaussian bandwidth will be
bounded by the polytope bandwidth. We constrain our estimated Gaussian kernels
to have diagonal covariance matrix and set the off-diagonal elements of
$\hat{\Sigma}_{r}$ to $0$.

###  5.2 Sample Size Ratio Estimation

For a high dimensional dataset with low training sample size, the polytopes
are sparsely populated with training samples. For improving the estimate of
the ratio $\frac{n_{ry}}{n_{y}}$ in Equation 5, we incorporate the samples
from other polytopes $Q_{s}$ based on the similarity $w_{rs}$ between $Q_{r}$
and $Q_{s}$ as:

$\displaystyle\hat{\left(\frac{n_{ry}}{n_{y}}\right)}$
$\displaystyle=\frac{\tilde{w}_{ry}}{\tilde{w}_{y}}=\frac{\tilde{w}_{ry}}{\sum_{r\in\mathcal{P}}\tilde{w}_{ry}}$
(12)
$\displaystyle=\frac{\sum_{s\in\mathcal{P}}\sum_{i=1}^{n}w_{rs}\mathbbm{1}(\mathbf{x}_{i}\in
Q_{s})\mathbbm{1}(y_{i}=y)}{\sum_{r\in\mathcal{P}}\sum_{s\in\mathcal{P}}\sum_{i=1}^{n}w_{rs}\mathbbm{1}(\mathbf{x}_{i}\in
Q_{s})\mathbbm{1}(y_{i}=y)}.$

As $n\to\infty$ and $h_{n}\to 0$, the estimated weights $w_{rs}$ should
satisfy the condition:

(13) $w_{rs}\to\begin{cases}0,&\text{if }Q_{r}\neq Q_{s}\\\ 1,&\text{if
}Q_{r}=Q_{s}.\end{cases}$

For simplicity, we will describe the estimation procedure for $w_{rs}$ in the
next sections. Note that if we satisfy Condition 13, then we have
$\frac{\tilde{w}_{ry}}{\tilde{w}_{y}}\to\frac{n_{ry}}{n_{y}}$ as $n\to\infty$.
Therefore, we modify Equation 5 as:

(14)
$\tilde{f}_{y}(\mathbf{x})=\frac{1}{\tilde{w}_{y}}\sum_{r\in\mathcal{P}}\tilde{w}_{ry}\mathcal{G}(\mathbf{x};\hat{\mu}_{r},\hat{\Sigma}_{r})\mathbbm{1}(r=\hat{r}^{*}_{\mathbf{x}}),$

where
$\hat{r}^{*}_{\mathbf{x}}=\operatornamewithlimits{argmin}_{r}\|\hat{\mu}_{r}-\mathbf{x}\|$.
Now we use $\tilde{f}_{y}(\mathbf{x})$ estimated using (5) in Equation (6),
(7) and (8), respectively. Below, we describe how we estimate $w_{rs}$ for KDF
and KDN .

Figure 1: Simulation datasets, Classification error, Hellinger distance from
true posteriors, mean max confidence or posterior for the simulation
experiments. The simulation experiments were repeated $45$ times. The median
performance is shown as a dark curve with shaded region as error bars showing
the $25$-th and the $75$-th percentile. KDF (Left block) and KDN (Right block)
improve both in- and out-of-distribution calibration of their respective
parent algorithms while maintaining nearly similar classification accuracy on
the simulation datasets. For Trunk simulation sample size was fixed at $5000$
for gradually higher number of dimensions and OOD test samples were sampled
uniformly from a hyper-sphere with radius fixed at $20$. Trunk simulation
dataset is visualized for the first two dimensions.

###  5.3 Kernel Density Forest

Consider $T$ number of decision trees in a random forest trained on $n$ i.i.d
training samples $\\{(\mathbf{x}_{i},y_{i})\\}_{i=1}^{n}$. Each tree $t$
partitions the feature space into $p_{t}$ polytopes resulting in a set of
polytopes: $\\{\\{Q_{t,r}\\}_{r=1}^{p_{t}}\\}_{t=1}^{T}$. The intersection of
these polytopes gives a new set of polytopes $\\{Q_{r}\\}_{r=1}^{p}$ for the
forest. For any point $\mathbf{x}_{r}\in Q_{r}$, we push every other sample
point $\mathbf{x}_{s}\in Q_{s}$ down the trees. Now, we define the kernel
$\mathcal{K}(r,s)$ as:

(15) $\mathcal{K}(r,s)=\frac{t_{rs}}{T},$

where $t_{rs}$ is the total number trees $\mathbf{x}_{r}$ and $\mathbf{x}_{s}$
end up in the same leaf node. Note that $0\leq\mathcal{K}(r,s)\leq 1$. If the
two samples end up in the same leaves in all the trees, they belong to the
same polytope, i.e. $Q_{r}=Q_{s}$.

In short, $\mathcal{K}(r,s)$ is the fraction of total trees where the two
samples follow the same path from the root to a leaf node. We exponentiate
$\mathcal{K}(r,s)$ so that Condition 13 is satisfied:

(16) $w_{rs}=\mathcal{K}(r,s)^{k\log n}.$

We choose $k$ using grid search on a hold-out dataset.

###  5.4 Kernel Density Network

Consider a fully connected ReLU-net trained on $n$ $iid$ training samples
$\\{(\mathbf{x}_{i},y_{i})\\}_{i=1}^{n}$. We have the set of all nodes denoted
by $\mathcal{N}_{l}$ at a particular layer $l$. We can randomly pick a node
$n_{l}\in\mathcal{N}_{l}$ at each layer $l$, and construct a sequence of nodes
starting at the input layer and ending at the output layer which we call an
activation path: $m=\\{n_{l}\in\mathcal{N}_{l}\\}_{l=1}^{L}$. Note that there
are $N=\Pi_{i=1}^{L}{|\mathcal{N}_{l}|}$ possible activation paths for a
sample in the ReLU-net. We index each path by a unique identifier number
$z\in\mathbb{N}$ and construct a sequence of activation paths as:
$\mathcal{M}=\\{m_{z}\\}_{z=1,\cdots,N}$. Therefore, $\mathcal{M}$ contains
all possible activation pathways from the input to the output of the network.

While pushing a training sample $\mathbf{x}_{i}$ through the network, we
define the activation from a ReLU unit at any node as ‘$1$’ when it has non-
negative input and ‘$0$’ otherwise. Therefore, the activation indicates on
which side of the affine function at each node the sample falls. The
activation for all nodes in an activation path $m_{z}$ for a particular sample
creates an activation mode $a_{z}\in\\{0,1\\}^{L}$. If we evaluate the
activation mode for all activation paths in $\mathcal{M}$ while pushing a
sample through the network, we get a sequence of activation modes:
$\mathcal{A}_{r}=\\{a_{z}^{r}\\}_{z=1}^{N}$. Here $r$ is the index of the
polytope where the sample falls in.

If the two sequences of activation modes for two different training samples
are identical, they belong to the same polytope. In other words, if
$\mathcal{A}_{r}=\mathcal{A}_{s}$, then $Q_{r}=Q_{s}$. This statement holds
because the above samples will lie on the same side of the affine function at
each node in different layers of the network. Now, we define the kernel
$\mathcal{K}(r,s)$ as:

(17)
$\mathcal{K}(r,s)=\frac{\sum_{z=1}^{N}\mathbbm{1}(a_{z}^{r}=a_{z}^{s})}{N}.$

Note that $0\leq\mathcal{K}(r,s)\leq 1$. In short, $\mathcal{K}(r,s)$ is the
fraction of total activation paths which are identically activated for two
samples in two different polytopes $r$ and $s$. We exponentiate the weights
using Equation 16.

Pseudocodes outlining the two algorithms are provided in Appendix C.

###  5.5 Geodesic Distance

Consider $\mathcal{P}_{n}=\\{Q_{1},Q_{2},\cdots,Q_{p}\\}$ as a partition of
$\mathbb{R}^{d}$ given by a random forest or a ReLU-net after being trained on
$n$ training samples. We measure distance between two points $\mathbf{x}\in
Q_{r},\mathbf{x}^{\prime}\in Q_{s}$ using the kernel introduced in Equation 15
and Equation 17 [22] and call it ‘Geodesic’ distance:

(18)
$d(r,s)=-\mathcal{K}(r,s)+\frac{1}{2}(\mathcal{K}(r,r)+\mathcal{K}(s,s))=1-\mathcal{K}(r,s)$

######  Corollary 5.1.

$(\mathcal{P}_{n},d)$ is a metric space.

######  Proof 5.2.

See Appendix A.5 for the proof.

As $d$ cannot distinguish between points within the same polytope, it has a
resolution similar to the size of the polytope. For discriminating between two
points within the same polytope, we fit a Gaussian kernel within the polytope
(described above). As $h_{n}\to 0$, the resolution for Geodesic distance
improves. Therefore, when we use geodesic distance in Equation 4, Theorem 4.1,
4.3 hold asymptotically for our approach. In the next section, we will
empirically show that using Geodesic distance scales better with higher
dimension compared to that of Euclidean distance.

Figure 2:  Performance summary of KDF and KDN on OpenML-CC18 data suite. The
dark curve in the middle shows the median of performance on $45$ datasets. The
shaded region shows the error bar consisting of the $25$-th and the $75$-th
percentile of the performance statistics. We conducted Wilcoxon rank-sum test
using improvements in different datasets at the final sample size as test
statistics. Each panel shows a p-value pair for each basline approach where
the first and second p-values correspond to testing whether baseline model
performs worse and better than our approach respectively. Left: KDF and KDN
maintains or improves performance of their parent algorithms for
classification. Middle: KDF and KDN significantly improve the in-distribution
calibration for random forest and ReLU-net. Right: Both of the proposed
approaches yield highly calibrated confidence in the OOD region.

##  6 Experimental Results

We conduct several experiments on simulated, OpenML-CC18 [23]
111https://www.openml.org/s/99 and vision benchmark datasets to gain insights
on the finite sample performance of KDF and KDN. The details of the simulation
datasets and hyperparameters used for all the experiments are provided in
Appendix B. For Trunk simulation dataset, we follow the simulation setup
proposed by Trunk [24] which was designed to demonstrate ‘curse of
dimensionality’. In the Trunk simulation, a binary class dataset is used where
each class is sampled from a Gaussian distribution with higher dimensions
having increasingly less discriminative information. We use both Euclidean and
Geodesic distance in (4) on simulation datasets and use only Geodesic distance
for benchmark datasets. For the simulation setups, we use classification
error, hellinger distance [25, 26] from the true class conditional posteriors
and mean max confidence or posterior [4] as performance statistics. While
measuring in-distribution calibration for the datasets in OpenML-CC18 data
suite, as we do not know the true distribution, we used maximum calibration
error as defined by Guo et al. [18] with a fixed bin number of $R=15$ across
all the datasets. Given $n$ OOD samples, we define OOD calibration error (OCE)
to measure OOD performance for the benchmark datasets as:

(19)
$\text{OCE}=\frac{1}{n}\sum_{i=1}^{n}\left|\max_{y\in\mathcal{Y}}(\hat{P}_{Y|X}(y|\mathbf{x}_{i}))-\max_{y\in\mathcal{Y}}(\hat{P}_{Y}(y))\right|.$

For the tabular and the vision datasets, we have used ID calibration
approaches, such as Isotonic [15, 16] and Sigmoid [17] regression, as
baselines. Additionally, for the vision benchmark dataset, we provide results
with an OOD calibration approach ACET [4], and report traditionally used OOD
detection statistics such as area under curve (AUROC) and false positive rate
at $95\%$ precision (FPR@95). For each approach, $90\%$ of the training data
was used to fit the model and the rest of the data was used to calibrate the
model.

Table 1: Performance Summary on different Vision Datasets for Different
Algorithms When Trained on CIFAR-10

Dataset | Statistics | Parent | KDN | Isotonic | Sigmoid | ACET
---|---|---|---|---|---|---
CIFAR-10 | Accuracy$(\%)$ | $83.70$ | $83.60$ | $83.68$ | $83.70$ | $80.17$
| MCE | $0.13$ | $0.02$ | $0.02$ | $0.04$ | $0.17$
CIFAR-100 | OCE | $0.78$ | $0.52$ | $0.53$ | $0.69$ | $0.83$
| AUROC (FPR@95) | $0.78(0.83)$ | $0.69(0.74)$ | $0.79(0.77)$ | $0.67(0.83)$ | $0.72(0.86)$
SVHN | OCE | $0.73$ | $0.52$ | $0.43$ | $0.66$ | $0.78$
| AUROC (FPR@95) | $0.84(0.74)$ | $0.86(0.81)$ | $0.87(0.59)$ | $0.67(0.75)$ | $0.83(0.74)$
Noise | OCE | $0.84$ | $0.00$ | $0.68$ | $0.72$ | $0.00$
| AUROC (FPR@95) | $0.71(0.92)$ | $0.98(0.0)$ | $0.71(0.96)$ | $0.40(0.92)$ | $1.00(0.0)$

###  6.1 Simulation Study

Figure 1 leftmost column shows $10000$ training samples with $5000$ samples
per class sampled within the region $[-1,1]\times[-1,1]$ from the six
simulation setups described in Appendix B. Therefore, the empty annular region
between $[-1,1]\times[-1,1]$ and $[-2,2]\times[-2,2]$ is the low density or
OOD region in Figure 1. Figure 1 quantifies the performance of the algorithms
which are visually represented in Appendix Figure 3. KDF and KDN maintain
similar classification accuracy to those of their parent algorithms. We
measure hellinger distance from the true distribution for increasing training
sample size within $[-1,1]\times[-1,1]$ region as a statistics for in-
distribution calibration. Column $3$ of left and Column $2$ of right block in
Figure 1 show KDF and KDN are better at estimating the high density region of
training distribution compared to their parent methods. In all of the
simulations, using geodesic distance measure results in better performance
compared to those while using Euclidean distance. For measuring OOD
performance, we keep the training sample size fixed at $1000$ and normalize
the training data by the maximum of their $l2$ norm so that the training data
is confined within a unit circle. For inference, we sample $1000$ inference
points uniformly from a circle where the circles have increasing radius and
plot mean max posterior for increasing distance from the origin. Therefore,
for distance up to $1$ we have in-distribution samples and distances farther
than $1$ can be considered as OOD region. As shown in Column $3$ of Figure 1,
mean max posteriors or confidence for KDF and KDN converge to the maximum of
the class priors, i.e., $0.5$ as we go farther away from the training data
origin.

Row $6$ of Figure 1 shows KDF-Geodesic and KDN-Geodesic scale better with
higher dimensions compared to their Euclidean counterpart algorithms
respectively.

###  6.2 OpenML-CC18 Data Study

We used OpenML-CC18 data suite for tabular benchmark dataset study. We exclude
any dataset which contains categorical features or NaN values 222We also
excluded the dataset with dataset id $23517$ as we could not achieve better
than chance accuracy using RF and DN on that dataset. and conduct our
experiments on $45$ datasets with varying dimensions and sample sizes. For the
OOD experiments, we follow a similar setup as that of the simulation data. We
normalize the training data by their maximum $l2$ norm and sample $1000$
testing samples uniformly from each hypersphere where the hyperspheres have
increasing radius starting from $1$ to $5$. For each dataset, we measure
improvement with respect to the parent algorithm as:

(20) $\frac{\mathcal{E}_{p}-\mathcal{E}_{M}}{\mathcal{E}_{p}},$

where $\mathcal{E}_{p}=$classification error, MCE or OCE for the parent
algorithm and $\mathcal{E}_{M}$ represents the performance of the approach in
consideration. Note that positive improvement implies the corresponding
approach performs better than the parent approach. We report the median of
improvement on different datasets along with the error bar in Figure 2. The
extended results for each dataset is shown separately in the appendix. Figure
2 left column shows on average KDF and KDN has nearly similar or better
classification accuracy to their respective parent algorithm whereas Isotonic
and Sigmoid regression have lower classification accuracy most of the cases.
However, according to Figure 2 middle column, KDF and KDN have similar in-
distribution calibration performance to the other baseline approaches. Most
interestingly, Figure 2 right column shows that KDN and KDF improves OOD
calibration of their respective parent algorithms by a huge margin while the
baseline approaches completely fails to address the OOD calibration problem.

###  6.3 Vision Benchmark Datasets Study

For vision experiments, we trained a Resnet-20 model on CIFAR-10 dataset and
calibrated it using different ID and OOD calibration approaches. As shown in
Table 1, ID calibration, i.e., Isotonic and Sigmoid regression have relatively
better calibration for near OOD datasets like CIFAR-100 and SVHN compared to
OOD calibration approach ACET. However, ID calibration approaches have worse
far OOD (noise) performance compared to ACET. Interestingly, KDN achieves a
compromise between ID and OOD calibration approaches by treating the
calibration problem as a continuum rather than two separate problems.

As the results suggest, KDN effectively solves ID and far OOD calibration
problem. However, for improving near OOD performance, it requires further
work. For example, as shown by Fort et al. [27], near OOD performance for KDN
may be improved further by using a few known outliers. Another possible
improvement can be achieved by using other kernels which can capture local
correlation between the samples within a polytope better than the Gaussian
kernel. However, for keeping the focus of the paper concise, we will pursue
near OOD calibration in our future work.

##  7 Discussion

In this paper, we convert deep discriminative network to kernel density
network by replacing the affine function over the polytopes in the
discriminative networks with a Gaussian kernel. This replacement of affine
function results in better in- and out-of-distribution calibration for our
proposed approaches while having classification accuracy similar to or better
than the parent algorithm. Theoretically, we show under certain conditions our
approaches asymptotically converge to the true training distribution and this
establishes confidence calibration for learning algorithms in in- and out-of-
distribution regions as a continuum rather than two different problems.

## Acknowledgements

The authors thank the support of the NSF-Simons Research Collaborations on the
Mathematical and Scientific Foundations of Deep Learning (NSF grant 2031985)
and THEORINET. This work is graciously supported by the Defense Advanced
Research Projects Agency (DARPA) Lifelong Learning Machines program through
contracts FA8650-18-2-7834 and HR0011-18-2-0025. Research was partially
supported by funding from Microsoft Research and the Kavli Neuroscience
Discovery Institute.

## References

  * Guo et al. [2017a] Chuan Guo, Geoff Pleiss, Yu Sun, and Kilian Q. Weinberger.  On calibration of modern neural networks.  In Doina Precup and Yee Whye Teh, editors, _Proceedings of the 34th International Conference on Machine Learning_ , volume 70 of _Proceedings of Machine Learning Research_ , pages 1321–1330. PMLR, 06–11 Aug 2017a.
  * Kristiadi et al. [2020] Agustinus Kristiadi, Matthias Hein, and Philipp Hennig.  Being bayesian, even just a bit, fixes overconfidence in ReLU networks.  In Hal Daumé III and Aarti Singh, editors, _Proceedings of the 37th International Conference on Machine Learning_ , volume 119 of _Proceedings of Machine Learning Research_ , pages 5436–5446. PMLR, 13–18 Jul 2020.
  * Xu et al. [2021a] Haoyin Xu, Kaleab A. Kinfu, Will LeVine, Sambit Panda, Jayanta Dey, Michael Ainsworth, Yu-Chung Peng, Madi Kusmanov, Florian Engert, Christopher M. White, Joshua T. Vogelstein, and Carey E. Priebe.  When are Deep Networks really better than Decision Forests at small sample sizes, and how?  _arXiv preprint arXiv:2108.13637_ , 2021a.
  * Hein et al. [2019] Matthias Hein, Maksym Andriushchenko, and Julian Bitterwolf.  Why relu networks yield high-confidence predictions far away from the training data and how to mitigate the problem.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 41–50, 2019.
  * Meinke et al. [2021] Alexander Meinke, Julian Bitterwolf, and Matthias Hein.  Provably robust detection of out-of-distribution data (almost) for free.  _arXiv preprint arXiv:2106.04260_ , 2021.
  * Liang et al. [2017] Shiyu Liang, Yixuan Li, and Rayadurgam Srikant.  Enhancing the reliability of out-of-distribution image detection in neural networks.  _arXiv preprint arXiv:1706.02690_ , 2017.
  * Nandy et al. [2020] Jay Nandy, Wynne Hsu, and Mong Li Lee.  Towards maximizing the representation gap between in-domain & out-of-distribution examples.  _Advances in Neural Information Processing Systems_ , 33:9239–9250, 2020.
  * Wan et al. [2018] Weitao Wan, Yuanyi Zhong, Tianpeng Li, and Jiansheng Chen.  Rethinking feature distribution for loss functions in image classification.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 9117–9126, 2018.
  * DeVries and Taylor [2018] Terrance DeVries and Graham W Taylor.  Learning confidence for out-of-distribution detection in neural networks.  _arXiv preprint arXiv:1802.04865_ , 2018.
  * Hendrycks et al. [2018] Dan Hendrycks, Mantas Mazeika, and Thomas Dietterich.  Deep anomaly detection with outlier exposure.  _arXiv preprint arXiv:1812.04606_ , 2018.
  * Ren et al. [2019] Jie Ren, Peter J Liu, Emily Fertig, Jasper Snoek, Ryan Poplin, Mark Depristo, Joshua Dillon, and Balaji Lakshminarayanan.  Likelihood ratios for out-of-distribution detection.  _Advances in neural information processing systems_ , 32, 2019.
  * Kingma et al. [2019] Diederik P Kingma, Max Welling, et al.  An introduction to variational autoencoders.  _Foundations and Trends® in Machine Learning_ , 12(4):307–392, 2019.
  * Goodfellow et al. [2020] Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative adversarial networks.  _Communications of the ACM_ , 63(11):139–144, 2020.
  * Nalisnick et al. [2018] Eric Nalisnick, Akihiro Matsukawa, Yee Whye Teh, Dilan Gorur, and Balaji Lakshminarayanan.  Do deep generative models know what they don’t know?  _arXiv preprint arXiv:1810.09136_ , 2018.
  * Zadrozny and Elkan [2001] Bianca Zadrozny and Charles Elkan.  Obtaining calibrated probability estimates from decision trees and naive bayesian classifiers.  In _Icml_ , volume 1, pages 609–616, 2001.
  * Caruana [2004] R Caruana.  Predicting good probabilities with supervised learning.  In _Proceedings of NIPS 2004 Workshop on Calibration and Probabilistic Prediction in Supervised Learning_ , 2004.
  * Platt et al. [1999] John Platt et al.  Probabilistic outputs for support vector machines and comparisons to regularized likelihood methods.  _Advances in large margin classifiers_ , 10(3):61–74, 1999.
  * Guo et al. [2017b] Chuan Guo, Geoff Pleiss, Yu Sun, and Kilian Q Weinberger.  On calibration of modern neural networks.  In _International conference on machine learning_ , pages 1321–1330. PMLR, 2017b.
  * Guo et al. [2019] Richard Guo, Ronak Mehta, Jesus Arroyo, Hayden Helm, Cencheng Shen, and Joshua T Vogelstein.  Estimating information-theoretic quantities with uncertainty forests.  _arXiv_ , pages arXiv–1907, 2019.
  * Kull et al. [2019] Meelis Kull, Miquel Perello Nieto, Markus Kängsepp, Telmo Silva Filho, Hao Song, and Peter Flach.  Beyond temperature scaling: Obtaining well-calibrated multi-class probabilities with dirichlet calibration.  _Advances in neural information processing systems_ , 32, 2019.
  * Xu et al. [2021b] Haoyin Xu, Kaleab A Kinfu, Will LeVine, Sambit Panda, Jayanta Dey, Michael Ainsworth, Yu-Chung Peng, Madi Kusmanov, Florian Engert, Christopher M White, et al.  When are deep networks really better than decision forests at small sample sizes, and how?  _arXiv preprint arXiv:2108.13637_ , 2021b.
  * Schölkopf [2000] Bernhard Schölkopf.  The kernel trick for distances.  _Advances in neural information processing systems_ , 13, 2000.
  * Bischl et al. [2017] Bernd Bischl, Giuseppe Casalicchio, Matthias Feurer, Pieter Gijsbers, Frank Hutter, Michel Lang, Rafael G Mantovani, Jan N van Rijn, and Joaquin Vanschoren.  Openml benchmarking suites.  _arXiv preprint arXiv:1708.03731_ , 2017.
  * Trunk [1979] Gerard V Trunk.  A problem of dimensionality: A simple example.  _IEEE Transactions on pattern analysis and machine intelligence_ , (3):306–307, 1979.
  * Kailath [1967] Thomas Kailath.  The divergence and bhattacharyya distance measures in signal selection.  _IEEE transactions on communication technology_ , 15(1):52–60, 1967.
  * Rao [1995] C Radhakrishna Rao.  A review of canonical coordinates and an alternative to correspondence analysis using hellinger distance.  _Qüestiió: quaderns d’estadística i investigació operativa_ , 1995.
  * Fort et al. [2021] Stanislav Fort, Jie Ren, and Balaji Lakshminarayanan.  Exploring the limits of out-of-distribution detection.  _Advances in Neural Information Processing Systems_ , 34:7068–7081, 2021.
  * Perry et al. [2019] Ronan Perry, Ronak Mehta, Richard Guo, Eva Yezerets, Jesús Arroyo, Mike Powell, Hayden Helm, Cencheng Shen, and Joshua T Vogelstein.  Random forests for adaptive nearest neighbor estimation of information-theoretic quantities.  _arXiv preprint arXiv:1907.00325_ , 2019.
  * Athey et al. [2019] Susan Athey, Julie Tibshirani, and Stefan Wager.  Generalized random forests.  2019\.
  * Hornik et al. [1989] Kurt Hornik, Maxwell Stinchcombe, and Halbert White.  Multilayer feedforward networks are universal approximators.  _Neural networks_ , 2(5):359–366, 1989.

##  Appendix A Proofs

###  A.1 Proof of Theorem 4.1

Consider the density estimator for the samples populating the polytopes
$\hat{f}$. Let $n$ be the total number of samples and $n_{r}$ be the number of
data points within polytope $Q_{r}$. As stated in the theorem statement, we
make the following assumptions:

  1. 1.

The polytope bandwidth $h_{n}\to 0$ as $n\to\infty$.

  2. 2.

$n$ grows faster than the shrinkage of $h_{n}$, i.e., $n\cdot h_{n}\to\infty$
as $h_{n}\to 0$ in probability.

For simplicity, we first explore the one-dimensional distribution. The
derivation can be readily extended to multi-dimensional scenarios. We consider
any Gaussian kernel $\mathcal{G}^{\prime}()$ with parameters chosen
independently of the data satisfying two conditions. In lemma A.1 and A.3, we
will show that the aforementioned class of estimators is consistent. The
conditions for choosing the Gaussian kernel parameters are:

  1. 1.

The center of the kernel can be any point $z_{r}$ within the polytope $Q_{r}$
as $n\to\infty$,

  2. 2.

The kernel bandwidth $\sigma_{r}$ is any non-negative number always bounded by
the polytope bandwidth $h_{n}$ as $n\to\infty$, i.e., $\sigma_{r}=C_{r}h_{n}$,
where $0<C_{r}<1$.

Now the class conditional density estimate at a point $x$ can be written as:

(21)
$\hat{f}(x)=\frac{1}{n}\sum_{r\in\mathcal{P}}n_{r}\mathcal{G}^{\prime}(x;\mu_{r},\Sigma_{r})\mathbbm{1}(r=r^{*}_{x}).$

######  Lemma A.1.

The class of estimators in (21) is an asymptotically unbiased class of
estimators of the true density $f$.

######  Proof A.2.

The polytope sample counts $n_{r}$ can be considered as binomially
distributed: $n_{r}\sim B(n,P_{r})$, where $P_{r}=\int_{Q_{r}}dF$ is the
probability of finding a training sample within the polytope $r$ and $F$ is
the cumulative density function associated with the density $f$. This allows
us to write : $\mathbb{E}[n_{r}]=n\leavevmode\nobreak\ P_{r}$. Using the mean
value theorem, we have that $P_{r}=h_{n}f(q_{r})$, for some $q_{r}\in Q_{r}$.
Note that, if we consider the multi-dimensional scenario for the mean value
theorem here, the proof can be easily generalized for multi-dimensional case.

Now consider the expectation of $\hat{f}$ with respect to the training
distribution:

(22) $\displaystyle\mathbb{E}[\hat{f}(x)]$
$\displaystyle=\sum_{r\in\mathcal{P}}\frac{n\leavevmode\nobreak\
P_{r}\leavevmode\nobreak\
\mathcal{G}^{\prime}(x;\mu_{r},\Sigma_{r})\mathbbm{1}(r=r^{*}_{x})}{n}$ (23)
$\displaystyle=\sum_{r\in\mathcal{P}}h_{n}f(q_{r})\mathcal{G}^{\prime}(x;\mu_{r},\Sigma_{r})\mathbbm{1}(r=r^{*}_{x}).$

Note that we are given the partitions of the feature space and we choose
Gaussian parameters such that they are independent of the training data. Now,
as $n\to\infty$, $h_{n}\to 0$, and $h_{n}$ can be considered an infinitesimal
measure $dz_{r}$. Furthermore, as the bandwidth of the Gaussian is limited by
the polytope bandwidth and the area under the Gaussian is $1$, the kernel
summation
$\sum_{r\in\mathcal{P}}\mathcal{G}^{\prime}(x,\mu_{r},\Sigma_{r})\mathbbm{1}(r=r^{*}_{x})$
becomes a dirac delta function evaluated at $x$ as $h_{n}\to 0$. Therefore, in
the limiting conditions,
$\sum_{r\in\mathcal{P}}\mathcal{G}^{\prime}(x,\mu_{r},\Sigma_{r})\mathbbm{1}(r=r^{*}_{x})\to\delta(z_{r}-x)$
and $f(q_{r})\to f(z_{r})$. Here $z_{r}$ is a point such that $z_{r}\in Q_{r}$
for all $n$ and $\sigma_{r}<h_{n}$. Therefore, there exists a small positive
number $\epsilon$ for which we can rewrite 23 as:

(24) $\displaystyle\mathbb{E}[\hat{f}(x)]$
$\displaystyle=\int_{x-\varepsilon}^{x+\varepsilon}f(z_{r})\delta(z_{r}-x)dz_{r}$
$\displaystyle=f(x)$

Therefore, $\hat{f}(x)$ is an asymptotically unbiased estimator of $f(x)$.

######  Lemma A.3.

The variance of the class of estimators in (21) asymptotically goes to $0$.

######  Proof A.4.

For binomially distributed samples $n_{r}$, we can write,
$\mathbb{V}ar[n_{r}]=nP_{r}(1-P_{r})$. Therefore, we can estimate the variance
of $\hat{f}(x)$ as:

(25)
$\displaystyle\mathbb{V}ar[\hat{f}(x)]=\sum_{r\in\mathcal{P}}\frac{nP_{r}(1-P_{r})}{n^{2}}\\{\mathcal{G}^{\prime}(x;\mu_{r},\Sigma_{r})\mathbbm{1}(r=r^{*}_{x})\\}^{2}$
$\displaystyle=\frac{h_{n}f(z_{r^{*}_{x}})(1-h_{n}f(z_{r^{*}_{x}}))}{n}\\{\frac{1}{\sqrt{2\pi}\sigma_{r^{*}_{x}}}\exp{(-\frac{(x-\mu_{r^{*}_{x}})^{2}}{2\sigma_{r^{*}_{x}}^{2}})}\\}^{2}$
$\displaystyle=\frac{h_{n}f(z_{r^{*}_{x}})(1-h_{n}f(z_{r^{*}_{x}}))}{nC_{r^{*}_{x}}^{2}h_{n}^{2}}\\{\frac{1}{\sqrt{2\pi}}\exp{(-\frac{(x-\mu_{r^{*}_{x}})^{2}}{2\sigma_{r^{*}_{x}}^{2}})}\\}^{2}$
$\displaystyle\leq\frac{f(z_{r^{*}_{x}})(1-h_{n}f(z_{r^{*}_{x}}))}{2\pi
C_{r^{*}_{x}}^{2}(nh_{n})}$

Equation 25 offers several interesting insights:

  1. 1.

As $h_{n}\to 0$, $z_{r^{*}_{x}}\to x$ and
$f(z_{r^{*}_{x}})(1-h_{n}f(z_{r^{*}_{x}}))\to f(x)$. Therefore, the variance
becomes directly dependent on $f(x)$. There is more variability at regions
with higher density.

  2. 2.

The variance is also higher if $C_{r^{*}_{x}}^{2}\ll 1$. Therefore,
$C_{r^{*}}$ should be as close to $1$ as possible and $nh_{n}$ should be as
high as possible for lower variance. Moreover, variance can be reduced at the
expense of higher bias with $C_{r^{*}_{x}}>1$. Most importantly note that, the
effect of $C_{r^{*}_{x}}$ cancels out in the numerator and the denominator
while estimating the posteriors in Equation 7.

  3. 3.

As $h_{n}\to 0$ the estimation becomes unbiased (see Equation 23), but the
estimation variance becomes unbounded. Therefore, for bounded and decreasing
variance the condition $nh_{n}\to\infty$ is necessary.

Lemma A.1 and A.3 together completes the proof of Theorem 4.1.

###  A.2 Proof of Theorem 4.3

We first expand $\hat{g}_{y}(\mathbf{x})$:

$\displaystyle\hat{g}_{y}(\mathbf{x})$
$\displaystyle=\frac{\hat{f}_{y}(\mathbf{x})\hat{P}_{Y}(y)}{\sum_{k=1}^{K}\hat{f}_{k}(x)\hat{P}_{Y}(k)}$
$\displaystyle=\frac{\tilde{f}_{y}(\mathbf{x})\hat{P}_{Y}(y)+\frac{b}{\log(n)}\hat{P}_{Y}(y)}{\sum_{k=1}^{K}(\hat{f}_{k}(\mathbf{x})\hat{P}_{Y}(k)+\frac{b}{\log(n)}\hat{P}_{Y}(k))}$

As the inference point $\mathbf{x}$ becomes more distant from training samples
(and more distant from all of the Gaussian centers), we have that
$\mathcal{G}(\mathbf{x},\hat{\mu}_{r},\hat{\Sigma}_{r})$ becomes smaller.
Thus, $\forall y,\tilde{f}_{y}(\mathbf{x})$ shrinks. More formally, $\forall
y$,

$\lim_{d_{\mathbf{x}}\to\infty}\tilde{f}_{y}(\mathbf{x})=0$

We can use this result to then examine the limiting behavior of our posteriors
as the inference point $\mathbf{x}$ becomes more distant from the training
data:

$\displaystyle\lim_{d_{\mathbf{x}}\to\infty}\hat{g}_{y}(\mathbf{x})$
$\displaystyle=\lim_{d_{\mathbf{x}}\to\infty}\frac{\tilde{f}_{y}(\mathbf{x})\hat{P}_{Y}(y)+\frac{b}{\log(n)}\hat{P}_{Y}(y)}{\sum_{k=1}^{K}(\tilde{f}_{k}(\mathbf{x})\hat{P}_{Y}(k)+\frac{b}{\log(n)}\hat{P}_{Y}(k))}$
$\displaystyle=\frac{(\lim_{d_{\mathbf{x}}\to\infty}\tilde{f}_{y}(\mathbf{x}))\hat{P}_{Y}(y)+\frac{b}{\log(n)}\hat{P}_{Y}(y)}{\sum_{k=1}^{K}(\lim_{d_{\mathbf{x}}\to\infty}\tilde{f}_{k}(\mathbf{x}))\hat{P}_{Y}(k)+\frac{b}{\log(n)}\hat{P}_{Y}(k))}$
$\displaystyle=\frac{\hat{P}_{Y}(y)}{\sum_{k=1}^{K}\hat{P}_{Y}(k)}$
$\displaystyle=\hat{P}_{Y}(y)$

###  A.3 Proof of Corollary 4.7

Consider a class of random forest algorithms and a training distribution which
satisfies all the conditions mentioned in Perry et al. [28], Athey et al.
[29]. The above class of algorithms partitions the feature space so that as
$n\to\infty$, the average polytope bandwidth $h_{n}\to 0$ and
$nh_{n}\to\infty$.

###  A.4 Proof of Corollary 4.9

Consider a training distribution that satisfies all the conditions mentioned
in [28, 29]. Using the universal approximation theorems for neural networks
[30], we can say the number of parameters for deep-net should increase as
$n\to\infty$ for an unbiased estimation of the posteriors. In other words, the
number of affine function should increase resulting in a finer partition of
the feature space, i.e., $h_{n}\to 0$ and $nh_{n}\to\infty$ for the polytopes.

###  A.5 Proof of Corollary 5.1

For proving that $d$ is a valid distance metric for $\mathcal{P}_{n}$, we need
to prove the following four statements:

  1. 1.

$d(r,s)=0$ when $r=s$.
Proof: By definition, $\mathcal{K}(r,s)=1$ and $d(r,s)=0$ when $r=s$.

  2. 2.

$d(r,s)>0$ when $r\neq s$.
Proof: By definition, $0\leq\mathcal{K}(r,s)<1$ and $d(r,s)>0$ for $r\neq s$.

  3. 3.

$d$ is symmetric, i.e., $d(r,s)=d(s,r)$.
Proof: By definition, $\mathcal{K}(r,s)=\mathcal{K}(s,r)$ which implies
$d(r,s)=d(s,r)$.

  4. 4.

$d$ follows the triangle inequality, i.e., for any three polytopes
$Q_{r},Q_{s},Q_{t}\in\mathcal{P}_{n}$: $d(r,t)\leq d(r,s)+d(s,t)$.
Proof: Let $\mathcal{A}_{r}$ denote the set of activation modes in a ReLU-net
and the set of leaf nodes in a random forest for a particular polytope $r$.
$N$ is the total number of possible activation paths in a ReLU-net or total
trees in a random forest. We can write:

(26) $\displaystyle N$ $\displaystyle\geq
n((\mathcal{A}_{r}\cap\mathcal{A}_{s})\cup(\mathcal{A}_{s}\cap\mathcal{A}_{t}))$
$\displaystyle=n(\mathcal{A}_{r}\cap\mathcal{A}_{s})+n(\mathcal{A}_{s}\cap\mathcal{A}_{t})-n(\mathcal{A}_{r}\cap\mathcal{A}_{s}\cap\mathcal{A}_{t})$
$\displaystyle\geq
n(\mathcal{A}_{r}\cap\mathcal{A}_{s})+n(\mathcal{A}_{s}\cap\mathcal{A}_{t})-n(\mathcal{A}_{r}\cap\mathcal{A}_{t})$

Rearranging the above equation, we get:

$\displaystyle N-n(\mathcal{A}_{r}\cap\mathcal{A}_{t})\leq
N-n(\mathcal{A}_{r}\cap\mathcal{A}_{s})+N-n(\mathcal{A}_{s}\cap\mathcal{A}_{t})$
$\displaystyle\implies 1-\frac{n(\mathcal{A}_{r}\cap\mathcal{A}_{t})}{N}\leq
1-\frac{n(\mathcal{A}_{r}\cap\mathcal{A}_{s})}{N}+1$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ -\frac{n(\mathcal{A}_{s}\cap\mathcal{A}_{t})}{N}$ (27)
$\displaystyle\implies d(r,t)\leq d(r,s)+d(s,t)$

##  Appendix B Simulations

Figure 3:  Visualization of true and estimated posteriors for class 0 from
five binary class simulation experiments. Column 1: 10,000 training points
with 5,000 samples per class sampled from 6 different simulation setups for
binary class classification. Trunk simulation is shown for two dimensional
case. The class labels are indicated by yellow and blue colors. Column 2-8:
True and estimated class conditional posteriors from different approaches. The
posteriors estimated from KDN and KDF are better calibrated for both in- and
out-of-distribution regions compared to those of their parent algorithms.

We construct six types of binary class simulations:

  * •

Gaussian XOR is a two-class classification problem with equal class priors.
Conditioned on being in class 0, a sample is drawn from a mixture of two
Gaussians with means $\pm[0.5,-0.5]^{\top}$ and standard deviations of $0.25$.
Conditioned on being in class 1, a sample is drawn from a mixture of two
Gaussians with means $\pm[0.5,-0.5]^{\top}$ and standard deviations of $0.25$.

  * •

Spiral is a two-class classification problem with the following data
distributions: let $K$ be the number of classes and $S\sim$
multinomial$(\frac{1}{K}\vec{1}_{K},n)$. Conditioned on $S$, each feature
vector is parameterized by two variables, the radius $r$ and an angle
$\theta$. For each sample, $r$ is sampled uniformly in $[0,1]$. Conditioned on
a particular class, the angles are evenly spaced between
$\frac{4\pi(k-1)t_{K}}{K}$ and $\frac{4\pi(k)t_{K}}{K}$, where $t_{K}$
controls the number of turns in the spiral. To inject noise along the spirals,
we add Gaussian noise to the evenly spaced angles
$\theta^{\prime}:\theta=\theta^{\prime}+\mathcal{N}(0,0.09)$. The observed
feature vector is then $(r\;\cos(\theta),r\;\sin(\theta))$.

  * •

Circle is a two-class classification problem with equal class priors.
Conditioned on being in class 0, a sample is drawn from a circle centered at
$(0,0)$ with a radius of $r=0.75$. Conditioned on being in class 1, a sample
is drawn from a circle centered at $(0,0)$ with a radius of $r=1$, which is
cut off by the region bounds. To inject noise along the circles, we add
Gaussian noise to the circle radii
$r^{\prime}:r=r^{\prime}+\mathcal{N}(0,0.01)$.

  * •

Sinewave is a two-class classification problem based on sine waves.
Conditioned on being in class 0, a sample is drawn from the distribution
$y=\cos(\pi x)$. Conditioned on being in class 1, a sample is drawn from the
distribution $y=\sin(\pi x)$. We inject Gaussian noise to the sine wave
heights $y^{\prime}:y=y^{\prime}+\mathcal{N}(0,0.01)$.

  * •

Polynomial is a two-class classification problem with the following data
distributions: $y=x^{a}$. Conditioned on being in class 0, a sample is drawn
from the distribution $y=x^{1}$. Conditioned on being in class 1, a sample is
drawn from the distribution $y=x^{3}$. Gaussian noise is added to variables
$y^{\prime}:y=y^{\prime}+\mathcal{N}(0,0.01)$.

  * •

Trunk is a two-class classification problem with gradually increasing
dimension and equal class priors. The class conditional probabilities are
Gaussian:

$\displaystyle P(X|Y=0)$ $\displaystyle=\mathcal{G}(\mu_{1},I),$
$\displaystyle P(X|Y=1)$ $\displaystyle=\mathcal{G}(\mu_{2},I),$

where $\mu_{1}=\mu,\mu_{2}=-\mu$, $\mu$ is a $d$ dimensional vector whose
$i$-th component is $(\frac{1}{i})^{1/2}$ and $I$ is $d$ dimensional identity
matrix.

Table 2: Hyperparameters for RF and KDF. Hyperparameters | Value
---|---
n_estimators | $500$
max_depth | $\infty$
min_samples_leaf | $1$
$b$ (bias) | $\exp(-10^{\sqrt{d}})$
$\lambda$ | $1\times 10^{-6}$
Table 3: Hyperparameters for ReLU-net and KDN. Hyperparameters | Value
---|---
number of hidden layers | $4$
nodes per hidden layer | $1000$
optimizer | Adam
learning rate | $3\times 10^{-4}$
$b$ (bias) | $\exp(-10^{\sqrt{d}})$
$\lambda$ | $1\times 10^{-6}$

##  Appendix C Pseudocodes

We provide the pseudocode for our porposed algorithms in Algorithm 1, 2 and 3.

Algorithm 1 Fit a KGX model.

1:

2:(1) $\theta$ $\triangleright$ Parent learner (random forest or deep network
model)

3:(2) $\mathcal{D}_{n}=(\mathbf{X},\mathbf{y})\in\mathbb{R}^{n\times
d}\times\\{1,\ldots,K\\}^{n}$ $\triangleright$ Training data

4:$\mathcal{G}$ $\triangleright$ a KGX model

5:function KGX.fit($\theta,\mathbf{X},\mathbf{y}$)

6: for $i=1,\ldots,n$ do $\triangleright$ Iterate over the dataset to
calculate the weights

7: for $j=1,\ldots,n$ do

8: $w_{ij}\leftarrow$ computeWeights($\mathbf{x}_{i},\mathbf{x}_{j},\theta$)

9: end for

10: end for

11:

12:

13: $\\{Q_{r}\\}_{r=1}^{\tilde{p}}\leftarrow$ getPolytopes($\mathbf{w}$)
$\triangleright$ Identify the polytopes by clustering the samples with similar
weight

14:

15: $\mathcal{G}.\\{\tilde{w}_{k}\\}_{k=1}^{K}\leftarrow 0$ $\triangleright$
Initialize the counts for each class

16: for $r=1,\ldots,\tilde{p}$ do $\triangleright$ Iterate over each polytope

17: for $k=1,\ldots,K$ do

18: $\mathcal{G}.\tilde{w}_{rk}\leftarrow$
countWeights($\\{\mathbf{w}_{rs}\\}_{s=1}^{\tilde{p}},k$) $\triangleright$
$w_{rk}$ is the number of weighted input samples in $Q_{r}$ with label $k$

19:
$\mathcal{G}.\tilde{w}_{k}\leftarrow\mathcal{G}.\tilde{w}_{k}+\mathcal{G}.\tilde{w}_{rk}$
$\triangleright$ Update the total count for each class

20: end for

21: $\mathcal{G}.\hat{\mu}_{r},\mathcal{G}.\hat{\Sigma}_{r}\leftarrow$
estimateParameters($\mathbf{X},\\{\mathbf{w}_{rs}\\}_{s=1}^{\tilde{p}}$)
$\triangleright$ Fit Gaussians using weighted MLE

22: end for

23: return $\mathcal{G}$

24:end function

Algorithm 2 Computing weights in KDF

1:

2:(1) $\mathbf{x}_{i},\mathbf{x}_{j}\in\mathbb{R}^{1\times d}$
$\triangleright$ two input samples to be weighted

3:(2) $\theta$ $\triangleright$ parent random forest with $T$ trees

4:$w_{ij}\in[0,1]$ $\triangleright$ compute similarity between $i$ and $j$-th
samples.

5:function computeWeights($\mathbf{x}_{i},\mathbf{x}_{j},\theta$)

6: $\mathcal{I}_{i}\leftarrow$ pushDownTrees($\mathbf{x}_{i},\theta$)
$\triangleright$ push $\mathbf{x}_{i}$ down $T$ trees and get the leaf numbers
it end up in.

7: $\mathcal{I}_{j}\leftarrow$ pushDownTrees($\mathbf{x}_{j},\theta$)
$\triangleright$ push $\mathbf{x}_{j}$ down $T$ trees and get the leaf numbers
it end up in.

8: $l\leftarrow$ countMatches($\mathcal{I}_{i},\mathcal{I}_{j}$)
$\triangleright$ count the number of times the samples end up in the same leaf

9: $w_{ij}\leftarrow\frac{l}{T}$

10: return $w_{ij}$

11:end function

Algorithm 3 Computing weights in KDN

1:

2:(1) $\mathbf{x}_{i},\mathbf{x}_{j}\in\mathbb{R}^{1\times d}$
$\triangleright$ two input samples to be weighted

3:(2) $\theta$ $\triangleright$ parent deep-net model

4:$w_{ij}\in[0,1]$ $\triangleright$ compute similarity between $i$ and $j$-th
samples.

5:function computeWeights($\mathbf{x}_{i},\mathbf{x}_{j},\theta$)

6: $\mathcal{A}_{i}\leftarrow$ pushDownNetwork($\mathbf{x}_{i},\theta$)
$\triangleright$ get activation modes $\mathcal{A}_{i}$

7: $\mathcal{A}_{j}\leftarrow$ pushDownNetwork($\mathbf{x}_{j},\theta$)
$\triangleright$ get activation modes $\mathcal{A}_{j}$

8: $l\leftarrow$ countMatches($\mathcal{A}_{i}$, $\mathcal{A}_{j}$)
$\triangleright$ count the number of times the two samples activate the
activation paths in a similar way

9: $w_{ij}\leftarrow\frac{l}{N}$ $\triangleright$ $N$ is the total number of
activation paths

10: return $w_{ij}$

11:end function

##  Appendix D Hardware and Software Configurations

  * •

Operating System: Linux (ubuntu 20.04), macOS (Ventura 13.2.1)

  * •

VM Size: Azure Standard D96as v4 (96 vcpus, 384 GiB memory)

  * •

GPU: Apple M1 Max

  * •

Software: Python 3.8, scikit-learn $\geq$ 0.22.0, tensorflow-macos$\leq$2.9,
tensorflow-metal $\leq$ 0.5.0.

Figure 4:  Extended results on OpenML-CC18 datasets. Left: Performance
(classification error, MCE and mean max confidence) of KDF on different
Openml-CC18 datasets. Right: Performance (classification error, MCE and mean
max confidence) of KDN on different Openml-CC18 datasets.  Figure 5:  Extended
results on OpenML-CC18 datasets (continued). Left: Performance (classification
error, MCE and mean max confidence) of KDF on different Openml-CC18 datasets.
Right: Performance (classification error, MCE and mean max confidence) of KDN
on different Openml-CC18 datasets.  Figure 6:  Extended results on OpenML-CC18
datasets (continued). Left: Performance (classification error, MCE and mean
max confidence) of KDF on different Openml-CC18 datasets. Right: Performance
(classification error, MCE and mean max confidence) of KDN on different
Openml-CC18 datasets.  Figure 7:  Extended results on OpenML-CC18 datasets
(continued). Left: Performance (classification error, MCE and mean max
confidence) of KDF on different Openml-CC18 datasets. Right: Performance
(classification error, MCE and mean max confidence) of KDN on different
Openml-CC18 datasets.