Learning to Defer to Multiple Experts: Consistent Surrogate Losses,
Confidence Calibration, and Conformal Ensembles

Rajeev Verma∗ Daniel Barrejón∗ Eric Nalisnick

University of Amsterdam Universidad Carlos III de Madrid *Equal contribution,
order determined by coin flip University of Amsterdam

###### Abstract

We study the statistical properties of learning to defer (L2D) to multiple
experts. In particular, we address the open problems of deriving a consistent
surrogate loss, confidence calibration, and principled ensembling of experts.
Firstly, we derive two consistent surrogates—one based on a softmax
parameterization, the other on a one-vs-all (OvA) parameterization—that are
analogous to the single expert losses proposed by Mozannar and Sontag (2020)
and Verma and Nalisnick (2022), respectively. We then study the frameworks’
ability to estimate
$\mathbb{P}({\textnormal{m}}_{j}={\textnormal{y}}|{\bm{x}})$, the probability
that the $j$th expert will correctly predict the label for ${\bm{x}}$. Theory
shows the softmax-based loss causes mis-calibration to propagate between the
estimates while the OvA-based loss does not (though in practice, we find there
are trade offs). Lastly, we propose a conformal inference technique that
chooses a subset of experts to query when the system defers. We perform
empirical validation on tasks for galaxy, skin lesion, and hate speech
classification.

##  1 INTRODUCTION

Solving complex problems often requires the involvement of multiple experts
(Fay et al., 2006). These experts may have non-overlapping specialties, such
as in a large construction project that requires the advice of engineers,
architects, geologists, lawyers, etc. Or perhaps the difficulty of the task
requires multiple opinions, as when a team of doctors consults on a difficult
medical diagnosis. Thus, modern hybrid intelligent (HI) systems (Kamar, 2016;
Dellermann et al., 2019; Akata et al., 2020)—so called because they combine
computational models with human decision makers—need to support the
participation of multiple experts.

Learning to defer (L2D) (Madras et al., 2018) presents an elegant framework
for implementing HI systems. A rejector model acts as a meta-classifier,
predicting whether the downstream classifier or human is more likely to make
the correct decision for a given input. Yet existing L2D frameworks do not
obviously accommodate additional experts. For instance, the rejector’s job
becomes more challenging when there are multiple experts. It has two decisions
to make: _when_ to defer and _to which_ expert. The latter decision is equally
important, and thus verifying that the rejector can accurately monitor expert
quality is essential for safe and effective deployment.

In this paper, we develop the statistical foundations of multi-class L2D with
multiple experts. Specifically, we address the following open problems:
deriving a consistent surrogate loss, studying whether these systems are
confidence calibrated, and developing a principled technique for ensembling
expert decisions. The first and second contributions ensure the soundness of
the optimization problem and resulting downstream decision making. Our third
contribution, expert ensembles, follows from our study of calibration, as we
propose a conformal inference procedure for selecting a subset of experts.
Querying multiple experts may be necessary either from a standpoint of safety
or efficiency. By the latter, we mean that if the system has already incurred
the efficiency cost of deferring, then the cost of asking additional experts
may be marginal. We empirically validate our findings on the tasks of image
classification, galaxy categorization, skin lesion diagnosis, and hate speech
detection.

##  2 LEARNING TO DEFER

Mozannar and Sontag (2020) and Verma and Nalisnick (2022) proposed the only
known consistent (surrogate) loss functions for multiclass L2D. Thus we focus
on their formulations. We provide the technical background of L2D in this
section before moving on to the multi-expert setting in the subsequent
section.

##### Data

We first define the data for multiclass L2D with one expert. Let $\mathcal{X}$
denote the feature space, and let $\mathcal{Y}$ denote the label space, which
we assume to be a categorical encoding of multiple ($K$) classes.
${\mathbf{x}}_{n}\in\mathcal{X}$ denotes a feature vector, and
${\textnormal{y}}_{n}\in\mathcal{Y}$ denotes the associated class defined by
$\mathcal{Y}$ (1 of $K$). The L2D problem also assumes that we have access to
(human) expert demonstrations. Denote the expert’s prediction space as
$\mathcal{M}$, which is usually taken to be equal to the label space:
$\mathcal{M}$ = $\mathcal{Y}$. Expert demonstrations are denoted
${\textnormal{m}}_{n}\in\mathcal{M}$ for the associated features
${\mathbf{x}}_{n}$. The combined N-element training sample is
$\mathcal{D}=\\{{\bm{x}}_{n},y_{n},m_{n}\\}_{n=1}^{N}$.

##### Models

Mozannar and Sontag (2020)’s and Verma and Nalisnick (2022)’s L2D frameworks
compose two models: a classifier and a rejector (Cortes et al., 2016a, b).
Denote the classifier as $h:\mathcal{X}\rightarrow\mathcal{Y}$ and the
rejector as $r:\mathcal{X}\rightarrow\\{0,1\\}$. When $r({\mathbf{x}})=0$, the
classifier makes the decision, and when $r({\mathbf{x}})=1$, the classifier
abstains and defers the decision to the human.

##### Learning

The learning problem requires fitting both the rejector and classifier. When
the classifier makes the prediction, then the system incurs a loss
$\ell(h({\bm{x}}),y)$. When the human makes the prediction (i.e.
$r({\bm{x}})=1$), then the system incurs a loss $\ell_{\text{exp}}(m,y)$.
Using the rejector to combine these losses, we have the overall classifier-
rejector loss:

$\begin{split}L(&h,r)=\\\
&\mathbb{E}_{{\mathbf{x}},{\textnormal{y}},{\textnormal{m}}}\left[(1-r({\mathbf{x}}))\
\ell(h({\mathbf{x}}),{\textnormal{y}})\ +\ r({\mathbf{x}})\
\ell_{\text{exp}}({\textnormal{m}},{\textnormal{y}})\right]\end{split}$ (1)

where the rejector is acting as an indicator function that controls which loss
to use. While this formulation is valid for general losses, the canonical
$0-1$ loss is of special interest for classification tasks:

$\begin{split}&L_{0-1}(h,r)=\\\ &\ \ \ \ \
\mathbb{E}_{{\mathbf{x}},{\textnormal{y}},{\textnormal{m}}}\left[(1-r({\mathbf{x}}))\
\mathbb{I}[h({\mathbf{x}})\neq{\textnormal{y}}]\ +\ r({\mathbf{x}})\
\mathbb{I}[{\textnormal{m}}\neq{\textnormal{y}}]\right]\end{split}$ (2)

where $\mathbb{I}$ denotes an indicator function that checks if the prediction
and label are equal. Upon minimization, the resulting Bayes optimal classifier
and rejector satisfy:

$\begin{split}h^{*}({\bm{x}})&=\operatorname*{arg\,max}_{y\in\mathcal{Y}}\
\mathbb{P}({\textnormal{y}}=y|{\bm{x}}),\\\
r^{*}({\bm{x}})&=\mathbb{I}\left[\mathbb{P}({\textnormal{m}}={\textnormal{y}}|{\bm{x}})\geq\max_{y\in\mathcal{Y}}\mathbb{P}({\textnormal{y}}=y|{\bm{x}})\right],\end{split}$
(3)

where $\mathbb{P}({\textnormal{y}}|{\bm{x}})$ is the probability of the label
under the data generating process, and
$\mathbb{P}({\textnormal{m}}={\textnormal{y}}|{\bm{x}})$ is the probability
that the expert is correct. The expert likely will have additional knowledge
not available to the classifier, and this assumption allows the expert to
possibly outperform the Bayes optimal classifier.

##### Softmax Surrogate

The key innovation of Mozannar and Sontag (2020) is the proposal of a
consistent surrogate loss for $L_{0-1}$. They accomplish this by first
unifying the classifier and rejector via an augmented label space that
includes the rejection option. Formally, this label space is defined as
$\mathcal{Y}^{\bot}=\mathcal{Y}\cup\\{\bot\\}$ where $\bot$ denotes the
rejection option. Secondly, Mozannar and Sontag (2020) use a reduction to cost
sensitive learning that ultimately resembles the cross-entropy loss for a
softmax parameterization. Let $g_{k}:\mathcal{X}\mapsto\mathbb{R}$ for
$k\in[1,K]$ where $k$ denotes the class index, and let
$g_{\bot}:\mathcal{X}\mapsto\mathbb{R}$ denote the rejection ($\bot$) option.
These $K+1$ functions are then combined in the following softmax-parameterized
surrogate loss:

$\begin{split}{\upphi}_{\text{SM}}(g_{1}&,\ldots,g_{K},g_{\bot};{\bm{x}},y,m)=\\\
&-\log\left(\frac{\exp\\{g_{y}({\bm{x}})\\}}{\sum_{y^{\prime}\in\mathcal{Y}^{\bot}}\exp\\{g_{y^{\prime}}({\bm{x}})\\}}\right)\\\
&-\mathbb{I}[m=y]\
\log\left(\frac{\exp\\{g_{\bot}({\bm{x}})\\}}{\sum_{y^{\prime}\in\mathcal{Y}^{\bot}}\exp\\{g_{y^{\prime}}({\bm{x}})\\}}\right).\end{split}$
(4)

The intuition is that the first term maximizes the function $g_{k}$ associated
with the true label. The second term then maximizes the rejection function
$g_{\bot}$ but only if the expert’s prediction is correct. At test time, the
classifier is obtained by taking the maximum over $k\in[1,K]$:
$\hat{y}=h({\bm{x}})=\operatorname*{arg\,max}_{k\in[1,K]}g_{k}({\bm{x}})$. The
rejection function is similarly formulated as
$r({\bm{x}})=\mathbb{I}[g_{\bot}({\bm{x}})\geq\max_{k}g_{k}({\bm{x}})]$. The
minimizers $g^{*}_{1},\ldots,g^{*}_{K},g^{*}_{\bot}$ of ${\upphi}_{\text{SM}}$
also uniquely minimize $L_{0-1}(h,r)$, the $0-1$ loss from Equation 2.

##### One-vs-All Surrogate

Verma and Nalisnick (2022) propose an alternative consistent surrogate for
multi-class L2D based on a one-vs-all (OvA) formulation. Again assume we have
$K+1$ functions
$g_{1}({\mathbf{x}}),\ldots,g_{K}({\mathbf{x}}),g_{\bot}({\mathbf{x}})$ such
that $g:\mathcal{X}\mapsto\mathbb{R}$. Their one-vs-all (OvA) surrogate loss
has the point-wise form:

$\begin{split}&{\uppsi}_{\text{OvA}}(g_{1},\ldots,g_{K},g_{\bot};{\bm{x}},y,m)=\\\
&\ \ \ {\upphi}[g_{y}({\bm{x}})]+\sum_{y^{\prime}\in\mathcal{Y},y^{\prime}\neq
y}{\upphi}[-g_{y^{\prime}}({\bm{x}})]\ \ +\\\ &\ \ \
{\upphi}[-g_{\bot}({\bm{x}})]+\mathbb{I}[m=y]\left({\upphi}[g_{\bot}({\bm{x}})]-{\upphi}[-g_{\bot}({\bm{x}})]\right)\end{split}$
(5)

where ${\upphi}:\\{\pm 1\\}\times\mathbb{R}\mapsto\mathbb{R}_{+}$ is a binary
surrogate loss. For instance, when ${\upphi}$ is the logistic loss, we have
${\upphi}[f({\bm{x}})]=\log(1+\exp\\{-f({\bm{x}})\\})$. The classifier and
rejector decisions are computed exactly the same as for the softmax loss:
$h({\bm{x}})=\operatorname*{arg\,max}_{k\in[1,K]}g_{k}({\bm{x}})$,
$r({\bm{x}})=\mathbb{I}[g_{\bot}({\bm{x}})\geq\max_{k}g_{k}({\bm{x}})]$.

##  3 L2D TO MULTIPLE EXPERTS

We now turn to the multi-expert setting, deriving two consistent surrogate
losses that are analogous to Mozannar and Sontag (2020)’s and Verma and
Nalisnick (2022)’s single-expert loss functions.

##### Data

Again let ${\mathbf{x}}_{n}\in\mathcal{X}$ and
${\textnormal{y}}_{n}\in\mathcal{Y}$ be the feature and label respectively, as
defined above. Now let there be $J$ experts, and denote each expert’s
prediction space as $\mathcal{M}_{j}$ (which again we will assume is equal to
the label space: $\mathcal{M}_{j}$ = $\mathcal{Y}$ $\ \forall j$). The expert
demonstrations are denoted ${\textnormal{m}}_{n,j}\in\mathcal{M}_{j}$ for the
associated features ${\mathbf{x}}_{n}$. The combined N-element training sample
is then denoted
$\mathcal{D}=\\{{\bm{x}}_{n},y_{n},m_{n,1},\ldots,m_{n,J}\\}_{n=1}^{N}$.

##### Models

Again we use the classifier-rejector formulation (Cortes et al., 2016a, b). In
L2D with one expert, the rejector makes a binary decision—to defer or not—but
in multi-expert L2D, the rejector also must choose _to which_ expert to assign
the instance. Hence let the rejector be denoted
$r:\mathcal{X}\rightarrow\\{0,1,\ldots,J\\}$. When $r({\mathbf{x}})=0$, the
classifier makes the decision in the typical way. When $r({\mathbf{x}})=j$,
the classifier abstains and defers the decision to the $j$th expert. The
classifier in this setting ($h$) is unchanged from the single expert case.

##### Learning

Again the learning objective is the $0-1$ loss. We can re-write Equation 2 for
the multi-expert setting as:

$\begin{split}&L_{0-1}(h,r)=\mathbb{E}_{{\mathbf{x}},{\textnormal{y}},\\{{\textnormal{m}}_{j}\\}_{j=1}^{J}}\Bigg{[}\mathbb{I}[r({\mathbf{x}})=0]\
\mathbb{I}[h({\mathbf{x}})\neq{\textnormal{y}}]\\\ &\ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ +\
\sum_{j=1}^{J}\mathbb{I}[r({\mathbf{x}})=j]\
\mathbb{I}[{\textnormal{m}}_{j}\neq{\textnormal{y}}]\Bigg{]}\end{split}$ (6)

The corresponding Bayes optimal classifier and rejector are:

$\begin{split}h^{*}({\bm{x}})&=\operatorname*{arg\,max}_{y\in\mathcal{Y}}\
\mathbb{P}({\textnormal{y}}=y|{\bm{x}}),\\\
r^{*}({\bm{x}})&=\begin{cases}&0\text{ if
}\mathbb{P}({\textnormal{y}}=h^{*}({\bm{x}})|{\bm{x}})>\mathbb{P}({\textnormal{m}}_{j^{\prime}}={\textnormal{y}}|{\bm{x}})\
\ \forall j^{\prime}\\\
&\operatorname*{arg\,max}_{j\in[1,J]}\mathbb{P}({\textnormal{m}}_{j}={\textnormal{y}}|{\bm{x}})\
\ \text{ otherwise},\end{cases}\end{split}$ (7)

where $\mathbb{P}({\textnormal{y}}|{\bm{x}})$ is again the probability of the
label under the data generating process and
$\mathbb{P}({\textnormal{m}}_{j}={\textnormal{y}}|{\bm{x}})$ is the
probability that the $j$th expert is correct.

###  3.1 Softmax Surrogate Loss

Given the preceding definitions, we can now define the multi-expert analog of
the softmax-based surrogate loss. Define the augmented label space as
$\mathcal{Y}^{\bot}=\mathcal{Y}\cup\\{\bot_{1},\ldots,\bot_{J}\\}$ where
$\bot_{j}$ denotes the decision to defer to the $j$th expert. Let the
classifier be composed of $K$ functions: $g_{k}:\mathcal{X}\mapsto\mathbb{R}$
for $k\in[1,K]$ where $k$ denotes the class index. Then let the rejector be
implemented with $J$ functions: $g_{\bot,j}:\mathcal{X}\mapsto\mathbb{R}$ for
$j\in[1,J]$ where $j$ is the expert index. We propose to combine these $K+J$
functions via the following softmax-parameterized surrogate loss:

$\begin{split}&{\upphi}_{\text{SM}}^{J}\left(g_{1},\ldots,g_{K},g_{\bot,1},\ldots,g_{\bot,J};{\bm{x}},y,m_{1},\ldots,m_{J}\right)=\\\
&-\log\left(\frac{\exp\\{g_{y}({\bm{x}})\\}}{\sum_{y^{\prime}\in\mathcal{Y}^{\bot}}\exp\\{g_{y^{\prime}}({\bm{x}})\\}}\right)\\\
&-\sum_{j=1}^{J}\mathbb{I}[m_{j}=y]\
\log\left(\frac{\exp\\{g_{\bot,j}({\bm{x}})\\}}{\sum_{y^{\prime}\in\mathcal{Y}^{\bot}}\exp\\{g_{y^{\prime}}({\bm{x}})\\}}\right).\end{split}$
(8)

The intuition is that the first term maximizes the function $g_{k}$ associated
with the true label. The second term maximizes the rejection function
$g_{\bot,j}$ but only if the $j$th expert’s prediction is correct. At test
time, the classifier is obtained by taking the maximum over $k\in[1,K]$:
$\hat{y}=h({\bm{x}})=\operatorname*{arg\,max}_{k\in[1,K]}g_{k}({\bm{x}})$. The
rejection function is similarly formulated as

$r({\bm{x}})=\begin{cases}&0\ \ \text{ if }\
g_{h({\bm{x}})}>g_{\bot,j^{\prime}}\ \ \forall j^{\prime}\in[1,J]\\\
&\operatorname*{arg\,max}_{j\in[1,J]}g_{\bot,j}({\bm{x}})\ \ \text{
otherwise}.\end{cases}$

Our proof of the soundness of Equation 8 follows the same approach that
Mozannar and Sontag (2020) employed—specifically, a reduction to cost-
sensitive learning.

######  Theorem 3.1.

${\uppsi}_{\text{SM}}^{J}$ (Equation 8) is a convex (in $g$), calibrated
surrogate loss for the $0-1$ multi-expert learning to defer loss (Equation 6).

The complete proof is in Appendix A.1. The result guarantees that the
minimizers $g^{*}_{1},\ldots,g^{*}_{K},g^{*}_{\bot,1},\ldots,g^{*}_{\bot,J}$
correspond to the Bayes optimal classifier and rejector given in Equation 7.

###  3.2 One-vs-All Surrogate Loss

We next turn to the OvA surrogate loss. Let the label space
$\mathcal{Y}^{\bot}$ and the functions
$g_{1},\ldots,g_{K},g_{\bot,1},\ldots,g_{\bot,J}$ be defined just as above for
the softmax case. The OvA-based multi-expert L2D surrogate is then:

$\begin{split}&{\uppsi}_{\text{OvA}}^{J}(g_{1},\ldots,g_{K},g_{\bot,1},\ldots,g_{\bot,J};{\bm{x}},y,m_{1},\ldots,m_{J})=\\\
&{\upphi}[g_{y}({\bm{x}})]+\sum_{y^{\prime}\in\mathcal{Y},y^{\prime}\neq
y}{\upphi}[-g_{y^{\prime}}({\bm{x}})]\ \
+\sum_{j=1}^{J}{\upphi}[-g_{\bot,j}({\bm{x}})]\\\
&+\sum_{j=1}^{J}\mathbb{I}[m_{j}=y]\left({\upphi}[g_{\bot,j}({\bm{x}})]-{\upphi}[-g_{\bot,j}({\bm{x}})]\right)\end{split}$
(9)

where ${\upphi}:\\{\pm 1\\}\times\mathbb{R}\mapsto\mathbb{R}_{+}$ is again a
binary surrogate loss. The classifier and rejector are computed exactly as in
the softmax case.

To prove consistency, we cannot construct our consistency proof in the same
direct manner as Mozannar and Sontag (2020). Like Verma and Nalisnick (2022),
we proceed by the method of error correcting output codes (ECOC) (Dietterich
and Bakiri, 1995; Langford et al., 2005; Allwein et al., 2001; Ramaswamy et
al., 2014), a general technique for reducing multiclass problems to multiple
binary problems.

We prove the consistency of ${\uppsi}_{\text{OvA}}^{J}$ by way of the
following two results.

######  Theorem 3.2.

For a strictly proper binary composite loss ${\upphi}$ with a well-defined
continuous inverse link function ${\upgamma}^{-1}$,
${\uppsi}_{\text{OvA}}^{J}$ (Equation 9) is a calibrated surrogate for the
$0-1$ multi-expert learning to defer loss (Equation 6).

The complete proof is in Appendix A.2. Assuming minimizability (Steinwart,
2007)—i.e. that our hypothesis class is sufficiently large (all measurable
functions)—the calibration result from Theorem 3.2 implies consistency.

######  Corollary 3.3.

Assume that $g\in\mathcal{F}$, where $\mathcal{F}$ is the hypothesis class of
all measurable functions. Minimizability (Steinwart, 2007) is then satisfied
for ${\uppsi}_{\text{OvA}}^{J}$, and it follows that
${\uppsi}_{\text{OvA}}^{J}$ is a consistent surrogate for the $0-1$ multi-
expert learning to defer loss (Equation 6).

Thus, the minimizers of ${\uppsi}_{\text{OvA}}^{J}$ (over all measurable
functions) agree with the Bayes optimal classifier and rejector (Equation 7).

##  4 CONFIDENCE CALIBRATION

We next turn to the calibration (Dawid, 1982) properties of multi-expert L2D.
While training with a consistent loss should produce models that are well-
calibrated, previous work on the single expert setting found that the
underlying parameterizations can strongly influence calibration in practice.
Specifically, Verma and Nalisnick (2022) show that the softmax formulation’s
estimators can be unbounded, resulting in ‘probability’ estimates above one.
As for the calibration of the classifier, Verma and Nalisnick (2022) found
that there is to systemic issue and can be improved with standard post-hoc
techniques if necessary. Their findings also apply to the multi-expert
scenario and thus we consider only the rejector going forward.

We are particularly interested in the rejector’s ability to estimate
$\mathbb{P}({\textnormal{m}}_{j}={\textnormal{y}}|{\bm{x}})$, the conditional
probability that the $j$th expert is correct. If the L2D system says that
$\mathbb{P}({\textnormal{m}}_{j}={\textnormal{y}}|{\bm{x}}_{0})=0.7$, then the
$j$th expert should be correct $70\%$ of the time for inputs very similar to
${\bm{x}}_{0}$. This quantity is crucial not only for the system’s ability to
correctly defer but is also useful for interpretability and safety—to quantify
what the model thinks that the human knows.

We next define the relevant notion of calibration. For an estimator of expert
correctness $t({\bm{x}}):\mathcal{X}\mapsto(0,1)$, we call $t$ calibrated if,
for any confidence level $c\in(0,1)$, the actual proportion of times the
expert is correct is equal to $c$:

$\begin{split}\mathbb{P}({\textnormal{m}}={\textnormal{y}}\ |\
t({\bm{x}})=c)=c.\end{split}$ (10)

This statement should hold for all possible instances ${\bm{x}}$ with
confidence $c$. Since expert correctness is a binary classification problem,
distribution calibration, confidence calibration, and classwise calibration
all coincide (Vaicenavicius et al., 2019). We can measure the degree of
calibration using expected calibration error (ECE). In this case, the relevant
ECE is defined as

$\text{ECE}(t)=\mathbb{E}_{{\mathbf{x}}}|\mathbb{P}\left({\textnormal{m}}={\textnormal{y}}\
|\ t({\mathbf{x}})=c\right)-c|,$

where $\mathbb{E}_{{\mathbf{x}}}$ is usually approximated with samples.

###  4.1 Softmax Parameterization

For the softmax formulation, the estimator of the probability that the $j$th
expert is correct can be derived as follows; see Appendix B.1 for a complete
derivation. The Bayes optimal functions $g^{*}_{1},\ldots,g^{*}_{\bot,J}$ have
the following relationship with the underlying probability of expert
correctness:

$\begin{split}&\frac{\mathbb{P}({\textnormal{m}}_{j}={\textnormal{y}}|{\bm{x}})}{1+\sum_{j^{\prime}=1}^{J}\mathbb{P}({\textnormal{m}}_{j^{\prime}}={\textnormal{y}}|{\bm{x}})}=\\\
&\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\underbrace{\frac{\exp\\{g_{\bot,j}^{*}\\}}{\sum_{y^{\prime}\in\mathcal{Y}^{\bot}}\exp\\{g^{*}_{y^{\prime}}({\bm{x}})\\}}}_{t_{\bot,j}^{*}({\bm{x}})}.\end{split}$
(11)

Denote the RHS of Equation 11 as $t_{\bot,j}^{*}({\bm{x}})$. Since we have $J$
equations, one for each expert, we can uniquely solve for
$\mathbb{P}({\textnormal{m}}_{j}={\textnormal{y}}|{\bm{x}})$ as:

$\mathbb{P}({\textnormal{m}}_{j}={\textnormal{y}}|{\bm{x}})=\frac{t_{\bot,j}^{*}({\bm{x}})}{1-\sum_{j^{\prime}=1}^{J}t_{\bot,j^{\prime}}^{*}({\bm{x}})}.$
(12)

Equation 12 exhibits the same pathology as the single expert setting: it is
unbounded from above. For $t_{\bot,j}({\bm{x}})>0$, as
$\sum_{j^{\prime}=1}^{J}t_{\bot,j^{\prime}}({\bm{x}})$ approaches one, the
estimate of $\mathbb{P}({\textnormal{m}}_{j}={\textnormal{y}}|{\bm{x}})$ will
go to infinity. Moreover, the estimator for the $j$th expert depends on the
estimators for the other experts. Thus, if one $t_{\bot,j}({\bm{x}})$ is mis-
calibrated, this error will likely propagate to the other estimators.

###  4.2 One-vs-All Parameterization

For the OvA formulation, the probability that the expert is correct is
directly modeled by the $j$th deferral function. For the logistic binary loss
${\upphi}$, we have:

$\begin{split}\mathbb{P}({\textnormal{m}}_{j}={\textnormal{y}}|{\bm{x}})\ &=\
{\upphi}\left(g^{*}_{\bot,j}({\bm{x}})\right)\\\ &=\
\frac{1}{1+\exp\\{-g^{*}_{\bot,j}({\bm{x}})\\}}.\end{split}$ (13)

This estimator has the correct range of $(0,1)$ for any setting of
$g_{\bot,j}\in\mathbb{R}$. Moreover, there is no dependence across expert
deferral functions $g_{\bot,1},\ldots,g_{\bot,J}$. We expect these properties
to result in better calibration in practice.

##  5 ENSEMBLING EXPERTS WITH CONFORMAL INFERENCE

Multi-expert L2D, as defined above, operates by selecting just one expert upon
deferral. This approach is sensible if querying each expert results in an
independent expense (such as a consulting fee). However, in other settings,
the cost incurred by deferring may just be that of time and efficiency (i.e. a
lack of automation). In this case, the cost of querying additional experts
would be negligible; for example, we could send multiple experts a message in
parallel asking for their decision. Given the estimators of
$\mathbb{P}({\textnormal{m}}_{j}={\textnormal{y}}|{\bm{x}})$ presented in the
previous section, it is then natural to ask how we might ensemble experts
according to these estimates of correctness. Below we present a methodology
based on conformal inference for obtaining dynamic, minimal ensembles of
experts.

##### Conformal Inference

Conformal inference (CI) (Shafer and Vovk, 2008) constructs a confidence
interval (or set) for predictive inference. In the traditional multi-class
classification setting, given a new observation ${\bm{x}}_{n+1}$, we wish to
determine the correct associated label ${\textnormal{y}}_{n+1}=y_{n+1}^{*}$,
where $y_{n+1}^{*}$ denotes the true class label. CI allows us to construct a
distribution-free confidence set $\mathsf{C}({\bm{x}}_{n+1})$ that will cover
the true label with _marginal_ probability $1-\alpha$:

$\mathbb{P}\left(y_{n+1}^{*}\notin\mathsf{C}({\bm{x}}_{n+1})\right)\ \leq\
\alpha\ \ \ \ \forall\ \mathbb{P}\in\mathfrak{P}$

where $\mathfrak{P}$ represents the space of all distributions—hence the
qualifier ‘distribution-free.’ Denote the test statistic as
$S({\bm{x}},y;\mathcal{D})$. It is known as a non-conformity function: a
higher value of $S$ means that $({\bm{x}},y)$ is less conforming to the
distribution represented (empirically) by $\mathcal{D}$. Despite this
guarantee, CI is only as good as its test statistic in practice. For instance,
the marginal coverage is naively satisfied if we construct the set randomly by
setting $\mathsf{C}({\bm{x}})=\mathcal{Y}$ with probability $1-\alpha$ and
returning the empty set otherwise. CI is implemented by calculating the non-
conformity function on a validation set and computing the empirical $1-\alpha$
quantile (with a finite sample correction). At test time, elements are added
to the set until the non-conformity function passes the previously-computed
quantile.

##### Conformal Sets of Experts

We propose applying CI to perform uncertainty quantification for the experts.
Thus, here, $\mathsf{C}({\bm{x}})$ represents a set of experts. Firstly, we
assume there is a best expert: for a new observation ${\bm{x}}_{n+1}$, let
$j_{n+1}^{*}$ denote the best expert such that

$\mathbb{P}({\textnormal{m}}_{j_{n+1}^{*}}={\textnormal{y}}|{\bm{x}}_{n+1})>\mathbb{P}({\textnormal{m}}_{e}={\textnormal{y}}|{\bm{x}}_{n+1})\
\forall e\neq j^{*}_{n+1}.$

We would then like to construct a set such that $j_{n+1}^{*}$ is covered with
marginal probability $1-\alpha$:

$\mathbb{P}\left(j_{n+1}^{*}\notin\mathsf{C}({\bm{x}}_{n+1})\right)\ \leq\
\alpha\ \ \ \ \forall\ \mathbb{P}\in\mathfrak{P}$

where $\mathsf{C}({\bm{x}}_{n+1})$ again denotes the conformal set and
$\mathfrak{P}$ is the same as above. Unlike in most applications of CI, we can
use the procedure to form an ensemble by aggregating the predictions of all
experts in the set. Moreover, this set will have a dynamic size that changes
with ${\bm{x}}$, ensuring our ensemble makes efficient use of expert queries.

##### Naive Statistic

We start by adapting a score function from multi-class classification. Let
$s_{j}({\bm{x}})$ denote the estimator that the $j$th expert is correct. For
the softmax case,
$s_{j}({\bm{x}})=t_{\bot,j}({\bm{x}})/(1-\sum_{j^{\prime}}t_{\bot,j^{\prime}}({\bm{x}}))$
(Equation 12), and for OvA, $s_{j}({\bm{x}})={\upphi}(g_{\bot,j}({\bm{x}}))$
(Equation 13). Let $\pi_{1},\ldots,\pi_{J}$ denote the indices for a
descending ordering of the estimators $s_{j}({\bm{x}})$, i.e. $s_{\pi_{1}}$ is
the expert who has the best chance of being correct (according to the
rejector). The resulting non-conformity function and test statistic are:

$S\left({\bm{x}},y,m_{1},\ldots,m_{J};\mathcal{D}\right)=\sum_{e=1}^{E}s_{\pi_{e}}({\bm{x}})$
(14)

where $\pi_{E}$ is the index of the expert who has the lowest score
$s_{\pi_{E}}$ of all correct experts ($m=y$). This expression means that we
will keep adding the correctness scores $s$ in descending order until we
include all experts who correctly predict the given instance. Hence $E=J$ only
when all experts are correct and $E<J$ otherwise.

##### Regularized Statistic

A problem with the statistic above is that multiple experts can be correct,
resulting in noise that obscures the identity of the best expert. In the
experiments, we show that this ‘naive’ statistic is not robust to noise,
resulting in inflated set sizes (which are sometimes vacuous). Similar
problems are discussed by Angelopoulos et al. (2020). To address this issue,
we employ conformal risk control (Angelopoulos et al., 2022) to directly
control the false negative rate. We create regularized prediction sets as
follows:

$C_{\lambda}({\bm{x}})=\left\\{j:s_{j}({\bm{x}})+\beta\left(s_{j}({\bm{x}})-k_{\text{reg}}\right)>1-\lambda\right\\}$
(15)

where $\beta$ and $k_{\text{reg}}$ are the parameters of the regularization.
We describe the procedure to choose these hyperparameters in Appendix C. The
general idea is to choose $k_{\text{reg}}$ so that confidences lower than this
threshold can happen with probability at most $\alpha$. We choose $\beta$ to
optimize for prediction set size.

Figure 1: System Accuracy with Increasing Experts. The plots above report the
classification accuracy of the overall system as the number of experts
increases from $2$ to $16$. The softmax (orange) and one-vs-all (blue) losses
are compared on three tasks: galaxy categorization (Galaxy-Zoo), skin lesion
diagnosis (HAM10000), and hate speech detection (Hate Speech). The results
show the losses offer roughly equivalent performance.

##  6 RELATED WORK

##### Learning to Defer

Learning to defer (L2D) was proposed in its modern form by Madras et al.
(2018). Yet classifiers with an option to reject or abstain date back at least
to Chow (1957). There have been two primary approaches to making the rejection
decision: confidence-based (Bartlett and Wegkamp, 2008; Yuan and Wegkamp,
2010; Jiang et al., 2018; Grandvalet et al., 2009; Ramaswamy et al., 2018; Ni
et al., 2019) and model-based (Cortes et al., 2016a, b). Our work, due to it
using a parameterized rejector, falls into the latter. The theoretical
properties of the classifier-rejector approach have been well-studied for
binary classification (Cortes et al., 2016a, b). The theory for multi-class
classification was first developed by Ni et al. (2019) and Charoenphakdee et
al. (2021). Mozannar and Sontag (2020) then built upon this work to establish
the first consistent surrogate loss for multi-class L2D. Previous L2D
extensions did not come with consistency guarantees (Raghu et al., 2019;
Wilder et al., 2020; Pradier et al., 2021; Okati et al., 2021; Liu et al.,
2022). Verma and Nalisnick (2022) proposed the second provably consistent
surrogate for multi-class L2D based on a one-vs-all formulation. Charusaie et
al. (2022) further studied the L2D optimization problem, proving results for
complementarity and active learning. Our work extends Mozannar and Sontag
(2020)’s and Verma and Nalisnick (2022)’s results to the mult-expert
setting—for which no one has yet to propose a consistent surrogate loss.

##### Calibration in L2D

Verma and Nalisnick (2022) motivate their OvA surrogate from the standpoint of
calibration and thus is the only other work that studies the confidence
calibration of L2D systems. We extend their work to the multi-expert setting.
Calibration has received much attention of late in the wider machine learning
literature (Guo et al., 2017; Kull et al., 2019; Gupta and Ramdas, 2022;
Vaicenavicius et al., 2019). The dominant methodology is to apply post-hoc
calibration: fitting additional parameters on validation data to re-calibrate
the formerly mis-calibrated model. These methods could potentially be applied
here—such as, by adding a temperature parameter to the per-expert terms in the
OvA loss—but we are primarily interested in the native, ‘out-of-the-box’
calibration properties of the losses.

##### Multi-Expert Models

There have been several works that use models to improve the decision making
of multiple experts (Benz and Rodriguez, 2022; Straitouri et al., 2022) and to
fuse decisions from models and humans (Keswani et al., 2021; Kerrigan et al.,
2021). Keswani et al. (2021)’s work is perhaps the most related to our work.
However, their model does not conform to the L2D paradigm as the classifier
and _all_ experts are queried to make _every_ prediction.

##### Conformal Inference for Human-AI Collaboration

We know of two works that have used CI for some form of human-AI
collaboration. Straitouri et al. (2022) apply CI to a classifier and then pass
the prediction set to a human to make the final decision. Babbar et al. (2022)
study a similar work flow (apply CI then pass to a human) and also propose
applying CI only to non-deferred samples, which results in smaller set sizes.
No previous work has applied CI to obtain sets of experts.

##  7 EXPERIMENTS

Our experimental setup closely follows that of Verma and Nalisnick (2022)—but
extended to multiple experts. For all runs, we report the mean and standard
error across $3$ random seeds. We perform three types of experiments. In the
first, we check the system accuracy of the derived consistent surrogate losses
in three consequential tasks (Subsection 7.1): galaxy classification, skin
lesion diagnosis, and hate speech detection. We find that both surrogates
achieve similar accuracies. Secondly, we investigate the confidence
calibration properties of the surrogates (Subsection 7.2). As hypothesized,
the OvA loss results in less calibration error in both of our simulations.
Lastly, we investigate the efficacy of our conformal ensembling procedure
(Subsection 7.3). For the naive statistic, the OvA loss’ superior calibration
results in appropriately smaller sets. For the regularized statistic, both
losses perform equally well.

(a) ECE as Experts Increase

(b) Acc. as Experts Increase

(c) ECE for Improving Experts

(d) ECE of Random Expert

Figure 2: Confidence Calibration Simulations. The figures above report the
results for confidence calibration simulations performed on CIFAR-10. The
first and second subfigures show calibration error (a) and system accuracy (b)
as the number of experts increases from $4$ to $20$. We see that the OvA
formulation (blue) has better calibration across all runs, but this translate
to better accuracy only for $16$ or fewer experts. The third and forth
subfigures show calibration error as experts’ abilities increase (a) and when
one expert is kept at random chance (b). OvA (blue) shows better calibration
in both metrics.

###  7.1 Overall System Accuracy

##### Data Sets

We report the overall system accuracy for three real-world data sets: HAM10000
(Tschandl et al., 2018) for skin lesions diagnosis, Galaxy-Zoo (Bamford et
al., 2009) for galaxy classification, and HateSpeech (Davidson et al., 2017)
for hate speech detection. The train-validation-test split is
$60\%-20\%-20\%$. Following Verma and Nalisnick (2022), we down-sample Galaxy-
Zoo to $10,000$ instances.

##### Simulated Experts

We simulate $16$ unique experts. We derive HAM10000’s expert predictions from
the predictions of an $8$-layer MLPMixer (Tolstikhin et al., 2021), which has
access to additional metadata such as age, gender, and diagnosis type. For
Galaxy-Zoo and HateSpeech, we sample predictions from the provided human
annotations. We construct additional experts using various perturbations of
these predictions. One is an expert who flips the given prediction with some
probability. Another has expertise across a subset of the classes. Another has
a uniformly random chance across all classes.

##### Models

We use a $34$-layer residual network (ResNet34) and a $50$-layer residual
network (ResNet50) as base models for HAM10000 and Galaxy-Zoo respectively.
For HateSpeech, we use a $100$ dimensional fasttext (Joulin et al., 2016)
representation of each input and a ConvNet (Kim, 2014) as the base model. We
refer the reader to Verma and Nalisnick (2022) for more details on the
training and hyperparameter selection, as we follow their setup.

##### Experimental Setup

We train the systems with an increasing number of experts, ranging from $2$ to
$16$. For each run, we enlarge the expert pool by randomly sampling an expert
from our set of simulated experts, and this process is repeated $3$ times with
different random seeds. We keep the base model fixed across these runs except
for the additional output dimension. Since there is no guarantee that an
additional expert is superior to the ones already present, the system accuracy
should not necessarily increase as the expert pool expands.

##### Results

Figure 1 reports the mean and standard error of the system accuracy as the
number of experts increases. The center of the circle denotes the mean, and
its radius is proportional to the standard error. The surrogate loss functions
appear to be comparable. OvA (blue) demonstrates slightly better performance
on Galaxy-Zoo when there are fewer experts ($2$ and $4$).

(a) Set Size, No Noise

(b) Set Size, with Noise

(c) CIFAR-10: Sys. Accuracy

(d) HAM10000: Sys. Accuracy

Figure 3: Conformal Sets of Experts. The figures above report our analysis of
the two statistics proposed in Section 5. Subfigures (a) and (b) show the
ability of the two statistics to select the correct number of experts as the
number of oracle experts increases—so optimal performance is the $y=x$ line.
Subfigure (a) reports the no-noise setting (so experts are either perfectly
incorrect or correct), and we see that the naive statistic (solid lines)
overestimates the set size. The problem is even worse in Subfigure (b).
However, the regularized statistic (dashed lines) is able to do well in both
cases. Subfigures (c) and (d) show how ensembling the set affects system
accuracy. The conformal approach is able to out-perform a fixed size of $5$
experts for a small number of oracles and is equivalent at higher numbers.

###  7.2 Confidence Calibration

In Section 4, we found that the two surrogates have very different estimators
of
$\mathbb{P}\left({\textnormal{m}}_{j}={\textnormal{y}}_{i}|{\mathbf{x}}_{i}\right)$,
the probability that the $j$th expert is correct. We now test if these
theoretical differences have consequences for practice. To ensure ECE is well-
defined for the softmax loss, we cap any confidences greater than $1$ at $1$.
We perform all simulations using the standard splits of CIFAR-10 (Krizhevsky,
2009) and a $28$-layer wide residual network (Zagoruyko and Komodakis, 2016),
following Verma and Nalisnick (2022).

##### Simulation #1: Increasing Experts

We perform a simulation to see how the methods perform under an increasing
number of experts. We generate a synthetic expert with a correctness
probability of $70\%$ over the first five classes and random across all other
classes. We then replicate that expert and add it to the expert pool, ranging
from $4$ to $20$ total experts. Figure 2(a) reports the average ECE across
experts as the pool increases. The OvA method (blue) is roughly stable at
about $2\%$ ECE as experts are added. The softmax method (orange) has roughly
double the ECE ($\sim 4.5\%$). In Figure 2(b), we report the overall system
accuracy to see if these calibration differences have an effect. We see some
positive effects, with OvA (blue) having a better accuracy for $12$ and fewer
experts. However, the softmax (orange) has the best accuracy at $20$ experts,
despite its calibration still being worse.

##### Simulation #2: Expert Dependence

We next perform a simulation to see how calibration error can propagate across
the estimators. We simulate four experts with one always being random and the
other three having a probability of correctness that increases from $20\%$ to
$95\%$ on the first five classes (random for others). We hypothesize that the
softmax’s ECE for the random expert will _increase_ when the probability of
correctness for the other three experts increases due to the tied
parameterization. Figures 2(c) and 2(d) report the results, with the former
reporting average ECE and the latter the ECE of just the random expert.
Firstly, from Figure 2(c), we see that again the OvA method is better
calibrated across all experimental settings. Then from Figure 2(d), we see
that our hypothesis is confirmed: OvA (blue) is able to model the random
expert well no matter the other experts’ abilities, but the softmax (orange)
is not. The softmax’s ECE increases almost in-step with the expert
correctness, except for some cancellation effect at $80\%$. This is clearly an
undesirable behavior from the standpoint of safety since any ECE above zero
means that the system is reporting that the expert is better than random and
thus misleading the user.

###  7.3 Conformal Ensembles

Lastly, we study our proposal of using CI to ensemble multiple experts. We
first analyze the two proposed statistics, demonstrating the regularized
version’s superior ability to recover the experts who are oracles. We then
report the downstream effect on the overall system accuracy, comparing
performance to that of a fixed-size ensemble of experts.

##### Experts and Setup

We experiment with two settings on CIFAR-10, each with $10$ total experts. In
the first (no noise), we synthesize experts such that they are an oracle on an
increasing subset of the classes and guaranteed to be wrong on the classes not
in that set. In the second (with noise), the experts are oracles in the same
way but now have a (uniformly) random chance of being correct for the non-
oracle classes. The theory of CI guarantees that the sets marginally cover all
oracle experts. Yet, ideally, we wish the sets to contain _only_ the experts
who are oracles. We use $\alpha=0.1$ in all experiments.

##### Expert Identification

The CI results are reported in Figures 3(a) and 3(b). The number of oracles is
on the y-axis, and the average set size is on the x-axis. Optimal performance
would be the $y=x$ line. The results for the no-noise setting are reported in
Figure 3(a). The naive statistic (solid lines) considerably inflates the set
size for both softmax and OvA. Yet OvA is much closer to $y=x$, which suggests
superior calibration leads to better CI. The performance of the regularized
statistic is shown by the dashed lines. Both softmax and OvA perform nearly
perfectly. Figure 3(b) reports the wish-noise setting, and we find that the
naive statistic performs terribly for both losses. The regularized statistic,
on the other hand, performs well for both softmax and OvA. Softmax
demonstrates slight superiority for $2-5$ oracles.

##### Overall System Accuracy

We next investigate using the conformal set as an ensembling strategy. Upon
deferral, we use majority voting across the set to generate the final
prediction. We compare this with the baseline of using a fixed ensemble
size—specifically, the top five ranked experts. We show the results for
CIFAR-10 in Figure 3(c). The crucial settings are for one and two oracles
since using a fixed ensemble size is guaranteed to fail—as is confirmed by the
plot ($<80\%$ accuracy). We see that the conformal ensembles are clearly
superior here, achieving around $90\%$ accuracy. For three or more oracles,
both methods have equal performance. This is expected since only three oracles
are needed to form a correct majority. Figure 3(d) reports system accuracy for
HAM10000 using a sub-sample of the same experts described in Section 7.1. All
methods perform within standard error of one another. Yet we emphasize that
CI’s adaptivity is still desirable so that the best experts can be identified
with transparency and queried efficiently.

##  8 CONCLUSIONS

We have extended the L2D framework to support multiple experts. We proposed
two optimization objectives and proved that they are both consistent. We also
studied their potential to be confidence calibrated, showing that the softmax-
based objective can result in mis-calibrated models in practice. Lastly, we
considered a principled procedure for selecting _minimal_ sets of experts to
ensemble. For future work, we aim to improve the data efficiency of our method
by extending the active learning results of Charusaie et al. (2022) to the
multi-expert setting.

#### Acknowledgements

This publication is part of the project Continual Learning under Human
Guidance (VI.Veni.212.203), which is financed by the Dutch Research Council
(NWO). This work was carried out on the Dutch national e-infrastructure with
the support of the SURF Cooperative.

## References

  * Akata et al. (2020) Zeynep Akata, Dan Balliet, Maarten De Rijke, Frank Dignum, Virginia Dignum, Guszti Eiben, Antske Fokkens, Davide Grossi, Koen Hindriks, Holger Hoos, et al.  A research agenda for hybrid intelligence: augmenting human intellect with collaborative, adaptive, responsible, and explainable artificial intelligence.  _Computer_ , 53(08):18–28, 2020.
  * Allwein et al. (2001) Erin L. Allwein, Robert E. Schapire, and Yoram Singer.  Reducing multiclass to binary: A unifying approach for margin classifiers.  _Journal of Machine Learning Research_ , 1:113–141, 2001\.
  * Angelopoulos et al. (2020) Anastasios Angelopoulos, Stephen Bates, Jitendra Malik, and Michael I Jordan.  Uncertainty sets for image classifiers using conformal prediction.  _arXiv preprint arXiv:2009.14193_ , 2020.
  * Angelopoulos et al. (2022) Anastasios N Angelopoulos, Stephen Bates, Adam Fisch, Lihua Lei, and Tal Schuster.  Conformal risk control.  _arXiv preprint arXiv:2208.02814_ , 2022.
  * Babbar et al. (2022) Varun Babbar, Umang Bhatt, and Adrian Weller.  On the utility of prediction sets in human-ai teams.  _International Joint Conference on Artificial Intelligence (IJCAI)_ , 2022.
  * Bamford et al. (2009) Steven P. Bamford, Robert C. Nichol, Ivan K. Baldry, Kate Land, Chris J. Lintott, Kevin Schawinski, Anže Slosar, Alexander S. Szalay, Daniel Thomas, Mehri Torki, Dan Andreescu, Edward M. Edmondson, Christopher J. Miller, Phil Murray, M. Jordan Raddick, and Jan Vandenberg.  Galaxy Zoo: the dependence of morphology and colour on environment*.  _Monthly Notices of the Royal Astronomical Society_ , 393(4):1324–1352, 02 2009.  ISSN 0035-8711.  doi: 10.1111/j.1365-2966.2008.14252.x.
  * Bartlett and Wegkamp (2008) Peter L. Bartlett and Marten H. Wegkamp.  Classification with a reject option using a hinge loss.  _Journal of Machine Learning Research_ , 9:1823–1840, jun 2008.  ISSN 1532-4435.
  * Benz and Rodriguez (2022) Nina Corvelo Benz and Manuel Gomez Rodriguez.  Counterfactual inference of second opinions.  _arXiv preprint arXiv:2203.08653_ , 2022.
  * Charoenphakdee et al. (2021) Nontawat Charoenphakdee, Zhenghang Cui, Yivan Zhang, and Masashi Sugiyama.  Classification with rejection based on cost-sensitive classification.  In _International Conference on Machine Learning_ , 2021.
  * Charusaie et al. (2022) Mohammad-Amin Charusaie, Hussein Mozannar, David Sontag, and Samira Samadi.  Sample efficient learning of predictors that complement humans.  In _International Conference on Machine Learning_ , pages 2972–3005. PMLR, 2022.
  * Chow (1957) C. K. Chow.  An optimum character recognition system using decision functions.  _IRE Trans. Electron. Comput._ , 6:247–254, 1957.
  * Cortes et al. (2016a) Corinna Cortes, Giulia DeSalvo, and Mehryar Mohri.  Learning with rejection.  In _International Conference on Algorithmic Learning Theory_ , 2016a.
  * Cortes et al. (2016b) Corinna Cortes, Giulia DeSalvo, and Mehryar Mohri.  Boosting with abstention.  In _Advances in Neural Information Processing Systems_ , 2016b.
  * Davidson et al. (2017) Thomas Davidson, Dana Warmsley, Michael W. Macy, and Ingmar Weber.  Automated hate speech detection and the problem of offensive language.  In _International AAAI Conference on Web and Social Media_ , 2017\.
  * Dawid (1982) A Philip Dawid.  The well-calibrated Bayesian.  _Journal of the American Statistical Association_ , 77(379):605–610, 1982.
  * Dellermann et al. (2019) Dominik Dellermann, Philipp Ebel, Matthias Söllner, and Jan Marco Leimeister.  Hybrid intelligence.  _Business & Information Systems Engineering_, 61(5):637–643, 2019.
  * Dietterich and Bakiri (1995) Thomas G. Dietterich and Ghulum Bakiri.  Solving multiclass learning problems via error-correcting output codes.  _Journal of Artificial Intelligence Research_ , 2(1):263–286, jan 1995.  ISSN 1076-9757.
  * Fay et al. (2006) Doris Fay, Carol Borrill, Ziv Amir, Robert Haward, and Michael A West.  Getting the most out of multidisciplinary teams: A multi-sample study of team innovation in health care.  _Journal of Occupational and Organizational Psychology_ , 79(4):553–567, 2006.
  * Grandvalet et al. (2009) Yves Grandvalet, Alain Rakotomamonjy, Joseph Keshet, and Stéphane Canu.  Support vector machines with a reject option.  In _Advances in Neural Information Processing Systems_ , 2009.
  * Guo et al. (2017) Chuan Guo, Geoff Pleiss, Yu Sun, and Kilian Q. Weinberger.  On calibration of modern neural networks.  In _International Conference on Machine Learning_ , 2017.
  * Gupta and Ramdas (2022) Chirag Gupta and Aaditya Ramdas.  Top-label calibration and multiclass-to-binary reductions.  In _International Conference on Learning Representations_ , 2022.
  * Jiang et al. (2018) Heinrich Jiang, Been Kim, Melody Y. Guan, and Maya Gupta.  To trust or not to trust a classifier.  In _Advances in Neural Information Processing Systems_ , 2018.
  * Joulin et al. (2016) Armand Joulin, Edouard Grave, Piotr Bojanowski, Matthijs Douze, Hervé Jégou, and Tomás Mikolov.  Fasttext.zip: Compressing text classification models.  _arXiv preprint arXiv:1612.03651_ , 2016.
  * Kamar (2016) Ece Kamar.  Directions in hybrid intelligence: Complementing ai systems with human intelligence.  In _International Joint Conference on Artificial Intelligence_ , pages 4070–4073, 2016.
  * Kerrigan et al. (2021) Gavin Kerrigan, Padhraic Smyth, and Mark Steyvers.  Combining human predictions with model probabilities via confusion matrices and calibration.  _Advances in Neural Information Processing Systems_ , 2021.
  * Keswani et al. (2021) Vijay Keswani, Matthew Lease, and Krishnaram Kenthapadi.  Towards unbiased and accurate deferral to multiple experts.  In _Proceedings of the 2021 AAAI/ACM Conference on AI, Ethics, and Society_ , pages 154–165, 2021.
  * Kim (2014) Yoon Kim.  Convolutional neural networks for sentence classification.  In _Conference on Empirical Methods in Natural Language Processing_ , 2014.
  * Krizhevsky (2009) Alex Krizhevsky.  Learning multiple layers of features from tiny images.  _Technical Report_ , 2009.
  * 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 multiclass probabilities with dirichlet calibration.  In _Advances in Neural Information Processing Systems_ , 2019.
  * Langford et al. (2005) John Langford, Tti-Chicago, Jl@hunch Net, and Alina Beygelzimer.  Sensitive error correcting output codes.  In _Conference on Learning Theory_ , 2005.
  * Liu et al. (2022) Jessie Liu, Blanca Gallego, and Sebastiano Barbieri.  Incorporating Uncertainty in Learning to Defer Algorithms for Safe Computer-Aided Diagnosis.  _Scientific reports_ , 12(1):1–9, 2022.
  * Madras et al. (2018) David Madras, Toniann Pitassi, and Richard Zemel.  Predict responsibly: Improving fairness and accuracy by learning to defer.  In _Advances in Neural Information Processing Systems_ , 2018.
  * Mozannar and Sontag (2020) Hussein Mozannar and David A. Sontag.  Consistent estimators for learning to defer to an expert.  In _International Conference on Machine Learning_ , 2020.
  * Ni et al. (2019) Chenri Ni, Nontawat Charoenphakdee, Junya Honda, and Masashi Sugiyama.  On the calibration of multiclass classification with rejection.  In _NeurIPS_ , 2019.
  * Okati et al. (2021) Nastaran Okati, Abir De, and Manuel Gomez-Rodriguez.  Differentiable learning under triage.  In _Advances in Neural Information Processing Systems_ , 2021.
  * Pradier et al. (2021) Melanie F Pradier, Javier Zazo, Sonali Parbhoo, Roy H Perlis, Maurizio Zazzi, and Finale Doshi-Velez.  Preferential mixture-of-experts: Interpretable models that rely on human expertise as much as possible.  _AMIA Summits on Translational Science Proceedings_ , 2021:525, 2021.
  * Raghu et al. (2019) Maithra Raghu, Katy Blumer, Greg Corrado, Jon Kleinberg, Ziad Obermeyer, and Sendhil Mullainathan.  The algorithmic automation problem: Prediction, triage, and human effort.  _arXiv preprint arXiv:1903.12220_ , 2019.
  * Ramaswamy et al. (2018) H. G. Ramaswamy, Ambuj Tewari, and Shivani Agarwal.  Consistent algorithms for multiclass classification with an abstain option.  _Electronic Journal of Statistics_ , 12:530–554, 2018.
  * Ramaswamy et al. (2014) Harish G. Ramaswamy, Balaji Srinivasan Babu, Shivani Agarwal, and Robert C. Williamson.  On the consistency of output code based learning algorithms for multiclass learning problems.  In Maria Florina Balcan, Vitaly Feldman, and Csaba Szepesvári, editors, _Proceedings of The 27th Conference on Learning Theory_ , volume 35 of _Proceedings of Machine Learning Research_ , pages 885–902, Barcelona, Spain, 13–15 Jun 2014. PMLR.  URL https://proceedings.mlr.press/v35/ramaswamy14.html.
  * Shafer and Vovk (2008) Glenn Shafer and Vladimir Vovk.  A tutorial on conformal prediction.  _Journal of Machine Learning Research_ , 9(3), 2008.
  * Steinwart (2007) Ingo Steinwart.  How to compare different loss functions and their risks.  _Constructive Approximation_ , 26:225–287, 2007.
  * Straitouri et al. (2022) Eleni Straitouri, Lequn Wang, Nastaran Okati, and Manuel Gomez Rodriguez.  Provably improving expert predictions with conformal prediction.  _arXiv preprint arXiv:2201.12006_ , 2022.
  * Tolstikhin et al. (2021) Ilya O. Tolstikhin, Neil Houlsby, Alexander Kolesnikov, Lucas Beyer, Xiaohua Zhai, Thomas Unterthiner, Jessica Yung, Daniel Keysers, Jakob Uszkoreit, Mario Lucic, and Alexey Dosovitskiy.  Mlp-mixer: An all-mlp architecture for vision.  In _Advances in Neural Information Processing Systems_ , 2021.
  * Tschandl et al. (2018) Philipp Tschandl, Cliff Rosendahl, and Harald Kittler.  The ham10000 dataset, a large collection of multi-source dermatoscopic images of common pigmented skin lesions.  _Scientific data_ , 5(1):1–9, 2018.
  * Vaicenavicius et al. (2019) Juozas Vaicenavicius, David Widmann, Carl Andersson, Fredrik Lindsten, Jacob Roll, and Thomas Schön.  Evaluating model calibration in classification.  In _Conference on Artificial Intelligence and Statistics_ , 2019.
  * Verma and Nalisnick (2022) Rajeev Verma and Eric Nalisnick.  Calibrated learning to defer with one-vs-all classifiers.  In Kamalika Chaudhuri, Stefanie Jegelka, Le Song, Csaba Szepesvari, Gang Niu, and Sivan Sabato, editors, _Proceedings of the 39th International Conference on Machine Learning_ , volume 162 of _Proceedings of Machine Learning Research_ , pages 22184–22202. PMLR, 17–23 Jul 2022.  URL https://proceedings.mlr.press/v162/verma22c.html.
  * Wilder et al. (2020) Bryan Wilder, Eric Horvitz, and Ece Kamar.  Learning to complement humans.  In _IJCAI_ , 2020.
  * Yuan and Wegkamp (2010) Ming Yuan and Marten Wegkamp.  Classification methods with reject option based on convex risk minimization.  _Journal of Machine Learning Research_ , 11(5):111–130, 2010.
  * Zagoruyko and Komodakis (2016) Sergey Zagoruyko and Nikos Komodakis.  Wide residual networks.  _arXiv preprint arXiv:1605.07146_ , 2016.

##  Appendix A Consistency Proofs

In this section, we provide proofs for the main results in the paper. First we
prove the Bayes optimal rule for L2D to multiple experts, then we prove that
the surrogate losses proposed in the paper are consistent. We continue the
notation from the main paper. For simplicity, we do not worry about measure-
theoretic considerations and assume that appropriate conditions hold that
allow us to interchange summations and integrals, for example.

Proof of the Bayes Rule for Learning to Defer with Multiple Experts

We have $J$ experts and a classifier, where the system either allows the
classifier to make the final prediction or defers to one of the $J$ experts.
When the classifier makes the prediction, the system incurs loss
$\ell_{\text{clf}}(\hat{{\textnormal{y}}},{\textnormal{y}})$ where
$\hat{{\textnormal{y}}}=h({\bm{x}})$. When the system defers to the
$j^{\text{th}}$ expert, it incurs a loss
$\ell_{\text{exp}}({\textnormal{m}}_{j},{\textnormal{y}})$. In what follows,
we frame the learning to defer problem as a general classification problem,
and aim to find a function
$g:\mathcal{X}\rightarrow\hat{\mathcal{Y}}:=\mathcal{Y}\cup\\{\bot_{1},\bot_{2},\ldots,\bot_{J}\\}$
and $|\hat{\mathcal{Y}}|=K+J$ with the minimum expected loss (also known as
risk). We consider $g$ as modeling a probabilistic decision rule
$\mathbf{\delta}(\hat{y}|{\bm{x}}):=[\delta(\hat{y}=1|{\bm{x}}),\delta(\hat{y}=2|{\bm{x}}),\ldots,\delta(\hat{y}=K+J|{\bm{x}})]$
where $\delta(\hat{y}=i|{\bm{x}})$ denotes the confidence in making the
$i^{\text{th}}$ decision for ${\bm{x}}\sim{\mathbf{x}}$. We write risk as
follows:

$\displaystyle\begin{split}\mathcal{R}_{\mathcal{D}}[\mathbf{\delta}\left(\hat{y}|{\mathbf{x}}\right)]={{\sum}}_{i=1}^{\text{K}}{{\sum}}_{j=1}^{\text{K}+\text{J}}\bigintsss_{{\mathbf{x}}}\delta(\hat{y}=j|{\mathbf{x}})\ell\left(\hat{{\textnormal{y}}}=j,{\textnormal{y}}=i\right)\mathbb{P}\left({\textnormal{y}}=i\right)\mathbb{P}\left({\mathbf{x}}|{\textnormal{y}}=i\right)d{\mathbf{x}},\end{split}$
(16)

where $\ell:(\hat{{\textnormal{y}}},{\textnormal{y}})\mapsto\mathbb{R}_{+}$ is
a general loss function, $i$ runs over the input label space, and $j$ runs
over the output prediction space (classifier and all the experts). We further
expand the risk in Equation A based on the definition of the loss function in
learning to defer

$\displaystyle\begin{split}\mathcal{R}_{\mathcal{D}}[\mathbf{\delta}\left(\hat{y}|{\mathbf{x}}\right)]&={{\sum}}_{i=1}^{\text{K}}\bigintss_{{\mathbf{x}}}\Biggr{(}{{\sum}}_{j=1}^{\text{K}}\delta(\hat{y}_{j}|{\mathbf{x}})\ell_{\text{clf}}\left(j,i\right)\\\
&+{{\sum}}_{j=\text{K}+1}^{\text{K}+\text{J}}\left({{\sum}}_{m=1}^{\text{K}}\delta(\hat{y}_{j}|{\mathbf{x}})\ell_{exp}\left(m_{j},{\textnormal{y}})\mathbb{P}(m_{j}|{\mathbf{x}},{\textnormal{y}}=i)\right)\right)\Biggr{)}\mathbb{P}\left({\textnormal{y}}=i\right)\mathbb{P}\left({\mathbf{x}}|{\textnormal{y}}=i\right)d{\mathbf{x}},\end{split}$

where we have used shorthand $\delta\left(\hat{y}_{j}|{\mathbf{x}}\right)$ to
denote $\delta\left(\hat{y}=j|{\mathbf{x}}\right)$, and
$\mathbb{P}\left(m_{j}|{\mathbf{x}},{\textnormal{y}}=i\right)=\mathbb{P}\left({\textnormal{m}}_{j}=m|{\mathbf{x}},{\textnormal{y}}=i\right)$.
Next, we denote:

$\displaystyle\begin{split}w_{i,j}&=\ell_{\text{clf}}(j,i)\\\
w_{i,\bot_{j}}&={{\sum}}_{m=1}^{\text{K}}\delta(\hat{y}_{j}|{\mathbf{x}})\ell_{exp}\left(m_{j},{\textnormal{y}})\mathbb{P}(m_{j}|{\mathbf{x}},{\textnormal{y}}=i)\right)\end{split}$

Thus,
$\mathcal{R}_{\mathcal{D}}[\mathbf{\delta}\left(\hat{y}|{\mathbf{x}}\right)]$
can be written as:

$\displaystyle\mathcal{R}_{\mathcal{D}}[\mathbf{\delta}\left(\hat{y}|{\mathbf{x}}\right)]$
$\displaystyle={{\sum}}_{i=1}^{\text{K}}\bigintsss_{{\mathbf{x}}}\left({{\sum}}_{j=1}^{\text{K}}\delta(\hat{y}_{j}|{\mathbf{x}})w_{i,j}+{{\sum}}_{j=\text{K}+1}^{\text{K}+\text{J}}\delta(\hat{y}_{j}|{\mathbf{x}})w_{i,\bot_{j}}\right)\mathbb{P}\left({\textnormal{y}}=i\right)\mathbb{P}\left({\mathbf{x}}|{\textnormal{y}}=i\right)d{\mathbf{x}}.$

Denote $w^{*}_{i,\bot}:=\min_{j\in[J]}\\{w_{i,\bot_{j}}\\}$, we have

$\displaystyle\begin{split}\mathcal{R}_{\mathcal{D}}[\mathbf{\delta}\left(\hat{y}|{\mathbf{x}}\right)]\geq{{\sum}}_{i=1}^{\text{K}}\bigintsss_{{\mathbf{x}}}\left({{\sum}}_{j=1}^{\text{K}}\delta(\hat{y}_{j}|{\mathbf{x}})w_{i,j}+{{\sum}}_{j=\text{K}+1}^{\text{K}+\text{J}}\delta(\hat{y}_{j}|{\mathbf{x}})w^{*}_{i,\bot}\right)\mathbb{P}\left({\textnormal{y}}=i\right)\mathbb{P}\left({\mathbf{x}}|{\textnormal{y}}=i\right)d{\mathbf{x}}.\end{split}$

We also denote
$\sum_{j=\text{K}+1}^{\text{K}+\text{J}}\delta(\hat{y}_{j}|{\mathbf{x}})$ as
$\delta(\hat{y}_{\bot}|{\mathbf{x}})$. Then the lower bound of
$\mathcal{R}_{\mathcal{D}}[\mathbf{\delta}\left(\hat{y}|{\mathbf{x}}\right)]$,
denoted as
$\bar{\mathcal{R}}_{\mathcal{D}}[\mathbf{\delta}\left(\hat{y}|{\mathbf{x}}\right)]$,
is

$\displaystyle\bar{\mathcal{R}}_{\mathcal{D}}[\mathbf{\delta}\left(\hat{y}|{\mathbf{x}}\right)]={{\sum}}_{i=1}^{\text{K}}\bigintsss_{{\mathbf{x}}}\left({{\sum}}_{j=1}^{\text{K}}\delta(\hat{y}_{j}|{\mathbf{x}})w_{i,j}+\delta(\hat{y}_{\bot}|{\mathbf{x}})w_{i,\bot}^{*}\right)\mathbb{P}\left({\textnormal{y}}=i\right)\mathbb{P}\left({\mathbf{x}}|{\textnormal{y}}=i\right)d{\mathbf{x}}.$

Since
$\sum_{j=1}^{\text{K}}\delta(\hat{y}_{j}|{\mathbf{x}})+\delta(\hat{y}_{\bot}|{\mathbf{x}})=1.0$
and
$\int_{{\mathbf{x}}}\mathbb{P}({\mathbf{x}}|{\textnormal{y}}=i)d{\mathbf{x}}=1.0$,
we can now decompose
$\bar{\mathcal{R}}_{\mathcal{D}}[\mathbf{\delta}\left(\hat{y}|{\mathbf{x}}\right)]$
in two terms:

$\bar{\mathcal{R}}_{\mathcal{D}}=\bar{\mathcal{R}}_{\mathcal{D}}^{\bot}+\bar{\mathcal{R}}_{\mathcal{D}}^{\mathbf{\delta}},$

where

$\displaystyle\begin{split}\bar{\mathcal{R}}_{\mathcal{D}}^{\bot}&=\sum_{i=1}^{\text{K}}\mathbb{P}({\textnormal{y}}=i)\cdot
w^{*}_{i,\bot},\\\
\bar{\mathcal{R}}_{\mathcal{D}}^{\mathbf{\delta}}&=\bigintss_{{\mathbf{x}}}\sum_{j=1}^{[\text{K}]\cup\\{\bot\\}}\delta(\hat{y}_{j}|{\mathbf{x}})Z_{j}({\mathbf{x}})d{\mathbf{x}},\text{and}\\\
Z_{j}({\mathbf{x}})&=\sum_{i=1}^{\text{K}}\left(w_{i,j}-w^{*}_{i,\bot}\right)\mathbb{P}({\mathbf{x}})\mathbb{P}({\textnormal{y}}=i|{\mathbf{x}}),\;j\in\\{1,2,\ldots,\text{K},\bot\\}.\end{split}$

To elaborate, we simplify the problem from deferring to multiple experts to
deferring to just the one expert with the minimum $w_{i,\bot_{j}}$ in
obtaining the lower bound $\bar{\mathcal{R}}_{\mathcal{D}}^{\bot}$. We also
observe that we have no control over $\bar{\mathcal{R}}_{\mathcal{D}}^{\bot}$.
However, we can control $\bar{\mathcal{R}}_{\mathcal{D}}^{\mathbf{\delta}}$ by
controlling the decision rule $\mathbf{\delta}$. We have
$Z_{\bot}({\mathbf{x}})=0$, and also it holds that

$\displaystyle\bar{\mathcal{R}}_{\mathcal{D}}^{\mathbf{\delta}}\geq\int_{{\mathbf{x}}}\min_{j}\left[Z_{j}\left({\mathbf{x}}\right)\right]d{\mathbf{x}},$

where the equality holds iff $\delta\left(\hat{y}_{k}|{\mathbf{x}}\right)=1.0$
for $k=\operatorname*{arg\,min}_{j}Z_{j}\left({\mathbf{x}}\right)$. Thus, the
optimal rule is to deterministically (i.e. with confidence 1.0) choose
$k\in\\{1,2,\ldots,\text{K},\bot\\}$ with the minimum
$Z_{j}\left({\mathbf{x}}\right)$. This means that choosing $j$ for which the
$Z_{j}({\mathbf{x}})$ is the smallest minimizes the risk. Given that
$Z_{\bot}=0$, this means that the classifier predicts when the minimum
$Z_{j}\left({\mathbf{x}}\right)$ is negative. Thus, deferral happens when
$Z_{j}({\mathbf{x}})$ is positive for all $j$, i.e., the optimal rejection
rule $r^{*}({\mathbf{x}})$ is:

$\displaystyle r^{*}({\mathbf{x}})$
$\displaystyle=\mathbb{I}\left[Z_{j}({\mathbf{x}})\geq 0;\forall
j\in\\{1,\ldots,\text{K}\\}\right].$ (17)

This rejection rule is similar to the learning to defer to one expert. Given
the definition of $Z_{j}({\mathbf{x}})$, the optimal behavior to choose which
expert to defer to is the one with minimum $w_{i,\bot{j}}$. We further
simplify the optimal rejection rule in the following proposition.

######  Proposition A.1.

The Bayes optimal rejection rule for L2D with multiple experts is given as:

$\displaystyle
r^{*}({\mathbf{x}})=\begin{cases}1\;\;\text{if}\;\;\mathbb{E}_{{\textnormal{y}}|{\mathbf{x}}}\left[\ell_{\text{clf}}(\hat{{\textnormal{y}}},{\textnormal{y}})\right]\geq\min_{j\in
J}\mathbb{E}_{{\textnormal{y}}|{\mathbf{x}}}\mathbb{E}_{{\textnormal{m}}|{\mathbf{x}},{\textnormal{y}}}\left[\ell_{\text{exp}}\left({\textnormal{m}}_{j},{\textnormal{y}}\right)\right]\;\;\forall\hat{{\textnormal{y}}}\in\mathcal{Y}\\\
0\;\;\text{otherwise}.\end{cases}$ (18)

Proof: The proof follows immediately from the definition of
$r^{*}({\mathbf{x}})$ in Equation 17.

In our work, we use the canonical $0$-$1$ loss for both $\ell_{\text{clf}}$
and $\ell_{\text{exp}}$. In this case, the rejection rule can trivially be
written as in the following corollary.

######  Corollary A.2.

For a misclassification $0$-$1$ loss, the optimal rejection rule is:

$\displaystyle r^{*}({\mathbf{x}})=\mathbb{I}\left[\max_{j\in
J}\mathbb{P}({\textnormal{m}}_{j}={\textnormal{y}}|{\mathbf{x}}={\bm{x}})\geq\max_{y\in\mathcal{Y}}\mathbb{P}({\textnormal{y}}=y|{\mathbf{x}}={\bm{x}})\right],$
(19)

where
$\mathbb{P}({\textnormal{m}}_{j}={\textnormal{y}}|{\mathbf{x}}={\bm{x}})$ is
the expert’s correctness probability for the $j^{\text{th}}$ expert, and
$\mathbb{P}({\textnormal{y}}=y|{\mathbf{x}}={\bm{x}})$ is the regular class
probability.

To sum it up, the Bayes optimal rule is to compare the confidences of the
experts and the classifier, and follow whosoever has the highest confidence.
The rule is analogous to the single expert setting proved in Mozannar and
Sontag (2020).

###  A.1 Proof of Theorem 3.1

Convexity of ${\uppsi}_{\text{SM}}^{J}$ is immediately clear. We provide the
proof for consistency below:

For simplicity, we denote

$\begin{split}&-\log\left(\frac{\exp\\{g_{y}({\bm{x}})\\}}{\sum_{y^{\prime}\in\mathcal{Y}^{\bot}}\exp\\{g_{y^{\prime}}({\bm{x}})\\}}\right)=\zeta_{y}({\bm{x}}),\\\
&-\log\left(\frac{\exp\\{g_{\bot,j}({\bm{x}})\\}}{\sum_{y^{\prime}\in\mathcal{Y}^{\bot}}\exp\\{g_{y^{\prime}}({\bm{x}})\\}}\right)=\zeta_{\bot,j}({\bm{x}}).\end{split}$
(20)

Then, ${\uppsi}_{\text{SM}}^{J}$ can be written as:

${\upphi}_{\text{SM}}^{J}\left(g_{1},\ldots,g_{K},g_{\bot,1},\ldots,g_{\bot,J};{\bm{x}},y,m_{1},\ldots,m_{J}\right)=\zeta_{y}({\bm{x}})+\sum_{j=1}^{J}\mathbb{I}[m_{j}=y]\cdot\zeta_{\bot,j}({\bm{x}}).$
(21)

We consider the pointwise risk $\mathcal{C}[{\uppsi}_{\text{SM}}^{J}]$ defined
as:

$\mathcal{C}[{\uppsi}_{\text{SM}}^{J}]=\mathbb{E}_{{\textnormal{y}}|{\mathbf{x}}={\bm{x}}}\mathbb{E}_{{\textnormal{m}}|{\mathbf{x}}={\bm{x}},{\textnormal{y}}=y}\left[{\uppsi}_{\text{SM}}^{J}(g_{1},\ldots,g_{K},g_{\bot,1},\ldots,g_{\bot,J};{\bm{x}},y,m_{1},\ldots,m_{J})\right],$
(22)

where ${\textnormal{m}}|{\mathbf{x}}={\bm{x}},{\textnormal{y}}=y$ is a compact
representation for each
${\textnormal{m}}_{j}|{\mathbf{x}}={\bm{x}},{\textnormal{y}}=y$. Our setup
assumes that each ${\textnormal{m}}_{j}$ is independent. Denote
$\eta_{y}({\bm{x}})=\mathbb{P}({\textnormal{y}}=y|{\mathbf{x}}={\bm{x}})$, we
expand the expectations:

$\displaystyle\mathcal{C}[{\uppsi}_{\text{SM}}^{J}]$
$\displaystyle=\sum_{y\in\mathcal{Y}}\eta_{y}({\bm{x}})\cdot\zeta_{y}({\bm{x}})+\sum_{j=1}^{J}\left(\sum_{y\in\mathcal{Y}}\eta_{y}({\bm{x}})\sum_{m_{j}\in\mathcal{M}}\mathbb{P}({\textnormal{m}}_{j}=m_{j}|{\mathbf{x}}={\bm{x}},{\textnormal{y}}=y)\mathbb{I}[m_{j}=y]\cdot\zeta_{\bot,j}({\bm{x}})\right).$
$\displaystyle\mathcal{C}[{\uppsi}_{\text{SM}}^{J}]$
$\displaystyle=\sum_{y\in\mathcal{Y}}\eta_{y}({\bm{x}})\cdot\zeta_{y}({\bm{x}})+\sum_{j=1}^{J}\left(\sum_{y\in\mathcal{Y}}\eta_{y}({\bm{x}})\sum_{m_{j}\in\mathcal{M}}\mathbb{P}({\textnormal{m}}_{j}=y|{\mathbf{x}}={\bm{x}},{\textnormal{y}}=y)\cdot\zeta_{\bot,j}({\bm{x}})\right).$
$\displaystyle\mathcal{C}[{\uppsi}_{\text{SM}}^{J}]$
$\displaystyle=\sum_{y\in\mathcal{Y}}\eta_{y}({\bm{x}})\cdot\zeta_{y}({\bm{x}})+\sum_{j=1}^{J}\mathbb{P}({\textnormal{m}}_{j}={\textnormal{y}}|{\mathbf{x}}={\bm{x}})\cdot\zeta_{\bot,j}({\bm{x}}).$

Next, we consider the minimizer of $\mathcal{C}[{\uppsi}_{\text{SM}}^{J}]$.
Since we have established convexity, we can analyze the minimizers of
$\mathcal{C}[{\uppsi}_{\text{SM}}^{J}]$ by taking the partial derivatives
w.r.t. $g_{y}\\{{\bm{x}}\\}$ and $g_{\bot,j}\\{{\bm{x}}\\}$ respectively and
set them to $0$.

Thus, w.r.t. $g_{y}\\{{\bm{x}}\\}$, we have

$\begin{split}\frac{\partial\mathcal{C}[{\uppsi}_{\text{SM}}^{J}]}{\partial
g_{y}\\{{\bm{x}}\\}}=0\implies\frac{\exp\\{g_{y}({\bm{x}})\\}}{\sum_{y^{\prime}\in\mathcal{Y}^{\bot}}\exp\\{g_{y^{\prime}}({\bm{x}})\\}}=\frac{\mathbb{P}({\textnormal{y}}=y|{\mathbf{x}}=x)}{1+\sum_{j=1}^{J}\mathbb{P}({\textnormal{m}}_{j}={\textnormal{y}}|{\mathbf{x}}={\bm{x}})}.\end{split}$
(23)

Similarly, w.r.t. $g_{\bot,j}\\{{\bm{x}}\\}$ we have

$\begin{split}\frac{\partial\mathcal{C}[{\uppsi}_{\text{SM}}^{J}]}{\partial
g_{\bot,j}\\{{\bm{x}}\\}}=0\implies\frac{\exp\\{g_{\bot,j}({\bm{x}})\\}}{\sum_{y^{\prime}\in\mathcal{Y}^{\bot}}\exp\\{g_{y^{\prime}}({\bm{x}})\\}}=\frac{\mathbb{P}({\textnormal{m}}_{j}={\textnormal{y}}|{\mathbf{x}}={\bm{x}})}{1+\sum_{j=1}^{J}\mathbb{P}({\textnormal{m}}_{j}={\textnormal{y}}|{\mathbf{x}}={\bm{x}})}.\end{split}$
(24)

The above equations hold true for optimal classifier and the rejector. Thus,
if we take the decision as in the main text, we are agreeing with the Bayes
solution (considering that denominators are same in both the above
conditions).

###  A.2 Proof of Theorem 3.2

The proof follows directly from Verma and Nalisnick (2022). However, for
completion we provide the full proof below.

For the surrogate prediction functions
$g_{1},\ldots,g_{K},g_{\bot,1},\ldots,g_{\bot,J}$, and the binary
classification surrogate loss
${\upphi}:\\{-1,1\\}\times\mathbb{R}\rightarrow\mathbb{R}_{+}$,
${\uppsi}_{\text{OvA}}^{J}$ takes the following pointwise-form:

$\begin{split}\mathcal{C}[{\uppsi}_{\text{OvA}}^{J}]={\uppsi}_{\text{OvA}}^{J}(g_{1},\ldots,g_{K},g_{\bot,1},\ldots,g_{\bot,J};{\bm{x}},y,m_{1},\ldots,m_{J})&={\upphi}[g_{y}({\bm{x}})]+\sum_{y^{\prime}\in\mathcal{Y},y^{\prime}\neq
y}{\upphi}[-g_{y^{\prime}}({\bm{x}})]\ \
+\sum_{j=1}^{J}{\upphi}[-g_{\bot,j}({\bm{x}})]\\\
&+\sum_{j=1}^{J}\mathbb{I}[m_{j}=y]\left({\upphi}[g_{\bot,j}({\bm{x}})]-{\upphi}[-g_{\bot,j}({\bm{x}})]\right).\end{split}$

We consider the pointwise inner ${\uppsi}_{\text{OvA}}$risk for some
${\mathbf{x}}={\bm{x}}$ written as follows:

$\mathcal{C}[{\uppsi}_{\text{OvA}}^{J}]=\mathbb{E}_{{\textnormal{y}}|{\mathbf{x}}={\bm{x}}}\mathbb{E}_{{\textnormal{m}}|{\mathbf{x}}={\bm{x}},{\textnormal{y}}=y}\left[{\uppsi}_{\text{OvA}}^{J}(g_{1},\ldots,g_{K},g_{\bot,1},\ldots,g_{\bot,J};{\bm{x}},y,m_{1},\ldots,m_{J})\right],$
(25)

We expand both the expectations one-by-one below:

$\begin{split}\mathcal{C}[{\uppsi}_{\text{OvA}}^{J}]&=\mathbb{E}_{{\textnormal{y}}|{\mathbf{x}}={\bm{x}}}\Biggr{[}{\upphi}[g_{y}({\bm{x}})]+\sum_{y^{\prime}\in\mathcal{Y},y^{\prime}\neq
y}{\upphi}[-g_{y^{\prime}}({\bm{x}})]\ \
+\sum_{j=1}^{J}{\upphi}[-g_{\bot,j}({\bm{x}})]\\\
&+\sum_{j=1}^{J}\left(\sum_{m_{j}\in\mathcal{M}}\mathbb{P}({\textnormal{m}}_{j}=m_{j}|{\mathbf{x}}={\bm{x}},{\textnormal{y}}=y)\mathbb{I}\left[m_{j}=y\right]\left({\upphi}[g_{\bot,j}({\bm{x}})]-{\upphi}[-g_{\bot,j}({\bm{x}})]\right)\right)\Biggr{]}.\end{split}$

Denote $\mathbb{P}({\textnormal{y}}=y|{\mathbf{x}}={\bm{x}})$ as
$\eta_{y}({\bm{x}})$, then

$\displaystyle\begin{split}\mathcal{C}[{\uppsi}_{\text{OvA}}^{J}]&=\mathbb{E}_{{\textnormal{y}}|{\mathbf{x}}={\bm{x}}}\sum_{y\in\mathcal{Y}}\left[{\upphi}[g_{y}({\bm{x}})]+\sum_{y^{\prime}\in\mathcal{Y},y^{\prime}\neq
y}{\upphi}[-g_{y^{\prime}}({\bm{x}})]\right]\\\
&+\sum_{j=1}^{J}\left({\upphi}[-g_{\bot,j}({\bm{x}})]+\sum_{y\in\mathcal{Y}}\eta_{y}({\bm{x}})\left(\sum_{m_{j}\in\mathcal{M}}\mathbb{P}({\textnormal{m}}_{j}=m_{j}|{\mathbf{x}}={\bm{x}},{\textnormal{y}}=y)\mathbb{I}\left[m_{j}=y\right]\left({\upphi}[g_{\bot,j}({\bm{x}})]-{\upphi}[-g_{\bot,j}({\bm{x}})]\right)\right)\right)\\\
\mathcal{C}[{\uppsi}_{\text{OvA}}^{J}]&=\sum_{y\in\mathcal{Y}}\eta_{y}({\bm{x}})\left[{\upphi}[g_{y}({\bm{x}})]+\sum_{y^{\prime}\in\mathcal{Y},y^{\prime}\neq
y}{\upphi}[-g_{y^{\prime}}({\bm{x}})]\right]\\\
&+\sum_{j=1}^{J}\left({\upphi}[-g_{\bot,j}({\bm{x}})]+\underbrace{\sum_{y\in\mathcal{Y}}\eta_{y}({\bm{x}})\sum_{m_{j}\in\mathcal{M}}\mathbb{P}({\textnormal{m}}_{j}=y|{\mathbf{x}}={\bm{x}},{\textnormal{y}}=y)}_{\text{$\mathbb{P}\left({\textnormal{m}}_{j}={\textnormal{y}}|{\mathbf{x}}={\bm{x}}\right)$}}\left({\upphi}[g_{\bot,j}({\bm{x}})]-{\upphi}[-g_{\bot,j}({\bm{x}})]\right)\right)\\\
\mathcal{C}[{\uppsi}_{\text{OvA}}^{J}]&=\sum_{y\in\mathcal{Y}}\eta_{y}({\bm{x}})\left[{\upphi}[g_{y}({\bm{x}})]+\sum_{y^{\prime}\in\mathcal{Y},y^{\prime}\neq
y}{\upphi}[-g_{y^{\prime}}({\bm{x}})]\right]\\\
&+\sum_{j=1}^{J}\bigg{[}{\upphi}[-g_{\bot,j}({\bm{x}})]+\mathbb{P}\left({\textnormal{m}}_{j}={\textnormal{y}}|{\mathbf{x}}={\bm{x}}\right)\left({\upphi}[g_{\bot,j}({\bm{x}})]-{\upphi}[-g_{\bot,j}({\bm{x}})]\right)\bigg{]}.\end{split}$

Denote
$\mathbb{P}({\textnormal{m}}_{j}={\textnormal{y}}|{\mathbf{x}}={\bm{x}})$ as
$p_{{\textnormal{m}}_{j}}$, then we have

$\begin{split}\mathcal{C}[{\uppsi}_{\text{OvA}}^{J}]&=\sum_{y\in\mathcal{Y}}\eta_{y}({\bm{x}})\left[{\upphi}[g_{y}({\bm{x}})]+\sum_{y^{\prime}\in\mathcal{Y},y^{\prime}\neq
y}{\upphi}[-g_{y^{\prime}}({\bm{x}})]\right]\\\
&+\sum_{j=1}^{J}\left[\left(1-p_{{\textnormal{m}}_{j}}\right){\upphi}[-g_{\bot,j}({\bm{x}})]+p_{{\textnormal{m}}_{j}}{\upphi}[g_{\bot,j}({\bm{x}})]\right]\end{split}$

We further simplify the above equation as follows:

$\begin{split}\mathcal{C}[{\uppsi}_{\text{OvA}}^{J}]&=\sum_{y\in\mathcal{Y}}\left[\eta_{y}({\bm{x}})\cdot{\upphi}[g_{y}({\bm{x}})]+\left(1-\eta_{y}({\bm{x}})\right)\cdot{\upphi}[-g_{y}({\bm{x}})]\right]+\sum_{j=1}^{J}\left[\left(1-p_{{\textnormal{m}}_{j}}\right){\upphi}[-g_{\bot,j}({\bm{x}})]+p_{{\textnormal{m}}_{j}}{\upphi}[g_{\bot,j}({\bm{x}})]\right]\end{split}$
(26)

Thus, we conclude from the above expression that we $\text{K}+\text{J}$ binary
classification problems where the pointwise risk (or inner risk) for the
$i^{th}$ binary classification problem is given as
$\eta_{y}\left({\bm{x}}\right){\upphi}\left(g_{y}\left({\bm{x}}\right)\right)+\left(1-\eta_{y}\left({\bm{x}}\right)\right){\upphi}\left(-g_{y}\left({\bm{x}}\right)\right)$
for $i\in[K]$ and
$p_{{\textnormal{m}}_{j}}\left({\bm{x}}\right){\upphi}\left(g_{\bot,j}\left({\bm{x}}\right)\right)+\left(1-p_{{\textnormal{m}}_{j}}\left({\bm{x}}\right)\right){\upphi}\left(-g_{\bot,j}\left(x\right)\right)$
when $i\in[J]\\}$. Thus, minimizer of the inner ${\uppsi}_{\text{OvA}}$ -risk
can be analyzed in terms of the pointwise minimizer of the inner
${\upphi}$-risk for each of the $\text{K}+\text{J}$ sub binary classification
problems. Denote the minimizer of pointwise inner ${\uppsi}_{\text{OvA}}$-risk
as $\bm{g}^{*}$, then the above decomposition means $g^{*}_{i}$ corresponds to
the minimizer of the inner ${\upphi}$-risk for the $i$th binary classification
problem.

We know that the Bayes solution for the binary classification problem is
$\operatorname{sign}\left(\eta({\bm{x}})-\frac{1}{2}\right)$ where
$\eta({\bm{x}})$ denotes $p({\textnormal{y}}=1|{\mathbf{x}}={\bm{x}})$. Now
when the binary surrogate loss ${\upphi}$ is a strictly proper composite loss
for binary classification, by the property of strictly proper composite
losses, we have $\operatorname{sign}(g^{*}_{y}({\bm{x}}))$ would agree with
the Bayes solution of the Binary classification, i.e. $g^{*}_{y}({\bm{x}})>0$
if $\eta_{y}\left({\bm{x}}\right)>\frac{1}{2}$. And similarly
$g^{*}_{\bot}\left({\bm{x}}\right)>0$ if
$p_{{\textnormal{m}}_{j}}\left({\bm{x}}\right)>\frac{1}{2}$. Furthermore, we
have the existence of a continuous and increasing inverse link function
${\upgamma}^{-1}$ for the binary surrogate ${\upphi}$ with the property that
${\upgamma}^{-1}\left(g^{*}_{y}\left({\bm{x}}\right)\right)$ would converge to
$\eta_{y}\left({\bm{x}}\right)$. Similarly,
${\upgamma}^{-1}\left(g^{*}_{\bot,j}\left({\bm{x}}\right)\right)$ would
converge to $p_{{\textnormal{m}}_{j}}\left({\bm{x}}\right)$.

Thus, when the binary surrogate loss ${\upphi}$ is a strictly proper composite
loss, and the classifier and the rejector are defined as in the main text, the
minimizer of the pointwise risk
$\mathcal{C}\left[{\uppsi}_{\text{OvA}}^{J}\right]$ agree with the Bayes
optimal solution. Thus, ${\uppsi}_{\text{OvA}}^{J}$ is a calibrated loss
function for L2D w.r.t. $0$-$1$ misclassification loss.

##  Appendix B Derivations

###  B.1 Derivation of Equation 11

Equation 11 is the same as Equation 24.

##  Appendix C Choice of Hyperparameters for Regularized Conformal Ensembles

We have two hyperparameters $k_{\text{reg}}$ and $\beta$ in the regularized
conformal ensemble procedure discussed in Section 5. Algorithm 1 below details
the procedure to choose $k_{\text{reg}}$.

Algorithm 1 Choice of $k_{\text{reg}}$

Data: Error rate: $\alpha$, $\\#$ Experts: J, Data:
$\\{\mathbf{s}^{i}({\bm{x}}),\mathbf{e}^{i}({\bm{x}})\\}_{i=1}^{\text{N}}$
where $\mathbf{s}^{i}({\bm{x}})\in[0,1]^{\text{J}}$ and
$\mathbf{e}^{i}({\bm{x}})\in\\{0,1\\}^{\text{J}}$

$B\leftarrow\Lbag\cdot\Rbag$ for  _i $\in$ $[1,N]$_ do

for  _j $\in$ $[1,J]$_ do

if  _$\mathbb{I}\left\\{e_{j}^{i}({\bm{x}})==1\right\\}$_ then

$B\leftarrow B\cup\\{s_{j}^{i}({\bm{x}})\\}$

end if

end for

$S\leftarrow\text{sort}\left(B\right)$ s.t. $u_{i},u_{j}\in S,i\leq j$, then
$u_{i}\geq u_{j}$ $k_{\text{reg}}^{*}\leftarrow 1-\alpha$ quantile of $S$

end for

We can employ corrections to account for finite sample size N on line $9$.
Given this procedure to choose $k^{*}_{\text{reg}}$, one may argue that our
choice of $k^{*}_{\text{reg}}$ can give us meaningful prediction sets by
designing a prediction set as:

$C_{2}({\bm{x}})=\\{j:s_{j}({\bm{x}})\geq k^{*}_{\text{reg}}\\}.$

However, our next proposition establishes that $C_{\lambda}({\bm{x}})$ results
in prediction sets at most as large as $C_{2}({\bm{x}})$.

(a) Set size.

(b) Sys. Acc, naive conformal.

(c) Sys. Acc, regularized conformal.

Figure 4: Gradual overlapping expertise results. The figures above report the
average set size (a) and system accuracies under the naive conformal method
(b) and regularized conformal method (c) for increasing expertise overlap for
CIFAR-10. From Figure (a) we see that the regularized conformal method is more
dynamic than the naive method for both OvA and softmax. System accuracies for
the naive conformal method (b) are slightly better than fixed-size ensemble
because we ensemble all experts, whereas for the regularized conformal method
(c) we see a slightly drop due to the adaptivity of the conformal sets

######  Proposition C.1.

Define the prediction sets
$C_{\lambda}({\bm{x}})=\\{j:s_{j}({\bm{x}})+\beta\left(s_{j}({\bm{x}})-k^{*}_{\text{reg}}\right)>1-\lambda\\}$
and $C_{2}({\bm{x}})=\\{j:s_{j}({\bm{x}})\geq k^{*}_{\text{reg}}\\}$, where
$k^{*}_{\text{reg}}$ is defined as above, $\beta\geq 0$, $0\leq\lambda\leq 1$,
then it trivially holds that

$C_{\lambda}({\bm{x}})\subseteq C_{2}({\bm{x}}).$

We tune $\beta$ in a grid-search manner to optimize for prediction set sizes.
We split the total number of deferred samples into two portions: one for
tuning hyperparameters $\beta$ and $k_{\text{reg}}$, another for the regular
conformal procedure.

##  Appendix D Additional Experiments and Results

###  D.1 Overlapping expertise among experts

Similar to the conducted experiments in Section 7.3 in the main manuscript, we
want to study ensembling multiple experts, this time under a different
experiment setup. We will have 10 experts for the CIFAR-10 dataset, each of
them being an oracle on a specific class out of the $10$ classes from the
dataset, and we will increase the overlapping probability of these experts
being correct on the other classes where the experts are not oracle from
$10\%$ to $95\%$ overlap. That is, we will vary from specialized experts to
fully overlapped experts. We hope to see that the average set size for
specialized experts is close to 1 and for fully overlapped experts close to
the total number of experts. We report the results in Figure 4.

##### Average Set Size

First of all, it is worth noticing from Figure 4(a) that the average set size
for the naive conformal method, both for softmax and OvA, is always close to
the total number of experts, for specialized and overlapped experts. If we
remember, the naive test statistic is calculated among all correct experts.
This is a very important point, because if an expert happens to be correct
outside of their expertise domain, this results in a very big non-conformity
score because of the low confidence of such expert. That is, imagine for
certain sample ${\bm{x}}$ and class $y=3$, for low overlapping probabilities,
we might have $E=3$, where $e=1$ could be the oracle for class $y=3$ and
$e=2,e=3$ two experts that were correct by chance. From Equation 14 in the
manuscript, we can expect that, best-case scenario
$s_{\pi_{1}}>s_{\pi_{2}}>s_{\pi_{3}}$, or even worse
$s_{\pi_{1}}>\dots>s_{\pi_{8}}>s_{\pi_{2}}>s_{\pi_{3}}$, because other
experts’ confidences could be also greater than confidences from correct
experts by chance. Therefore, we will obtain bigger test-statistics that
result in very larger conformal sets. This problem has already been addressed
in Angelopoulos et al. (2022). However, notice how the regularized test
statistic is capable of producing smaller sets. The idea is that now we
optimize additional parameters (described in Section 5 in the manuscript) to
ensure that confidences lower than a certain threshold are filtered out for
the calculation of the test statistic.

##### System Accuracies

In Figure 4(b) and 4(c) we report the system accuracies for the naive
conformal method and the regularized conformal method respectively. For the
naive conformal method we obtain better results than using a fixed-size
ensemble of experts of size 5. Because set sizes are almost always close to
the total number of experts, and we do majority voting, then as long as the is
a correct expert plus an expert correct by chance, we will predict correctly.
However, for the regularized conformal method we notice a drop in the system
accuracy for lower overlapping probabilities. Since for such cases now the set
sizes are smaller, we have smaller number of experts in the set and therefore
less chance of having correct experts by chance in the ensemble. Despite this
drop in the accuracy, we clearly have a more dynamic and less conservative
ensembling method.


* Roming et al. (2005) Roming, P. W. A., Kennedy, T. E., Mason, K. O., et al. 2005, Space Sci. Rev., 120, 95, doi: 10.1007/s11214-005-5095-4
  * Rosolowsky et al. (2021) Rosolowsky, E., Hughes, A., Leroy, A. K., et al. 2021, MNRAS, 502, 1218, doi: 10.1093/mnras/stab085
  * Ryder et al. (1999) Ryder, S. D., Walsh, W., & Malin, D. 1999, Publications of the Astronomical Society of Australia, 16, 84–88, doi: 10.1071/AS99084
  * Sahu et al. (2019) Sahu, S., Subramaniam, A., Simunovic, M., et al. 2019, ApJ, 876, 34, doi: 10.3847/1538-4357/ab11d0
  * Salim et al. (2007) Salim, S., Rich, R. M., Charlot, S., et al. 2007, ApJS, 173, 267, doi: 10.1086/519218
  * Salo et al. (2015) Salo, H., Laurikainen, E., Laine, J., et al. 2015, ApJS, 219, 4, doi: 10.1088/0067-0049/219/1/4
  * Scheuermann et al. (2023) Scheuermann, F., Kreckel, K., Barnes, A. T., et al. 2023, MNRAS, 522, 2369, doi: 10.1093/mnras/stad878
  * Schinnerer et al. (2023) Schinnerer, E., Emsellem, E., Henshaw, J. D., et al. 2023, ApJ, 944, L15, doi: 10.3847/2041-8213/acac9e
  * Schlegel et al. (1998) Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525, doi: 10.1086/305772
  * Singh et al. (2014) Singh, K. P., Tandon, S. N., Agrawal, P. C., et al. 2014, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9144, Space Telescopes and Instrumentation 2014: Ultraviolet to Gamma Ray, ed. T. Takahashi, J.-W. A. den Herder, & M. Bautz, 91441S, doi: 10.1117/12.2062667
  * Singh et al. (2021) Singh, S., Ashby, M. L. N., Vig, S., et al. 2021, MNRAS, 504, 4143, doi: 10.1093/mnras/stab1048
  * Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163, doi: 10.1086/498708
  * Stetson (1987) Stetson, P. B. 1987, PASP, 99, 191, doi: 10.1086/131977
  * Stuber et al. (2023) Stuber, S. K., Schinnerer, E., Williams, T. G., et al. 2023, A&A, 676, A113, doi: 10.1051/0004-6361/202346318
  * Subramaniam et al. (2016) Subramaniam, A., Sindhu, N., Tandon, S. N., et al. 2016, The Astrophysical Journal, 833, L27, doi: 10.3847/2041-8213/833/2/l27
  * Sun et al. (2020) Sun, J., Leroy, A. K., Ostriker, E. C., et al. 2020, ApJ, 892, 148, doi: 10.3847/1538-4357/ab781c
  * Sun et al. (2022a) Sun, J., Leroy, A. K., Rosolowsky, E., et al. 2022a, AJ, 164, 43, doi: 10.3847/1538-3881/ac74bd
  * Sun et al. (2022b) —. 2022b, arXiv e-prints, arXiv:2206.07055.  https://arxiv.org/abs/2206.07055
  * Tandon et al. (2017) Tandon, S. N., Subramaniam, A., Girish, V., et al. 2017, AJ, 154, 128, doi: 10.3847/1538-3881/aa8451
  * Tandon et al. (2020) Tandon, S. N., Postma, J., Joseph, P., et al. 2020, AJ, 159, 158, doi: 10.3847/1538-3881/ab72a3
  * Thilker et al. (2007) Thilker, D. A., Bianchi, L., Meurer, G., et al. 2007, ApJS, 173, 538, doi: 10.1086/523853
  * Tibshirani (1996) Tibshirani, R. 1996, Journal of the Royal Statistical Society. Series B (Methodological), 58, 267.  http://www.jstor.org/stable/2346178
  * Tilanus & Allen (1993) Tilanus, R. P. J., & Allen, R. J. 1993, A&A, 274, 707
  * Ujjwal et al. (2022) Ujjwal, K., Kartha, S. S., Subramanian, S., et al. 2022, MNRAS, 516, 2171, doi: 10.1093/mnras/stac2285
  * Vazdekis et al. (2016) Vazdekis, A., Koleva, M., Ricciardelli, E., Röck, B., & Falcón-Barroso, J. 2016, MNRAS, 463, 3409, doi: 10.1093/mnras/stw2231
  * Véron-Cetty & Véron (2010) Véron-Cetty, M. P., & Véron, P. 2010, A&A, 518, A10, doi: 10.1051/0004-6361/201014188
  * Villar et al. (2011) Villar, V., Gallego, J., Pérez-González, P. G., et al. 2011, ApJ, 740, 47, doi: 10.1088/0004-637X/740/1/47
  * Weisz et al. (2012) Weisz, D. R., Johnson, B. D., Johnson, L. C., et al. 2012, ApJ, 744, 44, doi: 10.1088/0004-637X/744/1/44
  * Whitmore et al. (2023) Whitmore, B. C., Chandar, R., Rodríguez, M. J., et al. 2023, ApJ, 944, L14, doi: 10.3847/2041-8213/acae94
  * Zurita & Pérez (2008) Zurita, A., & Pérez, I. 2008, A&A, 485, 5, doi: 10.1051/0004-6361:20079129