algorithmic [name=Theorem,numberwithin=section]theorem
[name=Lemma,numberwithin=section]lemma

# Conformal Prediction via
Regression-as-Classification

Etash Guha
RIKEN Center for AI Project, SambaNova Systems
<EMAIL_ADDRESS>&Shlok Natarajan
Salesforce
<EMAIL_ADDRESS>&Thomas Möllenhoff
RIKEN Center for AI Project
<EMAIL_ADDRESS>&Mohammad Emtiyaz Khan
RIKEN Center for AI Project
<EMAIL_ADDRESS>
Eugene Ndiaye
Apple
<EMAIL_ADDRESS>
###### Abstract

Conformal prediction (CP) for regression can be challenging, especially when
the output distribution is heteroscedastic, multimodal, or skewed. Some of the
issues can be addressed by estimating a distribution over the output, but in
reality, such approaches can be sensitive to estimation error and yield
unstable intervals. Here, we circumvent the challenges by converting
regression to a classification problem and then use CP for classification to
obtain CP sets for regression. To preserve the ordering of the continuous-
output space, we design a new loss function and make necessary modifications
to the CP classification techniques. Empirical results on many benchmarks
shows that this simple approach gives surprisingly good results on many
practical problems.

##  1 Introduction

Quantifying and estimating the uncertainty of machine-learning models is an
important task for many problems, especially mission-critical applications
where reliable predictions are required. Conformal Prediction (CP) (Vovk et
al., 2005) has recently gained popularity and has been used successfully in
applications such as breast cancer detection (Lambrou et al., 2009), stroke
risk prediction (Lambrou et al., 2010), and drug discovery (Cortés-Ciriano &
Bender, 2020). Under mild conditions, CP techniques aim to construct a
prediction set that, for given test inputs, is guaranteed to contain the true
(unknown) output with high probability. The set is built using a conformity
score, which, roughly speaking, indicates the similarity between a new test
example and the training examples. The conformal set merely gathers examples
that have large conformity scores. Despite its popularity, CP for regression
can be challenging, especially when the output distribution is
heteroscedastic, multimodal, or skewed (Lei & Wasserman, 2014). The main
challenge lies in the design of the conformity score. It is common to use a
simple choice for score functions such as distance to mean regressor, but such
choices may ignore the subtle features of the shape of the output
distribution. For instance, this could lead to symmetric intervals or ignoring
the heteroscedasticity. In theory, it is better to estimate the (conditional)
distribution over the output, for example, by using kernel density estimation
and directly using it to build a confidence interval. However, such estimation
approaches are also challenging, and estimates can be sensitive to the choice
of kernel and hyperparameters, which can yield unstable results.

We circumvent the challenges by exploiting the existing CP techniques for
classification. We proceed by first converting regression to a classification
problem and then using CP techniques for classification to obtain a conformal
set. Regression-as-classification approaches are popular for various
applications in computer vision and have led to more accurate training than
only-regression training (Stewart et al., 2023; Zhang et al., 2016; Fu et al.,
2018; Rothe et al., 2015; Van Den Oord et al., 2016; Diaz & Marathe, 2019). We
leverage them to construct a distribution-based conformal set that can
flexibly capture the shape of the output distribution while preserving the
simplicity and efficiency of CP for classification. First, we discretize the
output space into bins, treating each bin as a distinct class. Second, to
preserve the ordering of the continuous output space, we design an alternative
loss function that penalizes the density on bins far from the true output
value but also facilitates variability by using an entropy regularization. The
loss design is similar in spirit to Weigend & Srivastava (1995); Diaz &
Marathe (2019). The resulting method can adapt to heteroscedasticity,
bimodality, or both in the label distribution. We verify this on synthetic and
real datasets where we achieve the shortest intervals compared to other CP
baselines. See examples in Figure 1.

(a) Heteroskedastic

(b) Bimodal

Figure 1: We show two examples where the output distribution is
heteroskedastic (left) and bimodal (right). In both cases, our method is able
to change the interval (shaded gray region) adaptively as the input values $x$
are increased. Examples outside the gray regions (white dots) are deemed
different from those inside it (black dots).

##  2 Background on Conformal Prediction

Given a new input $x_{\rm new}$, CP techniques aim to construct a set that
contains the true but unknown output $y_{\rm new}$ with high probability.
Assuming that a pair of input-output variables $(x,y)$ has a joint density
$p(x,y)$ and a conditional density $p(y\mid x)$, oracle prediction sets (with
joint and conditional coverage) for the output $y$ can be constructed as

$\displaystyle\\{z\in\mathbb{R}:p(x,z)\geq\tau_{\alpha}\\}~{}\text{ or
}~{}\\{z\in\mathbb{R}:p(z\mid x)\geq\tau_{\alpha,x}\\},$ (1)

where the thresholds $\tau_{\alpha}$ and $\tau_{\alpha,x}$ are selected to
ensure that the corresponding sets have a probability mass that meets or
exceeds prescribed confidence level $1-\alpha\in(0,1)$. As the ground-truth
distribution is unknown, we rely solely on estimating these uncertainty sets
using the density estimators $\hat{p}(x,y)$ and $\hat{p}(y\mid x)$. The latter
can be inaccurate due to numerous sources of errors such as model
misspecification, small sample size, high optimization errors during training,
and overfitting. Without a stronger distribution assumption, the finite-sample
guarantee is typically not upheld.

Conformal Prediction has arisen as a method for yielding sets that do hold
finite-sample guarantees. Given a partially observed instance $(x_{\rm
new},y_{\rm new})$ where $y_{\rm new}$ is unknown, Conformal Prediction (CP)
(Vovk et al., 2005) constructs a set of values that contains $y_{\rm new}$
with high probability without knowing the underlying data distribution. Under
conformal prediction, this property is guaranteed under the mild assumption
that the data satisfies exchangeability. The set is called the conformal set
and is built using a conformity score, denoted by $\sigma(x,y)$, which
measures how appropriate an output value is for a given input example. There
are many ways to build the conformity score, but they all involve splitting
the data into a training set $\mathcal{D}_{\rm tr}$ and a calibration set
$\mathcal{D}_{\rm cal}$. Often, a prediction model $\mu_{\text{tr}}(x)$ is
built using the training set, and then a conformity score is obtained using
this model along with the calibration set (we will shortly give an example).
The conformal set merely gathers the points with larger conformity scores:

$\\{z\in\mathbb{R}:\,\sigma(x_{\rm new},z)\geq Q_{1-\alpha}(\mathcal{D}_{\rm
cal})\\},$

where $Q_{1-\alpha}(\mathcal{D}_{\rm cal})$ is the $(1-\alpha)$-quantile of
the conformity scores on the calibration data. This set provably contains
$y_{\rm new}$ with probability larger than $1-\alpha$ for any finite sample
size and without assumption on the ground-truth distribution.

There are many design choices for this conformity score. For example, one can
choose a prediction model $\mu_{\rm tr}(x)$ as an estimate of the conditional
expectation and measure conformity as the absolute value of the residual,
i.e., $\sigma(x,y)=-|y-\mu_{\rm tr}(x)|$. The corresponding conformal set is a
single interval centered around the prediction $\mu_{\rm tr}(x)$ and of
constant length $Q_{1-\alpha}(\mathcal{D}_{\rm cal})$ for any example $x_{\rm
new}$, without taking into account its variability. However, in situations
where the underlying data distribution demonstrates skewness or
heteroscedasticity, we may desire a more flexible conformity score. We wish to
generate an approach to approximate this density function that is versatile in
the distributions it can represent and avoids the difficulties seen in prior
methods. We explore established density estimation in Classification Conformal
Prediction that are already performing effectively.

#### Related Works

Since their introduction, a lot of work has been done to improve the set of
conformal predictions. As simple score function, distance to conditional mean
ie $\sigma(x,y)=|y-\mu(x)|$ where $\mu(x)$ is an estimate of $\mathbb{E}(y\mid
x)$ was prominently used (Papadopoulos et al., 2002; Lei et al., 2018; Guha et
al., 2023). Instead, (Romano et al., 2019) suggests estimating a conditional
quantile instead and a conformity score function based on the distance from a
trained quantile regressor, i.e.
$\sigma(x,y)=\max(\mu_{\alpha/2}(x)-y,y-\mu_{1-\alpha/2}(x))$ where
$\mu_{\alpha}$ are the $\alpha$-th quantile regressors. Within the literature,
our strategy is more closely related to the distribution-based methods (Lei et
al., 2013). Following a similar line of work (Chernozhukov et al., 2021)
argued that the conformal quantile regression score function might be less
adaptive since the distance of the quantile behaves similarly to the distance
of the mean estimate. Instead, they suggested estimating the cumulative
(conditional) distribution function and directly outputting
$\\{y:F(x,y)\in[\alpha/2,1-\alpha/2]\\}$. However, their method cannot account
for bimodality since it can only output a single interval by design. An
equally interesting approach for distribution-based prediction sets is based
on learning a Bayesian estimator, which, however, may require a well-specified
prior and can be computationally expensive (Fong & Holmes, 2021). Variants
based on estimating the conditional density of the response using histogram
regressors (Sesia & Romano, 2021) could detect a possible asymmetry of the
ground-truth distribution. However, our densities are estimated through
classification techniques, whereas their densities are learned through a
histogram of many regressors. Smoothness over the distribution is encoded in
our loss function, whereas they use a post hoc subinterval finding algorithm
through linear programming to prevent many disjoint intervals. Our techniques
have some overlap with ordinal regression and ordinal classification. Xu et
al. (2023) discusses various risk categories (similar to coverage) for ordinal
classification, while our work considers different score functions where
coverage serves as the loss function. Lu et al. (2022) discuss the adaptation
of APS style score functions to accommodate the ordinal structure of classes,
which was subsequently utilised in (Sesia & Romano, 2021).

##  3 CP via Regression-as-Classification

###  3.1 Classification Conformal Prediction

We aim to compute a conformity function that accurately predicts the
appropriateness of a label for a specific data point. Given that the
distribution function across labels can adopt diverse forms, such as being
bimodal, heavy-tailed, or heteroscedastic, our approach must effectively
factor in such shapes while upholding its coverage precision. A frequently
employed technique involves using the conditional label density as a
conformity function, leading to reliable results in the classification
context. Typically, practitioners perform conformal prediction for
classification with probability estimates from a Softmax neural network that
covers $K$ output logits using cross-entropy loss. Let us denote the
parametrized density

$q_{\theta}(\cdot\mid x)=\mathrm{softmax}(f_{\theta}(x)),\text{ where
}\mathrm{softmax}(v)_{j}=\frac{\exp(v_{j})}{\sum_{k=1}^{K}\exp(v_{k})},$

as the outputted discrete probability distribution over the labels of the
input $x$. Traditionally, we fit our neural network by minimizing the cross-
entropy loss on the training set:

$\hat{\theta}\in\underset{\theta}{\operatorname*{arg\,min}}\,\sum_{i=1}^{n}\operatorname{KL}(\delta_{y_{i}}\mid\mid
q_{\theta}(\cdot\mid x_{i}))\text{.}$

Here, $\delta_{y_{i}}$ is the Dirac Distribution with all of the
probabilitistic mass on $y_{i}$. Let us assume that we have trained and
acquired a $\hat{\theta}$ that has minimized the traditional cross-entropy
loss function on the training dataset. A natural conformity score is simply
the probability of a label according to the learned conditional distribution,
i.e., $\sigma(x,y)=q_{\hat{\theta}}(y\mid x)$. This approach is both
straightforward and efficient. The neural network’s flexibility allows it to
learn numerous label distributions with adaptivity across examples with less
explicit design of specific prior structure.

###  3.2 Regression to Classification Approach

Naturally, we strive to embody such practicality and effectiveness in
regression settings. However, the distribution of labels in the regression
scenario is continuous, and learning a continuous distribution directly using
a neural network is challenging (Rothfuss et al., 2019). Often, Bayesian or
kernel density estimators are employed to estimate this distribution. Other
techniques acquire knowledge of this distribution by training numerous
regressors and categorizing the regressors, for example, Conformal Histogram
Regression (Sesia & Romano, 2021).

However, these methods look different from the conformal prediction approaches
for classification. It would be desirable to be able to use similar methods
for both classification and regression conformal prediction. One method of
unifying classification and regression problems outside the conformal
prediction literature is known as Regression-as-Classification. We simply turn
a regression problem into a classification problem by binning the range space.
Specifically, we generate $K$ bins with $K$ equally spaced numbers covering
the interval $\mathcal{Y}=[y_{\min},y_{\max}]$, where $y_{\min}$ (or
$y_{\max}$) is the minimum (or maximum) value of the labels observed in the
training set. More explicitly, we define our discretization of the label space
as

$\hat{\mathcal{Y}}=\\{\hat{y}_{1},\dots,\hat{y}_{K}\\}\text{ where
}\hat{y}_{k+1}=\hat{y}_{k}+\frac{\hat{y}_{K}-\hat{y}_{1}}{K-1}\text{ with
}\hat{y}_{1}=y_{\min}\text{ and }\hat{y}_{K}=y_{\max}.$

These values $\hat{y}\in\hat{\mathcal{Y}}$ form the midpoints for each bin of
our discretization. Naturally, the $k$th bin is all the labels in the range
space $\mathcal{Y}$ closest to $\hat{y}_{k}$. Intuitively, we can treat each
bin as a class. Thus, we have turned a regression problem into a
classification problem. This method is simple but has yielded surprising
results. Some work has suggested that this form of binning results in more
stable training (Stewart et al., 2023) and gives significantly better results
for learning conditional expectations. To unify classification and regression
conformal prediction, a simple solution is to employ the Classification
Conformal Prediction model with discrete labels
$\tilde{y}_{i}=\operatorname*{arg\,min}_{\hat{y}\in\hat{\mathcal{Y}}}|y_{i}-\hat{y}|$.
This will aid in training the neural network with modified labels through
cross-entropy loss, resulting in a discrete distribution of
$q_{\theta}(\cdot\mid x)$, as outlined in the previous section.

To compute conformity scores for all labels, we employ linear interpolation
from the discrete probability function $q_{\theta}(\cdot\mid x)$ to generate
the continuous distribution $\bar{q}_{\theta}(\cdot\mid x)$. Nevertheless,
this approach encounters a critical issue when employed for regression.

###  3.3 Data Fitting

A critical problem with employing CrossEntropy loss in the classification
conformal prediction context is that any structural relationships between
classes are disregarded. This is not surprising given that in the
classification context, no structure exists between classes, and each class is
independent. Therefore, the CrossEntropy loss does not need to differentiate
between whether $q_{\theta}$ allocates probable mass far away or close to the
actual label. Instead, it only incentivizes the allocation of probable mass on
the correct label. However, within the regression setting, despite a multitude
of labels, the labels adhere to an ordinal structure. Thus, to enhance the
accuracy of labeling, it is imperative to devise a loss function that
incentivizes the allocation of probabilistic mass not only to the correct bin
but also to the neighboring bins. Formally, given an input and output pair
$(x,y)$, our goal is to determine a density estimate $q_{\theta}(\cdot\mid x)$
that assigns low (resp. high) probability to points that are far (resp. close)
to the true label $y$, i.e.,

$q_{\theta}(\hat{y}\mid x_{i})\text{ high (resp. small) when the loss
}\ell(\hat{y},y_{i})\text{ is small (resp. high).}$

Hence, a natural desideratum for learning the probability density function
$q_{\theta}$ is that their product $\ell(y,\hat{y})q(\hat{y}\mid x)$ is small
in expectation. We propose to find a distribution $q_{\theta}$ minimizing the
loss

$\displaystyle\mathbb{E}_{\hat{y}\sim q(\cdot\mid x)}[\ell(y,\hat{y})]$
$\displaystyle=\sum_{k=1}^{K}\ell(y,\hat{y}_{k})q(\hat{y}_{k}\mid x)$

Minimizing this loss in the space of all possible distributions
$q(\cdot\,|\,x)$ is equivalent to minimizing the original loss
$\ell(y,\cdot)$, where the minimizing distribution is a Dirac delta
$\delta_{\hat{y}_{\star}}$ at the minimizer $\hat{y}_{\star}$ of
$\ell(y,\cdot)$. Therefore, we expect an unregularized version of this loss to
share similarities with the typical empirical risk minimization on
$\ell(y,\cdot)$. However, a key difference is that when minimizing in a
restricted family of distributions (for example, those representable by a
neural network with a fixed architecture), the distributional output can
represent multi-modal or heavy-tailed label distributions. Minimizing the
original loss $\ell(y,\cdot)$, one would always be confined to a point
estimate.

###  3.4 Entropy Regularization

Although the proposed loss function better encodes the connection between
bins, it tends towards outputting Dirac distributions, as outlined in the
previous paragraph. Overconfidence in our neural network is a commonly
reported problem in the literature (Wei et al., 2022). Nevertheless,
smoothness has been a traditional requirement in density learning. As such, we
rely on a classical entropy-regularization technique for learning density
estimators (Wainwright & Jordan, 2008). Given the set of density estimators
that match the training label distribution well and put high probability mass
on the best bins, we prefer the probability distribution that maximizes the
entropy since this intuitively takes fewer assumptions on the data
distribution structure. Our choice is based on selecting density estimators
that effectively match the distribution of the training labels and assign a
higher probability to the best bins. Formally, we can calculate the entropy of
our probability distribution by using the Shannon entropy $\mathcal{H}$ of the
produced probability distribution $q(\cdot\,|\,x)$ as a penalty term as
follows:

$\mathcal{H}(q_{\theta}(\cdot\,|\,x))=\sum_{k=1}^{K}q_{\theta}(\hat{y}_{k}|x)\log
q_{\theta}(\hat{y}_{k}|x).$

#### Summary

We learn a distribution by minimizing the following expected loss over a given
training data $\mathcal{D}_{\rm tr}=\\{(x_{1},y_{1}),\ldots,(x_{n_{\rm
tr}},y_{n_{\rm tr}})\\}$ of size $n_{\rm tr}$:

$\mathcal{L}(\theta)=\sum_{i=1}^{{n_{\rm
tr}}}\sum_{k=1}^{K}\ell(y_{i},\hat{y}_{k})\,q_{\theta}(\hat{y}_{k}\,|\,x_{i})-\tau\mathcal{H}(q_{\theta}(\cdot\,|\,x_{i})),$
(2)

In particular, we choose $\ell(y_{i},\hat{y}_{k})$ as
$\ell(y_{i},\hat{y}_{k})=|y_{i}-\hat{y}_{k}|^{p}$ where $p>0$ is a
hyperparameter. This selection functions as a natural distance metric that
meets all required objectives. As a result, we can employ a technique similar
to the aforementioned Classification Conformal Prediction methods. Initially,
we train a Softmax neural network $f_{\theta}$ with $K$ logits, grounded on
the loss function Equation 2 on the training dataset. To calculate conformity
scores for the calibration set, we utilize the linearly interpolated
$\sigma(x,y)=\bar{q}_{\theta}(y\mid x)$ for a specific data point $(x,y)$.
Namely, for any $y$ between $\hat{y}_{k}$ and $\hat{y}_{k+1}$,
$\bar{q}_{\theta}(y\mid x)$ is defined as

$\displaystyle\bar{q}_{\theta}(y\mid x)$
$\displaystyle=\gamma_{k}q_{\theta}(\hat{y}_{k}\mid
x)+(1-\gamma_{k})q_{\theta}(\hat{y}_{k+1}\mid x)$ $\displaystyle\text{ where
}\gamma_{k}$
$\displaystyle=\frac{\hat{y}_{k+1}-y}{\hat{y}_{k+1}-\hat{y}_{k}}.$

Complete details of this procedure are outlined in Algorithm 1.

Algorithm 1 Regression to Classication Conformal Prediction (R2CCP).

1: Input:

  * •

Dataset $\mathcal{D}_{n}=\\{(x_{1},y_{1}),\ldots,(x_{n},y_{n})\\}$ and new
input $x_{n+1}$

  * •

Desired confidence level $\alpha\in(0,1)$

2: Hyperparameters: temperature $\tau>0$, $p>0$, number of bins $K>1$

3: Discretize the output space $[y_{\text{min}},y_{\text{max}}]$ into $K$
equidistant bins with midpoints $\\{\hat{y}_{1},\ldots,\hat{y}_{K}\\}$

4: Randomly split the dataset $\mathcal{D}_{n}$ in training $\mathcal{D}_{\rm
tr}$ and calibration $\mathcal{D}_{\rm cal}$

5: Find a distribution $q_{\hat{\theta}}(\cdot\,|\,x)$ by (approximately)
optimizing on the training set $\mathcal{D}_{\rm tr}$
$\hat{\theta}\in\arg\min_{\theta\in\mathbb{R}^{d}}~{}\sum_{i=1}^{{n_{\rm
tr}}}\sum_{k=1}^{K}|y_{i}-\hat{y}_{k}|^{p}q_{\theta}(\hat{y}_{k}\mid
x_{i})-\tau\mathcal{H}(q_{\theta}(\cdot\,|\,x_{i}))$ where
$q_{\theta}(\hat{y}_{k}|\,x)=\text{softmax}(f_{\theta}(x))_{k}$ for a model
(e.g., neural net) $f_{\theta}:\mathbb{R}^{d}\to\mathbb{R}^{K}$.

6: $\mathcal{S}\leftarrow\left\\{\bar{q}_{\hat{\theta}}(y\mid x)~{}\text{ for
}~{}(x,y)\in\mathcal{D}_{\rm cal}\right\\}$ #
COMMENT$\bar{q}_{\hat{\theta}}(\cdot\,|\,x)$ is linear interpolation of
softmax probabilities.

7: $Q_{\alpha}(\mathcal{D}_{\rm
cal})\leftarrow\operatorname{quantile}(\mathcal{S},\alpha)$

8: return Conformal Set
$\Gamma^{(\alpha)}(x_{n+1})=\\{z\in\mathbb{R}~{}|~{}\bar{q}_{\hat{\theta}}(z\mid
x_{n+1})\geq Q_{\alpha}(\mathcal{D}_{\rm cal})\\}$

##  4 Experiments

All code to run our method can be installed via

pip install r2ccp.

We investigate the empirical behavior of our R2CCP (Regression-to-
Classification Conformal Prediction) method, which we have explained in detail
in Algorithm 1. We have three sets of experiments. The first one is described
in Section 4.1 and presents empirical evidence of the algorithm’s ability to
produce narrower intervals by utilizing label density characteristics,
including heteroscedasticity, bimodality, or a combination. Section 4.2
demonstrates the effectiveness of our algorithm on synthetic and real data by
comparing it with various benchmarks from the Conformal Prediction literature
in terms of length and coverage. Section 4.3 evaluates the effect of different
loss functions on the final learned densities and their impact on the
intervals produced. All experiments were run over $5$ different seeds at a
coverage level of $90\%$, and the standard error over the experiments is
reported in the subscript. We do not tune the hyperparameters and keep values
of $K=50$, $p=0.5$, and $\tau=0.2$ constant across all experiments. For all
experiments, we report length, meaning the length of all the sets predicted,
and coverage, the percent of instances where the true label is contained in
the predicted intervals.

###  4.1 Specific Characteristics of Label Distribution

#### Heteroscedascity

We generate a toy dataset where the input is one-dimensional. It contains
samples from the following distribution: $y\sim\mathcal{N}(0,|x|)$. The label
distribution is heteroscedastic, meaning the variance of the labels changes as
the input changes. In traditional Conformal Prediction literature, many
existing algorithms fail to capture heteroscedasticity, resulting in wide
intervals for inputs $x$ where the label distribution of $y$ has low variance.
However, our learned algorithm can directly learn this relation and adjust the
outputted probability distribution accordingly. Thus, we see that the lengths
of the intervals will increase as the variance of the label distribution
increases, which is desirable. We see this relation in Figure 1a. Moreover, we
also use the dataset generated from (Lei & Wasserman, 2014) as discussed in
Appendix D. This dataset exhibits heteroskedascity and bimodality as $X$
passes $-0.5$. We see that our learned method can adjust intervals accordingly
to maintain coverage and length for all $X$. We plot this in Figure 1b. We
plot how the intervals (grey) change as the data distribution (black) changes.
As the variance of the labels increases as $x$ increases, the produced
intervals adaptively get wider, taking advantage of the heteroscedasticity.

#### Bimodality

We showcase our algorithm’s capability to address labels with a bimodal
distribution. Our bimodal dataset is generated by repeatedly (1) sampling two
sets of random features that are close geometrically and (2) giving one set of
features a label of $1$ plus some Gaussian noise and giving the other a label
of $-1$ plus some Gaussian noise. Therefore, our dataset is comprised of many
similar data points with bimodally distributed labels. This bimodal
distribution is particularly hard for many existing CP algorithms to solve
since it requires outputting two disjoint conformal sets to achieve low
length. CQR, for instance, cannot deal with this circumstance and will
generate a conformal set that covers the entire range space. However, our
method is flexible enough to assign a low probability to labels between $-1$
and $1$, and our resulting conformal set will not include these intermediate
labels. We see this in Figure 6d. In Figure 6d, our outputted probability
distribution has two modes around labels $-1$ and $1$ and assigns a low
probability value to the valley between the modes.

###  4.2 Comparison to other Conformal Prediction Algorithms

The crucial criteria for assessing a Conformal Prediction algorithm consist of
(1) coverage, representing the percentage of generated conformal sets that
incorporate accurate labels, and (2) the length of the generated conformal
sets. Our baseline techniques consist of the Kernel Density Estimator (KDE) as
proposed by Lei & Wasserman (2014), alongside a conformity score shown by the
estimated probability density. Furthermore, we take into consideration Fong &
Holmes (2021) (CB), which employs the likelihood of a posterior distribution
in their conformity function, and Sesia & Romano (2021) (CHR), which uses
quantile regressors on every bin of a histogram density estimator. We have
included the Conformal Quantile Regression as described by Romano et al.
(2019) (CQR), which employs conformity based on the labels’ distance from
quantile regressors. Moreover, we had the Distributional Conformal Prediction
which uses a Cumulative Distribution Function to form its intervals from
Chernozhukov et al. (2021) (DCP). Additionally, we have used the Lasso
Conformal Predictor with a distance to mean regressors, which is the most
straightforward option (LASSO). We use synthetic data exhibiting bimodally
distributed and log-normally distributed to illustrate particular weaknesses
of existing methods. We use real datasets also in Romano et al. (2019).
Specifically, these are several datasets from the UCI Machine Learning
repository (Bio, Blog, Concrete, Community, Energy, Forest, Stock, Cancer,
Solar, Parkinsons, Pendulum) (Nottingham et al., 2023) and the medical
expenditure panel survey number 19–21 (MEPS-19–21) (Cohen et al., 2009). These
regression datasets are commonly used to benchmark regression models. Our
approach yields tighter intervals on real datasets than some of the strongest
baselines.

#### Results

We report lengths and coverages results in Table 1 respectively. We added
figures depicting example probability distributions on these datasets in
Figure 6 in Appendix H. The intervals produced by our method are the shortest
over $10$ of the $16$ datasets. From Figure 6, we see that our method can
learn many different shaped distributions well, which accounts for this
significant improvement in intervals. Overall, the Kernel Density Estimation
and our method can accurately predict the best intervals on datasets where the
connection between the data and feature is simple, such as Bimodal or Log-
Normal. While the Kernel Density Estimator can fail when the labels are
complexly related to inputs, no such connection exists in these datasets since
the labels are independent of the features. Thus, the Kernel Density Method’s
simplicity allows it to learn the label density quickly. Our method also seems
to handle the case where there is no connection between the data and the
labels, as seen in Figure 6d in Appendix H. Moreover, on datasets where the
label density is smooth and close to the Laplace prior, such as on Figure 6c
and Figure 6i, Conformal Bayes and our method can accurately learn this
distribution since both methods can output smooth distributions. Moreover, on
very sharp datasets such as on Concrete in Figure 6j and Energy in Figure 6h,
both our method and CQR can capture the sharp and unnoisy distribution needed
to achieve strong length. Moreover, on complex distributions such as in
Pendulum in Figure 6e and Bio in Figure 6f, we see that both CHR and our
method have the flexibility to portray complex distributions resulting in the
most accurate intervals. Thus, the flexibility of our algorithm to smoothly
learn sharp, wide, complex, and simple conditional label densities results in
our method achieving the best length most consistently over the entire
dataset.

Dataset |  Bimodal |  LogNorm |  Concrete |  MEPS-19 |  MEPS-20 |  MEPS-21 |  Bio |  Community
---|---|---|---|---|---|---|---|---
CQR |  $2.14_{(0.01)}$ |  $1.58_{(0.03)}$ |  $\mathbf{0.39_{(0.01)}}$ |  $1.87_{(0.05)}$ |  $2.00_{(0.07)}$ |  $1.99_{(0.03)}$ |  $1.34_{(0.01)}$ |  $\mathbf{1.44_{(0.03)}}$
KDE |  $\mathbf{0.35_{(0.01)}}$ |  $\mathbf{1.40_{(0.04)}}$ |  $1.54_{(0.03)}$ |  $2.16_{(0.02)}$ |  $2.51_{(0.05)}$ |  $2.39_{(0.03)}$ |  $2.27_{(0.00)}$ |  $2.23_{(0.09)}$
Lasso |  $2.14_{(0.00)}$ |  $3.30_{(0.06)}$ |  $2.74_{(0.03)}$ |  $4.64_{(0.06)}$ |  $4.79_{(0.04)}$ |  $4.92_{(0.04)}$ |  $3.89_{(0.00)}$ |  $3.26_{(0.05)}$
CB |  $2.16_{(0.00)}$ |  $1.45_{(0.01)}$ |  $0.96_{(0.01)}$ |  $4.47_{(0.02)}$ |  $4.50_{(0.02)}$ |  $4.51_{(0.01)}$ |  $2.09_{(0.00)}$ |  $1.80_{(0.00)}$
CHR |  $2.14_{(0.00)}$ |  $1.52_{(0.05)}$ |  $0.47_{(0.02)}$ |  $2.60_{(0.08)}$ |  $2.53_{(0.03)}$ |  $2.75_{(0.03)}$ |  $1.59_{(0.01)}$ |  $\mathbf{1.49_{(0.05)}}$
DCP |  $2.14_{(0.01)}$ |  $1.74_{(0.07)}$ |  $0.47_{(0.01)}$ |  $68.64_{(0.15)}$ |  $66.71_{(0.40)}$ |  $67.56_{(0.33)}$ |  $1.74_{(0.01)}$ |  $1.59_{(0.03)}$
R2CCP (Ours) |  $0.46_{(0.01)}$ |  $1.96_{(0.03)}$ |  $\mathbf{0.38_{(0.01)}}$ |  $\mathbf{1.60_{(0.01)}}$ |  $\mathbf{1.70_{(0.03)}}$ |  $\mathbf{1.72_{(0.03)}}$ |  $\mathbf{1.11_{(0.01)}}$ |  $\mathbf{1.47_{(0.03)}}$

Dataset |  Diabetes |  Solar |  Parkinsons |  Stock |  Cancer |  Pendulum |  Energy |  Forest
---|---|---|---|---|---|---|---|---
CQR |  $1.30_{(0.05)}$ |  $1.98_{(0.22)}$ |  $\mathbf{0.42_{(0.01)}}$ |  $1.85_{(0.23)}$ |  $\mathbf{3.09_{(0.13)}}$ |  $2.25_{(0.31)}$ |  $\mathbf{0.19_{(0.01)}}$ |  $3.18_{(0.19)}$
KDE |  $1.34_{(0.05)}$ |  $\mathbf{0.50_{(0.01)}}$ |  $3.79_{(0.02)}$ |  $4.72_{(0.29)}$ |  $3.82_{(0.09)}$ |  $3.96_{(0.09)}$ |  $2.72_{(0.07)}$ |  $\mathbf{2.90_{(0.17)}}$
Lasso |  $3.01_{(0.04)}$ |  $3.54_{(0.12)}$ |  $3.46_{(0.03)}$ |  $1.39_{(0.04)}$ |  $3.55_{(0.14)}$ |  $3.99_{(0.07)}$ |  $1.29_{(0.03)}$ |  $3.97_{(0.29)}$
CB |  $\mathbf{1.19_{(0.01)}}$ |  $3.78_{(0.02)}$ |  $3.42_{(0.01)}$ |  $1.32_{(0.01)}$ |  $\mathbf{3.14_{(0.04)}}$ |  $3.71_{(0.03)}$ |  $1.26_{(0.01)}$ |  $3.75_{(0.03)}$
CHR |  $1.40_{(0.02)}$ |  $1.49_{(0.23)}$ |  $0.68_{(0.02)}$ |  $1.59_{(0.07)}$ |  $3.42_{(0.12)}$ |  $\mathbf{1.69_{(0.11)}}$ |  $0.23_{(0.01)}$ |  $\mathbf{3.03_{(0.15)}}$
DCP |  $1.29_{(0.02)}$ |  $15.69_{(0.06)}$ |  $0.83_{(0.04)}$ |  $1.69_{(0.10)}$ |  $3.57_{(0.07)}$ |  $1.76_{(0.10)}$ |  $0.23_{(0.01)}$ |  $6.00_{(0.02)}$
R2CCP (Ours) |  $1.34_{(0.02)}$ |  $3.80_{(2.61)}$ |  $0.50_{(0.00)}$ |  $\mathbf{0.92_{(0.02)}}$ |  $\mathbf{3.21_{(0.08)}}$ |  $\mathbf{1.60_{(0.07)}}$ |  $\mathbf{0.20_{(0.02)}}$ |  $3.80_{(0.26)}$

Table 1: This is the length results over all datasets. We see that our method
achieves the best length on $10$ of the $16$ datasets. Meanwhile, CQR is best
at $5$, CHR is best at $3$, CB is best at $1$, and KDE is the best at $3$. Our
method achieves the shortest intervals across these datasets.

(a) Original (b) No Entropy

(c) MLE (d) MLE with Entropy

Figure 2: The resulting density estimates with different loss functions. We
see that removing entropy from our loss function or using MLE as error terms
causes sharp density estimates. Moreover, adding in the entropy regularization
with MLE does not smooth the density estimate but instead raises the entire
distribution uniformly; which does not provide valuable information for CP.

###  4.3 Ablation Studies

Our loss function consists of an error and entropy terms. The error term
penalizes distributions that put weight far away from the true label, whereas
the entropy term acts as a regularizing term in the probability distribution
space. We will do ablations on both the error and entropy terms to illustrate
the importance of each part. For the error term, there are several notable
alternatives that a practitioner may use. An alternative is the log maximum
likelihood or cross-entropy formulations of the error term, which we denote as

$\mathcal{L}_{\rm
MLE}(\theta)=-\sum_{i=1}^{N}\log(q_{\theta}(\tilde{y}_{i}\mid x_{i})),$

where we remind the reader that
$\tilde{y}_{i}=\operatorname*{arg\,min}_{\hat{y}\in\hat{\mathcal{Y}}}|\hat{y}-y_{i}|$.
We note that both MLE and CE are equivalent formulations in this setting.
These two are standard error terms often used in practice. We will train our
models with this loss function over all datasets to see how the change of
error term affects the intervals’ length and the learned density. We will show
that our chosen error term is better for producing optimal intervals
empirically. Moreover, we will also demonstrate the importance of our entropy
term. We do this by retraining our models with the entropy part omitted. We
also test by combining the MLE error term with entropy regularization.
Therefore, we will perform $4$ different ablation experiments on the loss
function by retraining on different variations of the loss functions and
reporting the final length and coverage generated. Overall, the four loss
functions we will use are the original loss function $\mathcal{L}$, the
original loss function without entropy $\mathcal{L}_{\text{\rm NE}}$, the MLE
loss function $\mathcal{L}_{\text{\rm MLE}}$, and the MLE loss function with
entropy added $\mathcal{L}_{\text{\rm MLE + E}}$. We will demonstrate that our
chosen loss function delivers the best results across all the loss functions.
We report our results in Table 2.

#### Results

We see that our chosen loss function achieves the best length across most loss
functions in Table 2. The only datasets where our chosen loss function does
not achieve the best length are Energy, Solar, Bimodal, and Log-Normal. We now
discuss the individual differences between our loss function and each variant
loss function. When comparing our loss $\mathcal{L}$ with the variant without
entropy $\mathcal{L}_{\text{\rm NE}}$, we see the lengths constantly increase
except for the bimodal distribution. Without an entropy regularizing term, the
outputted probability distributions are not smooth, placing much mass on a
single data point. This overconfidence results in the $\alpha$th quantile of
the low probability of the true label as in Figure 2.

When looking at the differences between our loss function when changing the
error term to Maximum Log Likelihood, we see similar overconfidence. While MLE
loss works on several datasets, such as Energy, Solar, Bimodal, and Log-
Normal, it also performs poorly on other datasets. MLE works well for datasets
with less simple label distributions but fails otherwise due to similar
overconfidence. Even when adding entropy as a regularizer to the MLE loss as
in $\mathcal{L}_{\text{\rm MLE + E}}$, we see that the addition of entropy
does not improve the length. Since MLE loss does not treat nearby bins
differently than far away bins, regularization decreases the overconfidence
but does so uniformly across all bins. The result is a roughly uniform
distribution with a single spike. This does not utilize the structure of
regression where bins near the best bin are preferable. Thus, prioritizing
nearby bins is crucial for achieving strong length. Moreover, when adding
entropy regularization to our error term, we do not see a uniform increase in
probability across all bins but instead a smoothening of a direct. This
connection between entropy regularization and our distance-based loss function
appears to be a powerful synergy per the intervals produced. Thus, the error
term utilizing the regression structure by prioritizing nearby bins and the
regularizing entropy term preventing overconfidence seem crucial for producing
good intervals.

Another important insight is that the datasets where other loss functions
excel are the same ones where other CP methods perform better than ours. In
particular, on the Bimodal, Log-Normal, Parkinson’s, and Solar datasets, our
CP method is not superior, and other loss functions outperform it. This
indicates that the MLE’s sharp distributions result in greater accuracy on
these specific datasets. Hence, adjusting the weighting terms for entropy
$\tau$ and smoothing $p$ could enhance the smoothness of the learned
distribution and, thereby, the performance of our algorithm across these
datasets, leading to the best outcomes. Nevertheless, to maintain fairness, we
avoid this approach.

Dataset |  Bimodal |  Log-Normal |  Concrete |  MEPS-19 |  MEPS-20 |  MEPS-21 |  Bio |  Community
---|---|---|---|---|---|---|---|---
$\mathcal{L}_{\text{NE}}$ |  $0.43_{0.001}$ |  $3.25_{0.058}$ |  $1.89_{0.037}$ |  $3.73_{0.057}$ |  $3.79_{0.061}$ |  $8.78_{0.163}$ |  $2.01_{0.006}$ |  $3.64_{0.054}$
$\mathcal{L}_{\text{MLE}}$ |  $\mathbf{0.35_{0.001}}$ |  $\mathbf{1.76_{0.020}}$ |  $0.74_{0.016}$ |  $\mathbf{1.58_{0.008}}$ |  $\mathbf{1.59_{0.009}}$ |  $1.71_{0.034}$ |  $2.71_{0.001}$ |  $2.32_{0.047}$
$\mathcal{L}_{\text{MLE + E}}$ |  $0.36_{0.001}$ |  $3.52_{0.008}$ |  $1.91_{0.017}$ |  $\mathbf{1.60_{0.008}}$ |  $1.63_{0.027}$ |  $\mathbf{1.70_{0.008}}$ |  $2.32_{0.003}$ |  $3.77_{0.013}$
$\mathcal{L}$ |  $0.44_{0.002}$ |  $1.82_{0.034}$ |  $\mathbf{0.37_{0.004}}$ |  $\mathbf{1.60_{0.009}}$ |  $\mathbf{1.59_{0.011}}$ |  $\mathbf{1.69_{0.013}}$ |  $\mathbf{1.10_{0.004}}$ |  $\mathbf{1.50_{0.025}}$

Dataset |  Diabetes |  Solar |  Parkinsons |  Stock |  Cancer |  Pendulum |  Energy |  Forest
---|---|---|---|---|---|---|---|---
$\mathcal{L}_{\text{NE}}$ |  $1.92_{0.052}$ |  $2.20_{0.233}$ |  $4.76_{0.033}$ |  $10.17_{0.165}$ |  $3.26_{0.207}$ |  $12.82_{0.479}$ |  $3.30_{0.093}$ |  $3.33_{0.254}$
$\mathcal{L}_{\text{MLE}}$ |  $1.56_{0.031}$ |  $\mathbf{0.14_{0.005}}$ |  $0.34_{0.006}$ |  $4.48_{0.214}$ |  $3.26_{0.117}$ |  $3.00_{0.203}$ |  $0.25_{0.011}$ |  $3.04_{0.110}$
$\mathcal{L}_{\text{MLE + E}}$ |  $1.96_{0.012}$ |  $0.15_{0.005}$ |  $\mathbf{0.23_{0.001}}$ |  $9.73_{0.074}$ |  $3.95_{0.022}$ |  $13.40_{0.221}$ |  $0.25_{0.009}$ |  $5.24_{0.055}$
$\mathcal{L}$ |  $\mathbf{1.37_{0.037}}$ |  $0.66_{0.092}$ |  $0.46_{0.036}$ |  $\mathbf{1.95_{0.039}}$ |  $\mathbf{3.04_{0.142}}$ |  $\mathbf{1.62_{0.042}}$ |  $\mathbf{0.21_{0.025}}$ |  $\mathbf{2.92_{0.127}}$

Table 2:  We present the length results over all of the variant loss
functions. We find that our loss function delivers the best over $12$
datasets, demonstrating that our chosen loss function often generates the best
intervals. For datasets where our method does not deliver the best results, it
is likely that tuning the weight on the entropy $\tau$ and the smoothing term
$p$ would likely have improved the results, but we do not do this for the sake
of evaluation.

## References

  * Chernozhukov et al. (2021) Victor Chernozhukov, Kaspar Wüthrich, and Yinchu Zhu.  Distributional conformal prediction.  _Proceedings of the National Academy of Sciences_ , 118(48):e2107794118, 2021.
  * Cohen et al. (2009) Joel W Cohen, Steven B Cohen, and Jessica S Banthin.  The medical expenditure panel survey: a national information resource to support healthcare cost research and inform policy and practice.  _Medical Care_ , pp. 44–50, 2009.
  * Cortés-Ciriano & Bender (2020) Isidro Cortés-Ciriano and Andreas Bender.  Concepts and applications of conformal prediction in computational drug discovery.  _Artificial Intelligence in Drug Discovery_ , pp. 63–101, 2020\.
  * Diaz & Marathe (2019) Raul Diaz and Amit Marathe.  Soft labels for ordinal regression.  In _IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2019.
  * Fong & Holmes (2021) Edwin Fong and Chris C Holmes.  Conformal Bayesian computation.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , 2021.
  * Fu et al. (2018) Huan Fu, Mingming Gong, Chaohui Wang, Kayhan Batmanghelich, and Dacheng Tao.  Deep ordinal regression network for monocular depth estimation.  In _IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2018.
  * Ghosh et al. (2023) Subhankar Ghosh, Taha Belkhouja, Yan Yan, and Janardhan Rao Doppa.  Improving uncertainty quantification of deep classifiers via neighborhood conformal prediction: Novel algorithm and theoretical analysis.  _arXiv preprint arXiv:2303.10694_ , 2023.
  * Guha et al. (2023) Etash Kumar Guha, Eugene Ndiaye, and Xiaoming Huo.  Conformalization of sparse generalized linear models.  In _International Conference on Machine Learning (ICML)_ , 2023.  URL https://proceedings.mlr.press/v202/guha23b.html.
  * Izbicki et al. (2022) Rafael Izbicki, Gilson Shimizu, and Rafael B Stern.  Cd-split and hpd-split: Efficient conformal regions in high dimensions.  _Journal of Machine Learning Research_ , 23(87):1–32, 2022.
  * Krizhevsky et al. (2012) Alex Krizhevsky, Ilya Sutskever, and Geoffrey E Hinton.  Imagenet classification with deep convolutional neural networks.  _Advances in neural information processing systems_ , 25, 2012.
  * Lambrou et al. (2009) Antonis Lambrou, Harris Papadopoulos, and Alex Gammerman.  Evolutionary conformal prediction for breast cancer diagnosis.  In _2009 9th international conference on information technology and applications in biomedicine_ , pp. 1–4. IEEE, 2009.
  * Lambrou et al. (2010) Antonis Lambrou, Harris Papadopoulos, Efthyvoulos Kyriacou, Constantinos S Pattichis, Marios S Pattichis, Alexander Gammerman, and Andrew Nicolaides.  Assessment of stroke risk based on morphological ultrasound image analysis with conformal prediction.  In _Artificial Intelligence Applications and Innovations: 6th IFIP WG 12.5 International Conference, AIAI 2010, Larnaca, Cyprus, October 6-7, 2010. Proceedings 6_ , pp. 146–153. Springer, 2010.
  * Lei & Wasserman (2014) Jing Lei and Larry Wasserman.  Distribution-free prediction bands for non-parametric regression.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 76(1):71–96, 2014.
  * Lei et al. (2013) Jing Lei, James Robins, and Larry Wasserman.  Distribution-free prediction sets.  _Journal of the American Statistical Association_ , 108(501):278–287, 2013.
  * Lei et al. (2018) Jing Lei, Max G’Sell, Alessandro Rinaldo, Ryan J Tibshirani, and Larry Wasserman.  Distribution-free predictive inference for regression.  _Journal of the American Statistical Association_ , 113(523):1094–1111, 2018.
  * Lin et al. (2021) Zhen Lin, Shubhendu Trivedi, and Jimeng Sun.  Locally valid and discriminative prediction intervals for deep learning models.  _Advances in Neural Information Processing Systems_ , 34:8378–8391, 2021.
  * Lu et al. (2022) Charles Lu, Anastasios N Angelopoulos, and Stuart Pomerantz.  Improving trustworthiness of ai disease severity rating in medical imaging with ordinal conformal prediction sets.  In _International Conference on Medical Image Computing and Computer-Assisted Intervention_ , pp. 545–554. Springer, 2022.
  * Ndiaye (2022) Eugene Ndiaye.  Stable conformal prediction sets.  In _International Conference on Machine Learning_ , pp. 16462–16479. PMLR, 2022.
  * Ndiaye & Takeuchi (2023) Eugene Ndiaye and Ichiro Takeuchi.  Root-finding approaches for computing conformal prediction set.  _Machine Learning_ , 112(1):151–176, 2023.
  * Nottingham et al. (2023) Kolby Nottingham, Rachel Longjohn, and Markelle Kelly.  UCI Machine Learning Repository, 2023.  URL https://archive.ics.uci.edu/datasets.  Accessed: September, 2023.
  * Papadopoulos et al. (2002) Harris Papadopoulos, Kostas Proedrou, Volodya Vovk, and Alex Gammerman.  Inductive confidence machines for regression.  In _Machine Learning: ECML 2002: 13th European Conference on Machine Learning Helsinki, Finland, August 19–23, 2002 Proceedings 13_ , pp. 345–356. Springer, 2002.
  * Romano et al. (2019) Y. Romano, E. Patterson, and E. J. Candes.  Conformalized quantile regression.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , 2019.
  * Rothe et al. (2015) Rasmus Rothe, Radu Timofte, and Luc Van Gool.  Dex: Deep expectation of apparent age from a single image.  In _IEEE International Conference on Computer Vision Workshops (ICCVW)_ , 2015.
  * Rothfuss et al. (2019) Jonas Rothfuss, Fabio Ferreira, Simon Walther, and Maxim Ulrich.  Conditional density estimation with neural networks: Best practices and benchmarks.  _arXiv preprint arXiv:1903.00954_ , 2019.
  * Sesia & Romano (2021) Matteo Sesia and Yaniv Romano.  Conformal prediction using conditional histograms.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , 2021.
  * Stewart et al. (2023) Lawrence Stewart, Francis Bach, Quentin Berthet, and Jean-Philippe Vert.  Regression as classification: Influence of task formulation on neural network features.  In _International Conference on Artificial Intelligence and Statistics (AISTATS)_ , 2023.
  * Van Den Oord et al. (2016) Aäron Van Den Oord, Nal Kalchbrenner, and Koray Kavukcuoglu.  Pixel recurrent neural networks.  In _International Conference on Machine Learning (ICML)_ , 2016.
  * Vovk et al. (2005) V. Vovk, A. Gammerman, and G. Shafer.  _Algorithmic learning in a random world_.  Springer, 2005.
  * Wainwright & Jordan (2008) Martin J Wainwright and Michael I Jordan.  Graphical models, exponential families, and variational inference.  _Foundations and Trends® in Machine Learning_ , 1(1–2):1–305, 2008.
  * Wei et al. (2022) Hongxin Wei, Renchunzi Xie, Hao Cheng, Lei Feng, Bo An, and Yixuan Li.  Mitigating neural network overconfidence with logit normalization.  In _International Conference on Machine Learning_ , pp. 23631–23644. PMLR, 2022.
  * Weigend & Srivastava (1995) Andreas S Weigend and Ashok N Srivastava.  Predicting conditional probability distributions: A connectionist approach.  _International Journal of Neural Systems_ , 6(02):109–118, 1995.
  * Xu et al. (2023) Yunpeng Xu, Wenge Guo, and Zhi Wei.  Conformal risk control for ordinal classification.  In _Uncertainty in Artificial Intelligence_ , pp. 2346–2355. PMLR, 2023.
  * Zhang et al. (2016) Richard Zhang, Phillip Isola, and Alexei A Efros.  Colorful image colorization.  In _European Conference on Computer Vision (ECCV)_ , 2016.

##  Appendix A Author Contribution Statement

Etash Guha, Shlok Natarajan, Thomas Möllenhoff, Eugene Ndiaye, and Mohammad
Emtiyaz Khan were all responsible for the development of the main idea and
writing of the paper. Etash Guha and Shlok Natarjan were responsible for the
experiments and implementation of the algorithm. Etash Guha, Thomas
Möllenhoff, and Eugene Ndiaye were responsible for the design of the loss
function, binning, and entropic regularization. Etash Guha and Eugene Ndiaye
were responsible for writing the related works.

##  Appendix B Limitations

There are several limitations to our algorithm. The loss function we choose to
optimize needs larger representational power due to increased complexity of
the output. Therefore, the neural networks we used in practice to minimize
this loss function are larger than the ones used for CQR. Moreover, to achieve
strong training, we had to slightly vary the hyperparameters for a dataset
depending on its size. For example, larger datasets needed models to be
trained with smaller initial learning rate to avoid divergence. This was not
needed for Quantie Regressors. If the network does not train well or achieve
good training or validation loss, the intervals produced will be suboptimal,
so finding a training setup that minimizes the loss effectively is crucial for
our algorithm to produce meaningful intervals. Moreover, we do not directly
learn the label distribution and our produced probability distributions do not
directly mirror that of the ground truth distribution. One reason for this is
the use of entropy regularization. While using entropy regularization
introduces bias, we found it to be necessary to prevent the neural network
from being overconfident on a simple bin in practice.

##  Appendix C More Details on Conformal Prediction and Variants

Given an input observation $x$ without its label $y$, Conformal Prediction
method Vovk et al. (2005) allows to select set of values that its unobserved
response $y$ are likely to take. For an arbitrary value $z$, the rational is
to score how the example $(x_{n+1},z)$ is a typical/conformal observation
given a sequence of iid data $\mathcal{D}_{n}=\\{(x_{i},y_{i})\\}_{i\in[n]}$.
Intuitively, the scoring function, that we denote $\sigma$, measures how
different is the new example compared to examples that we have already
observed. We compute the score function $\sigma(x,y)$ for every example in the
augmented dataset $\mathcal{D}_{n+1}(z)=\mathcal{D}_{n}\cup(x,z)$ and report
the fraction of time previous example score better than the new one.

This define the so called conformity function $\pi(\cdot)\in[0,1]$ as

$\displaystyle\pi(z)=1-\frac{\mathrm{Rank}(z)}{n+1}\text{ where
}\mathrm{Rank}(z)=\sum_{(x,y)\in
D_{n+1}(z)}\mathds{1}\\{\sigma(x,y)\leq\sigma(x_{n+1},z)\\}.$ (3)

If $\pi(z)$ is small (resp. large) then $z$ is non conformal (resp.
conformal). Therefore, for a threshold level $\alpha\in(0,1)$, a simple way to
find our conformal set $\Gamma^{(\alpha)}(x_{n+1})$ for the response $y_{n+1}$
of the input $x_{n+1}$ is to select all the values $z$ whose conformity exceed
the threshold i.e.

$\Gamma^{(\alpha)}(x_{n+1})=\\{z\in\mathbb{R}:\,\pi(z)\geq\alpha\\}.$

The validity of the method rely on the fact that if we assume that order in
which the observations arrive is irrelevant and that the score function
$\sigma$ is invariant wrt permutation of the data, then the sequence of scores
$\\{\sigma(x_{1},y_{1}),\ldots,\sigma(x_{n},y_{n}),\sigma(x_{n+1},y_{n+1})\\}$
are also equally likely and thus the $\mathrm{Rank}(y_{n+1})$ is uniformly
distributed in $\\{1,\ldots,n+1\\}$. As such, for any confidence level
$\alpha\in(0,1)$,

$\mathbb{P}(y_{n+1}\in\Gamma^{(\alpha)}(x_{n+1}))\geq\mathbb{P}(\pi(y_{n+1})\geq\alpha)=\mathbb{P}(\mathrm{Rank}(y_{n+1})\leq(n+1)(1-\alpha))\geq
1-\alpha.$

#### Conditional validity

The previous coverage guarantee may not always be ideal because it is not
conditional on the particular input $x_{n+1}$. Unfortunately, it is not
possible to have such conditional validity when the distribution of the inputs
is continuous. Nevertheless, it can be approximated using some localisation
strategies based on a weighted conformity function; see (Lin et al., 2021;
Ghosh et al., 2023; Izbicki et al., 2022). It is interesting to note that some
of these approaches make use of (clever) partitioning of the feature space
(which can be high dimensional) whereas our approach partitions the output
space. Both strategies can be combined. This also might allows to further
reduce the size of the conformal prediction sets.

Validity can be conditional on a discrete attribute of the input, $x_{n+1}$,
through the use of Mondrian Conformal Prediction. This method involves a
slight modification of the conformity function. Let a taxonomy $\kappa$ be a
function that partitions the data into a finite number of categories eg
$\kappa(x,.)$ is a classification of the example $x$. This allows to modify
the pi-value function in order to achieve validity conditional on the category
of the test example :

$\pi(z)=\sum_{(x,y)\in\mathcal{D}_{n+1}(z)}\frac{\mathds{1}\\{\sigma(x,y)\leq\sigma(x_{n+1},z),\kappa(x,y)=\kappa(x_{n+1},z)\\}}{\sum_{(x,y)\in\mathcal{D}_{n+1}(z)}\mathds{1}\\{\kappa(x,y)=\kappa(x_{n+1},z)\\}}$

#### Computational issues

As mentioned in the introduction, the selection of the score function affects
the size and shape of the prediction sets. There are also computational issues
that arise when the full data set is used to fit the function. For example,
consider using the estimated distance to the conditional mean
$\sigma(x,y)=|y-\hat{E}[y|x]|$. The estimator $\hat{E}[\cdot|x]$ is learned on
the extended dataset $\mathcal{D}_{n+1}(z)$ for all possible values $z$, which
is computationally infeasible for most methods without stronger assumptions;
see (Ndiaye & Takeuchi, 2023), (Ndiaye, 2022). As such, most of the popular
strategies are based on the fitting of the prediction model to an left-out
independent set of data, as we did in this paper.

##  Appendix D Details on Dataset from Lei & Wasserman (2014)

For the dataset referenced in Figure 1b, we generate the dataset by sampling
many $(X,Y)$ pairs from the following distribution.

$\displaystyle X$ $\displaystyle\sim\operatorname{Uniform}([-1.5,1.5])$
$\displaystyle Y\mid X$
$\displaystyle\sim\frac{1}{2}\mathcal{N}\left(f(X)-g(X),\sigma^{2}(X)\right)+\frac{1}{2}\mathcal{N}\left(f(X)+g(X),\sigma^{2}(X)\right)$
$\displaystyle f(X)$ $\displaystyle=(X-1)^{2}(X+1),$ $\displaystyle g(X)$
$\displaystyle=4\sqrt{(X+1/2)\mathds{1}_{\\{X\geq-1/2\\}}},$
$\displaystyle\sigma^{2}(X)$ $\displaystyle=\frac{1}{4}+|X|$

This dataset demonstrates bimodality after $x$ passes the threshold of $-0.5$.
This dataset was similarly used by the Lei & Wasserman (2014).

##  Appendix E Coverage Data

We have presented the coverage data for the Conformal Prediction comparison as
well as the ablation experiments in Table 3 and Table 4, respectively. Since
all methods obey the classical Conformal Prediction framework, we expect
coverage at the guaranteed level of $1-\alpha$. This guarantee is indeed what
we see across all experiments. This confirms that our coverage guarantee
indeed holds in practice.

Dataset |  Bimodal |  LogNorm |  Concrete |  MEPS-19 |  MEPS-20 |  MEPS-21 |  Bio |  Community
---|---|---|---|---|---|---|---|---
CQR |  $0.90_{(0.00)}$ |  $0.91_{(0.00)}$ |  $0.91_{(0.00)}$ |  $0.90_{(0.00)}$ |  $0.90_{(0.00)}$ |  $0.90_{(0.00)}$ |  $0.90_{(0.00)}$ |  $0.90_{(0.00)}$
KDE |  $0.93_{(0.00)}$ |  $0.90_{(0.00)}$ |  $0.90_{(0.00)}$ |  $0.91_{(0.00)}$ |  $0.92_{(0.00)}$ |  $0.91_{(0.00)}$ |  $0.90_{(0.00)}$ |  $0.90_{(0.00)}$
Lasso |  $0.90_{(0.00)}$ |  $0.91_{(0.00)}$ |  $0.90_{(0.00)}$ |  $0.90_{(0.00)}$ |  $0.90_{(0.00)}$ |  $0.90_{(0.00)}$ |  $0.90_{(0.00)}$ |  $0.90_{(0.00)}$
CB |  $0.90_{(0.00)}$ |  $0.91_{(0.01)}$ |  $0.88_{(0.01)}$ |  $0.90_{(0.00)}$ |  $0.90_{(0.00)}$ |  $0.89_{(0.00)}$ |  $0.90_{(0.00)}$ |  $0.90_{(0.00)}$
CHR |  $0.90_{(0.00)}$ |  $0.91_{(0.00)}$ |  $0.91_{(0.00)}$ |  $0.90_{(0.00)}$ |  $0.90_{(0.00)}$ |  $0.90_{(0.00)}$ |  $0.90_{(0.00)}$ |  $0.90_{(0.00)}$
DCP |  $0.90_{(0.00)}$ |  $0.91_{(0.00)}$ |  $0.91_{(0.00)}$ |  $1.00_{(0.00)}$ |  $1.00_{(0.00)}$ |  $1.00_{(0.00)}$ |  $1.00_{(0.00)}$ |  $0.90_{(0.00)}$
R2CCP (Ours) |  $0.90_{(0.00)}$ |  $0.90_{(0.00)}$ |  $0.90_{(0.00)}$ |  $0.90_{(0.00)}$ |  $0.90_{(0.00)}$ |  $0.90_{(0.00)}$ |  $0.90_{(0.00)}$ |  $0.90_{(0.00)}$

Dataset |  Diabetes |  Solar |  Parkinsons |  Stock |  Cancer |  Pendulum |  Energy |  Forest
---|---|---|---|---|---|---|---|---
CQR |  $0.91_{(0.00)}$ |  $0.91_{(0.00)}$ |  $0.90_{(0.00)}$ |  $0.91_{(0.00)}$ |  $0.92_{(0.00)}$ |  $0.91_{(0.00)}$ |  $0.90_{(0.00)}$ |  $0.91_{(0.00)}$
KDE |  $0.91_{(0.00)}$ |  $0.71_{(0.15)}$ |  $0.89_{(0.00)}$ |  $0.90_{(0.00)}$ |  $0.92_{(0.01)}$ |  $0.90_{(0.01)}$ |  $0.91_{(0.00)}$ |  $0.81_{(0.08)}$
Lasso |  $0.91_{(0.00)}$ |  $0.90_{(0.00)}$ |  $0.90_{(0.00)}$ |  $0.91_{(0.00)}$ |  $0.92_{(0.00)}$ |  $0.90_{(0.00)}$ |  $0.90_{(0.00)}$ |  $0.90_{(0.00)}$
CB |  $0.89_{(0.00)}$ |  $0.92_{(0.01)}$ |  $0.90_{(0.00)}$ |  $0.88_{(0.01)}$ |  $0.87_{(0.02)}$ |  $0.88_{(0.01)}$ |  $0.89_{(0.01)}$ |  $0.90_{(0.01)}$
CHR |  $0.91_{(0.00)}$ |  $0.91_{(0.00)}$ |  $0.90_{(0.00)}$ |  $0.91_{(0.00)}$ |  $0.92_{(0.00)}$ |  $0.91_{(0.00)}$ |  $0.91_{(0.00)}$ |  $0.91_{(0.00)}$
DCP |  $0.91_{(0.00)}$ |  $1.00_{(0.00)}$ |  $0.90_{(0.00)}$ |  $0.91_{(0.00)}$ |  $0.92_{(0.00)}$ |  $0.90_{(0.00)}$ |  $0.91_{(0.00)}$ |  $1.00_{(0.00)}$
R2CCP (Ours) |  $0.90_{(0.00)}$ |  $0.92_{(0.02)}$ |  $0.96_{(0.00)}$ |  $0.92_{(0.02)}$ |  $0.90_{(0.00)}$ |  $0.90_{(0.00)}$ |  $0.91_{(0.01)}$ |  $0.89_{(0.00)}$

Table 3: This is the coverage data over all datasets. We see that all methods
achieve the roughly $1-\alpha$ coverage guaranteed by conformal prediction.

Dataset |  Bimodal |  LogNorm |  Concrete |  MEPS-19 |  MEPS-20 |  MEPS-21 |  Bio |  Community
---|---|---|---|---|---|---|---|---
$\mathcal{L}$ |  $0.90_{(0.02)}$ |  $0.90_{(0.02)}$ |  $0.90_{(0.02)}$ |  $0.90_{(0.01)}$ |  $0.90_{(0.01)}$ |  $0.90_{(0.01)}$ |  $0.90_{(0.00)}$ |  $0.90_{(0.01)}$
$\mathcal{L}_{\text{NE}}$ |  $0.90_{(0.02)}$ |  $0.81_{(0.03)}$ |  $0.84_{(0.03)}$ |  $0.90_{(0.01)}$ |  $0.90_{(0.01)}$ |  $0.90_{(0.01)}$ |  $0.68_{(0.01)}$ |  $0.88_{(0.02)}$
$\mathcal{L}_{\text{MLE}}$ |  $0.90_{(0.02)}$ |  $0.90_{(0.02)}$ |  $0.90_{(0.02)}$ |  $0.90_{(0.01)}$ |  $0.90_{(0.01)}$ |  $0.90_{(0.01)}$ |  $0.90_{(0.00)}$ |  $0.90_{(0.01)}$
$\mathcal{L}_{\text{MLE + E}}$ |  $0.90_{(0.02)}$ |  $0.90_{(0.02)}$ |  $0.90_{(0.02)}$ |  $0.90_{(0.01)}$ |  $0.90_{(0.01)}$ |  $0.90_{(0.01)}$ |  $0.92_{(0.00)}$ |  $0.90_{(0.01)}$

Dataset |  Diabetes |  Solar |  Parkinsons |  Stock |  Cancer |  Pendulum |  Energy |  Forest
---|---|---|---|---|---|---|---|---
$\mathcal{L}$ |  $0.90_{(0.03)}$ |  $0.90_{(0.02)}$ |  $0.95_{(0.01)}$ |  $0.90_{(0.03)}$ |  $0.90_{(0.05)}$ |  $0.90_{(0.03)}$ |  $0.90_{(0.03)}$ |  $0.89_{(0.03)}$
$\mathcal{L}_{\text{NE}}$ |  $0.94_{(0.02)}$ |  $0.90_{(0.02)}$ |  $0.73_{(0.01)}$ |  $0.89_{(0.03)}$ |  $0.90_{(0.05)}$ |  $0.90_{(0.03)}$ |  $0.90_{(0.02)}$ |  $0.77_{(0.04)}$
$\mathcal{L}_{\text{MLE}}$ |  $0.90_{(0.03)}$ |  $0.90_{(0.02)}$ |  $0.90_{(0.01)}$ |  $0.91_{(0.03)}$ |  $0.90_{(0.05)}$ |  $0.90_{(0.03)}$ |  $0.90_{(0.02)}$ |  $0.89_{(0.03)}$
$\mathcal{L}_{\text{MLE + E}}$ |  $0.90_{(0.03)}$ |  $0.90_{(0.02)}$ |  $0.90_{(0.01)}$ |  $0.90_{(0.03)}$ |  $0.90_{(0.05)}$ |  $0.90_{(0.03)}$ |  $0.90_{(0.03)}$ |  $0.90_{(0.03)}$

Table 4: This is the coverage data for different loss functions from the
ablation.

##  Appendix F More Experimental Details

In order to maintain fairness across all baselines, we use the same size
neural network for our method across all experiments. Specifically, we
discretize the range space into $K=50$ points, weight the entropy term by
$\tau=0.2$, use a $1000$ hidden dimension, use $4$ layers, use weight decay of
$1e-4$, use $p=.5$, and use AdamW as an optimizer. For most of the
experiments, we use learning rate $1e-4$ and batch size $32$. However, for
certain datasets, namely the MEPS datasets, we used a larger batch size of
$256$ to improve training time and used a smaller learning rate to prevent
training divergence. We did change any other parameter between all of our
runs. For the baselines of CQR and CHR, we use the neural network
configurations mentioned in the paper. We found that the parameterizations
mentioned by the authors in the papers achieved the best performance and
changing the parameterizations weakened their results.

##  Appendix G Additional Experiments

We add several additional experiments to help contextualize our results.
First, we add experiments tracking how our method performs over different
coverage levels relative to CHR and CQR across several datasets in Section
G.2. Second, in Section G.3, we analyze what percentage of the time the
predicted intervals are singeltons, meaning there is only one continuous
interval. Third, we analyze the correlation between the residual prediction
error and the predicted interval’s length with our method in Section G.4.

###  G.1 Performance with different Number of Bins

We set the number of bins at a constant $K=50$. A natural question is how the
number of bins affects the performance of our conformal prediction
methodology. To test this, we retrain several models with varying number of
bins on the Concrete dataset. We report the size of the predicted intervals
for each number of bins in Figure 3. We observe that as the number of bins
increases, the average length of the predicted intervals decreases and
stabilizes. This suggests that after a sufficient number of bins, the length
of the intervals becomes relatively constant. In practice, classification
problems with $1000$ bins or more have been routinely solved since $2012$
(Krizhevsky et al., 2012). This suggests that the increasing number of bins
does not make the learning task more difficult, explaining how performance
remains relatively unaffected as the number of bins is increased.

Figure 3: We present an ablation of how the number of bins affects the average
length of the generated intervals.

###  G.2 Different Coverage Levels

While the coverage level of $\alpha=.1$ is relatively standard across the
Conformal Prediction literature, we want to see how our method compares to
other methods across different coverage levels. We take the trained models
from the original experiments and compute the intervals across different
coverage levels for our method, CHR, and CQR. As expected, for all methods,
the length of the predicted intervals increases as the required coverage level
increases. We report our results in Figure 4.

|  |
---|---|---
|  |
|  |
|  |
|  |

(a) MEPS-20

(b) MEPS-21

(c) DIABETES

(d) BIMODAL

(e) PENDULUM

(f) BIO

(g) LOG NORMAL

(h) ENERGY

(i) BREAST CANCER

(j) CONCRETE

(k) STOCK

(l) COMMUNITY

(m) FOREST

(n) PARKINSON’S

(o) MEPS-21

Figure 4: We plot the width of the predicted intervals for our method, CHR,
and CQR across different coverage levels.

###  G.3 Number of Singletons

We also provide an analysis of the percentage of times that our predicted
intervals are singleton, meaning only one interval is predicted. We present
our results in Table 5. We see that for a majority of the datasets, our model
tends to predict singleton intervals. However, for both the Bimodal dataset
and the Bio dataset, we see that our method predicts far less singletons. This
observation is intuitive for the Bimodal dataset as two intervals are expected
to account for the different modes.

Dataset |  Bimodal |  Log-Normal |  Concrete |  MEPS-19 |  MEPS-20 |  MEPS-21 |  Bio |  Community
---|---|---|---|---|---|---|---|---
Number of Singletons |  $5$ |  $197$ |  $194$ |  $3157$ |  $3509$ |  $3132$ |  $6725$ |  $382$
Number of Examples |  $400$ |  $200$ |  $206$ |  $3157$ |  $3509$ |  $3132$ |  $9146$ |  $399$

Dataset |  Diabetes |  Solar |  Parkinsons |  Stock |  Cancer |  Pendulum |  Energy |  Forest
---|---|---|---|---|---|---|---|---
Number of Singletons |  81 |  174 |  1103 |  105 |  37 |  115 |  153 |  96
Number of Examples |  89 |  214 |  1175 |  108 |  39 |  126 |  154 |  104

Table 5:  We present an analysis of the number of singletons predicted by our
method over all datasets. For most datasets, we see that our method tends to
predict singletons.

###  G.4 Correlation between Prediction Error and Length of Predicted
Intervals

We want to see how the prediction error correlates with the predicted interval
lengths. We define prediction error as the absolute error residual between
predicted label and the true label. Given our model is taylored to generated
intervals, there are many ways to predict a single label for a given set of
features. For example, we can select the bin with the largest probability
predicted by our model

$\hat{y}=\hat{y}_{l}\text{ where
}k=\underset{k\in[K]}{\operatorname*{arg\,max}}q_{\theta}(\hat{y}_{k}|x)\text{.}$

We can also use the expected label according to our distribution as the
predicted label

$\hat{y}=\sum_{k=1}^{K}\hat{y}_{k}q_{\theta}(\hat{y}_{k}|x)\text{.}$

We choose the second method for the sake of these experiments. We compute the
absolute residual error as $\|\hat{y}-y\|$. We plot our results in Figure 5.
We see that the prediction error is not strongly correlated with the length of
the predicted interval.

|  |
---|---|---
|  |
|  |
|  |
|  |

(a) MEPS-20

(b) MEPS-21

(c) DIABETES

(d) BIMODAL

(e) PENDULUM

(f) BIO

(g) LOG NORMAL

(h) ENERGY

(i) BREAST CANCER

(j) CONCRETE

(k) STOCK

(l) COMMUNITY

(m) FOREST

(n) PARKINSON’S

(o) SOLAR

Figure 5: We plot the corrleation between absolute residual error and the
length of our predicted intervals over several datasets. We see that there is
not a strong correlation.

##  Appendix H Output Distribution

Here, we plot the distribution of lengths of all Conformal Prediction methods
over all datasets. We also plot example density functions learned by out
method, CHR, and KDE on all datasets. We note that KDE’s learned densities
seem to be relatively less informative. Moreover, the learned densities from
CHR are noisy and not smooth. Moreover, we plot several examples of only our
learned probability distribution. We use this when referencing the experiments
to demonstrate the different shapes of label distributions from the data.

|  |
---|---|---
|  |
|  |
|  |
|  |

(a) MEPS-20

(b) MEPS-21

(c) DIABETES

(d) BIMODAL

(e) PENDULUM

(f) BIO

(g) LOG NORMAL

(h) ENERGY

(i) BREAST CANCER

(j) CONCRETE

(k) STOCK

(l) COMMUNITY

(m) FOREST

(n) PARKINSON’S

(o) SOLAR

Figure 6: Example Probability Distributions outputted by our learned
methodology on different datapoints from different datasets. We can see there
are a variety of different shapes of distributions learned.

(a) Bimodal Length Distribution (b) Bio Length Distribution

(c) Breast Cancer Length Distribution (d) Community Length Distribution

(e) Concrete Length Distribution (f) Diabetes Length Distribution

(g) Energy Length Distribution (h) Forest Length Distribution

(i) Log Normal Length Distribution (j) MEPS 19 Length Distribution

(k) MEPS 20 Length Distribution (l) MEPS 21 Length Distribution

(m) Parkinsons Length Distribution (n) Pendulum Length Distribution

(o) Solar Length Distribution (p) Stock Length Distribution

(q) BIO

(r) BREAST CANCER

Figure 7: Plot of outputted density functions for ours, KDE, and CHR on for
BIO and BREAST CANCER

(a) BIMODAL

(b) COMMUNITY

Figure 8: Plot of outputted density functions for ours, KDE, and CHR on for
BIMODAL and COMMUNITY

(a) CONCRETE

(b) DIABETES

Figure 9: Plot of outputted density functions for ours, KDE, and CHR on for
CONCRETE and DIABETES

(a) ENERGY

(b) FOREST

Figure 10: Plot of outputted density functions for ours, KDE, and CHR on for
ENERGY and FOREST

(a) LOG NORMAL

(b) MEPS-21

Figure 11: Plot of outputted density functions for ours, KDE, and CHR on for
LOG NORMAL and MEPS-19

(a) PARKINSON’S

(b) PENDULUM

Figure 12: Plot of outputted density functions for ours, KDE, and CHR on for
PARKINSON’S and PENDULUM

(a) SOLAR

(b) STOCK

Figure 13: Plot of outputted density functions for ours, KDE, and CHR on for
SOLAR and STOCK


Integrals of tau functions I: Round dance tau function and multi-matrix
integrals

A. Yu<EMAIL_ADDRESS>

1Shirshov Institute of Oceanology, Russian Academy of Science, Nahimovskii
Prospekt 36, Moscow 117997, Russia
2A.I. Alikhanov Institute for Theoretical and Experimental Physics of NRC
Kurchatov Institute, B. Cheremushkinskaya, 25, Moscow, 117259, Russia

###### Abstract

The simplest nontrivial tau functions of the Toda lattice and the N-component
Toda lattice are compared in their applications to multimatrix integrals.

##  1 Introduction

I want to dedicate this article to Andrei Pogrebkov in connection with his
75th birthday. Andrei has made notable and varied contributions to the
development of the theory of integrable systems - the vertex operator (the
work by Pogrebkov and Sushko 1974, late published [1]), correctly defined
Hamiltonian formalism, the inverse problem of scattering theory, a detailed
study of a number of specific integrable equations, and a number of other
often very original approaches to integrability theory. The text below
represents a slightly different direction, but I am grateful for the always
inspiring discussions with Andrey on this topic.

##  2 The simplest nontrivial tau functions of the multicomponent TL

The well-known Cauchy-Littlewood relation can be written as

$e^{-\sum_{m>0}\frac{1}{m}p_{m}\tilde{p}_{m}}=1+\sum_{\kappa>0}\sum_{\alpha_{1}>\cdots>\alpha_{\kappa}\atop\beta_{1}>\cdots>\beta_{\kappa}}s_{(\alpha|\beta)}(\mathbf{p})s_{(\beta|\alpha)}(\tilde{\bf{p}})=\tau_{o}(\mathbf{p},\tilde{\bf{p}})$
(2.1)

Here $\mathbf{p}=(p_{1},p_{2},\dots)$ and
$\tilde{\mathbf{p}}=(\tilde{p}_{1},\tilde{p}_{2},\dots)$ are two infinite sets
of variables and $s_{(\alpha|\beta)}$ denotes the Schur polynomial [2] in the
power sum variables labeled by a partition $(\alpha|\beta)$ written in the
Frobenius coordinates $\alpha=(\alpha_{1},\dots,\alpha_{\kappa})$ and
$\beta=(\beta_{1},\dots,\beta_{\kappa})$. This is also the simplest example of
the Toda lattice (TL) [3] tau functions [4],[5] where the variables
$\mathbf{p}$ and $\tilde{\bf{p}}$ play the role of the so-called higher times.
111For higher times, other designations are sometimes used, for example,
$nx_{n}$ instead of $p_{n}$, the sign can be also different. To fix the
correspondence we write the Appendix. It is written in form of series in the
Schur functions as it was worked out in [6], [7]. The straightforward
generalization of this expression to the N-component Toda lattice (N-TL)
yields

$\tau_{o}(\mathbf{p}^{1},\tilde{\bf{p}}^{1},\dots,\mathbf{p}^{\textsc{N}},\tilde{\bf{p}}^{\textsc{N}})=1+\sum_{\kappa>0}\sum_{\alpha^{1}_{1}>\cdots>\alpha^{1}_{\kappa},\atop\beta^{1}_{1}>\cdots>\beta^{1}_{\kappa}}\cdots\sum_{\alpha^{\textsc{N}}_{1}>\cdots>\alpha^{\textsc{N}}_{\kappa},\atop\beta^{\textsc{N}}_{1}>\cdots>\beta^{\textsc{N}}_{\kappa}}\,\prod_{1\leq
i\leq\textsc{N}}^{\circlearrowleft}s_{(\alpha^{i}|\beta^{i})}(\mathbf{p}^{i})s_{(\beta^{i}|\alpha^{i+1})}(\tilde{\bf{p}}^{i})$
(2.2)

where each set of the higher times is denoted by $\mathbf{p}^{i}$ and
$\tilde{\bf{p}}^{i}$, $i=1,\dots,\textsc{N}$. The symbol $\circlearrowleft$
above the symbol of the product means that we suppose the condition
$\alpha^{\textsc{N}+1}=\alpha^{1}$, which closes the chain of Schur functions
in a circle which explains the motivation of the name “round dance”. I like
this expression and it somehow reminds me of Matisse’s painting “Dance”. There
is a beautiful expression for series (2.2), similar to the left side in (2.1),
see formula (38) in [8].

##  3 Graphs and matrix models

Consider the set of complex $GL_{N}$-matrices $Z_{1},\dots,Z_{n}$ and their
Hermitian conjugates $Z_{-1},\dots,Z_{-n}$, where $Z_{-i}:=Z_{i}^{\dagger}$.
We multiply each $Z_{i}$ by $C_{i}\in GL_{N}$ on the right,
$Z_{i}\to(Z_{i}C_{i})$, where $i=\pm 1,\dots,\pm n$. The matrices $Z_{i}$ will
be called random, and the matrices $C_{i}$ will be called source ones. The
elements of random matrices are distributed according to the Gauss law as
follows:

$\langle(Z_{a})_{ij}(Z_{a^{\prime}})_{i^{\prime}j^{\prime}}\rangle=\frac{1}{N}\delta_{a+a^{\prime},0}\delta_{ii^{\prime}}\delta_{jj^{\prime}}$
(3.1)

Such correlation functions are obtained for the so-called Ginibre complex
multi-matrix ensembles [9],[11],[10]. To calculate various correlation
functions, we use Wick’s rule. We will study special correlation functions,
which are obtained as follows:

We denote the collection of random matrices by Z. We consider the set of
$GL_{N}$ matrices

$W_{1}(\textsc{Z}),\dots,W_{\textsc{v}}(\textsc{Z})$

where each $W_{a}(\textsc{Z})$ is the product of $Z_{i}C_{i}$ pairs. Each
product $W_{a}(\textsc{Z})$ is defined up to a cyclic permutation. We set the
condition: each pair enters one and only once to one of the members of the set
$\\{W_{a}(\textsc{Z}),a=1,\dots,\textsc{v}\\}$. In this case, there is a two-
dimensional orientable surface $\Omega$ and an embedded graph222For embedded
graphs, there are other names: ”thick graphs” and ”ribbon graphs” in physical
literature and ”maps” in mathematical literature with $n$ numbered edges and
with f numbered vertices, drawn on $\Omega$ as follows. Each edge labeled
$|i|$ ($|i|=1,\dots,n$) consists of two half-edges with numbers $i$ and $-i$
and associate the half-edge $i$ with the matrix $Z_{i}$. Each vertex labeled
$a$ is connected to $W_{a}(Z)$ as follows: if we move clockwise, then the
outgoing half-edges will have the same numbers as the random matrices that are
included in the product $W_{a}(\textsc{Z})$, if we read it from left to right.
Then we number the corners of the faces of the graph as follows: when
traversing a vertex clockwise, each half-edge with the number $i$ ($i=\pm
1,\dots,\pm n$) is followed by an angle with the same number. Then in corner
$i$ we place $C_{i}$. We will use the notation $W_{a}(\textsc{I})$ which
denote $W_{a}(\textsc{Z})$ where we put all matrices $Z_{i}$ to be the
identity matrix. We see that $W_{a}(\textsc{I})$ can be called the monodromy
of the vertex $a$ and $W_{a}(\textsc{Z})$ the monodromy dressed by random
matrices. Suppose that the Euler characteristic of $\Omega$ is e,
$\textsc{e}=\textsc{f}-n+\textsc{v}$, here f is the number of faces of
$\Gamma$. Now we will number all the faces of $\Gamma$. We also assign a
monodromy to each face which is a matrix defined up to the cyclic order.
Monodromy $W_{b}^{*}(\textsc{I})$ of the face with number $b$ is introduced as
the product of all source matrices $C_{i}$ at the corners of the face in the
order in which we encounter them when going around the face counterclockwise.
The dressed monodromy of the face $W_{a}(\textsc{Z})$ is obtained from
$W_{a}(\textsc{I})$, when we replace each matrix $C_{i}$ that enters
$W_{a}(I)$ by $Z_{i}C_{i}$.

It can be seen that the dressed monodromies of faces of $\Gamma$ is the
dressed monodromy of the vertices of the dual graph $\Gamma^{*}$. This
construction was analyzed in the works [13],[14],[15],[16].

### Integrals of products of the Schur functions

For any sets of partitions $\\{\lambda^{a},a=1,\dots,\textsc{v}\\}$ and
$\\{\lambda^{b},b=1,\dots,\textsc{v}$} we have

$\langle\prod_{a=1}^{\textsc{v}}s_{\lambda^{a}}(W_{a}(\textsc{Z})\rangle=c\delta_{\lambda}\prod_{b=1}^{\textsc{f}}s_{\lambda}(W_{b}^{*}(\textsc{I}))$
(3.2)
$\langle\prod_{b=1}^{\textsc{f}}s_{\lambda^{b}}(W^{*}_{b}(\textsc{Z}))\rangle=c\delta_{\lambda}\prod_{a=1}^{\textsc{v}}s_{\lambda}(W_{a}(\textsc{I}))$
(3.3) $c=\left(\frac{\dim\lambda}{|\lambda|!}\right)^{-n}N^{-n|\lambda|}$
(3.4)

where $\dim\lambda$ is the dimension of the irreducable representation of the
symmetric group $S_{d}$ ($d=|\lambda|!$) labeled by $\lambda$. The symbol
$\delta_{\lambda}$ is equal to $1$ if all partitions are equal to each other:
$\lambda^{1}=\lambda^{2}=\cdots=:\lambda$, otherwise $\delta_{\lambda}=0$.
Notice the diagonalization property for the products of the Schur function
when we evaluate the expectation $\langle\rangle$.

######  Remark 3.1.

These nice relations can be compared with the relations for embedded graphs
[12]

$\sigma\alpha=\varphi$ (3.5) $\sigma\varphi=\alpha$ (3.6)

where, $\sigma,\alpha,\varphi\in S_{2n}$ and $\sigma$ is the involution
without fix points, $\alpha$ is the product of cycles related to vertices and
$\varphi$ is the product of cycles

For a given matrix $X$ and a partition $\mu=(\mu_{1},\dots,\mu_{\ell})$ we
introduce the notation

$\mathbf{p}_{\mu}(X)=\mathrm{tr}(X)^{\mu_{1}}\cdots\mathrm{tr}(X)^{\mu_{\ell}}$
(3.7)

One can associate $\mathbf{p}_{\mu}(X)$ with a set of polygons consisting of
$\mu_{1}$-gons,…,$\mu_{\ell}$-gons. Suppose that an $i$-gon is included in
this set $m_{i}$ times. Introduce number
$\texttt{Aut}(\mu):=\prod_{i>0}m_{i}!i^{m_{i}}$, which can be called the
number of automorphisms of this collection: number of permutations polygons
with the same number of edges multiplied by the number of rotations of each
polygon.

One can derive (3.2)-(3.3) in different ways, in particular [13], starting
from the relation

$\langle\prod_{b=1}^{\textsc{f}}\frac{\mathbf{p}_{\mu^{b}}(W^{*}_{b}(\textsc{Z}))}{|\texttt{Aut}(\mu^{b})|}\rangle=\sum_{\nu^{1},\dots,\nu^{\textsc{v}}}{\cal{H}}_{\Omega}\left(\mu^{1},\dots,\mu^{\textsc{f}},\nu^{1},\dots,\nu^{\textsc{v}}\right)\prod_{a=1}^{\textsc{v}}\mathbf{p}_{\nu^{a}}(W_{a}(\textsc{I}))$
(3.8)

where all partitions
$\mu^{1},\dots,\mu^{\textsc{f}},\nu^{1},\dots,\nu^{\textsc{v}}$ are of the
same weight, say $d$, $d\leq N$. Here ${\cal{H}}_{\Omega}$ is the Hurwitz
number, $\Omega$ the base surface and partitions
$\mu^{1},\dots,\mu^{\textsc{f}},\nu^{1},\dots,\nu^{\textsc{v}}$ are the
ramification profiles.

The meaning of the formula (3.8) is as follows. For all partitions
$\mu^{1},\nu^{1},\cdots$ equal to the partition $(1)$, this formula describes
the base surface, glued from polygons associated with
$\\{\mathrm{tr}(W^{*}_{b}(\textsc{Z})),b=1,\dots,\textsc{f}\\}$, see (3.6).
Expression for $\mathbf{p}_{\mu}$ corresponds to the set of polygons on the
covering surface, which correspond to the set of profiles $\mu^{b}$ in the
center (in the ”capital”) of the polygon $b$. Wick’s rule gives all possible
gluings of the collection of polygons given by $\mathbf{p}_{\nu}$. Each way of
gluing leads to a set $\\{\nu^{a},a=1,\dots,\textsc{v}\\}$ of branching
profiles at the vertices. The number of ways to glue polygons (up to the
automorphisms) for given sets of all ramification profiles
$\mu^{1},\dots,\mu^{\textsc{f}},\nu^{1},\dots,\nu^{\textsc{v}}$ is the Hurwitz
number in its geometric meaning.

Then from the Frobenius formula for the Hurwitz numbers and from the
orthogonality relations of the characters symmetric group, we get (3.3) and
(3.2). Another way to get (3.3) and (3.2) - consistent use of well-known
formulas for integrating Schur functions (using analytic continuation of
parameters), see for instance [15].

If instead of the Ginibra ensemble we consider ensembles of unitary matrices
(circular ensembles), then the relations (3.2) and (3.3) will change slightly,
in particular, the dimension of the representation of the symmetric group will
have to be replaced by the dimension of the representation of the linear
group. On the related topics see also [17] and [18].

Actually we need only formula (3.2) for the next part.

##  4 Round dance tay function and multimatrix Ginibre ensemble

Now consider the construction of the previous section and in formala (2.2)
(where we replace N by v to avoid a mess in the usage of the capital N) we
choose the times $\tilde{\bf{p}}^{a}$ as follows

$\tilde{p}^{(a)}_{k}=\mathrm{tr}(W_{a}(\textsc{Z}))^{k},\quad
a=1,\dots,\textsc{v}$ (4.1)

We get

$\langle\tau(\mathbf{p}^{1},\tilde{\bf{p}}^{1}),\dots,\mathbf{p}^{\textsc{v}},\tilde{\bf{p}}^{\textsc{v}})\rangle=\sum_{\lambda}\left(\frac{\dim\lambda}{|\lambda|!}\right)^{-n}N^{-n|\lambda|}\prod_{a=1}^{\textsc{v}}s_{\lambda}(\mathbf{p}^{a})\prod_{b=1}^{\textsc{f}}s_{\lambda}(W_{b}^{*}(\textsc{I}))$
(4.2)

The series (4.2) appeared in the works [19],[20],[21],[22], [23],[24] as a
generating function for Hurwitz numbers. To achieve these series we chose the
condition (4.1).

Note that, apparently, the right-hand side can only be a tau function of the
hypergeometric type of the TL hierarchy [25],[26],[27],[28]. It’s not an exact
statement, but I don’t see any other possibilities at the moment. That is,
there may be a tau function of the form

$\sum_{\lambda}s_{\lambda}(\mathbf{p}^{1})s_{\lambda}(\mathbf{p}^{2})\prod_{(i,j)\in\lambda}r(j-i)$
(4.3)

with some function $r$.

What are our chances of getting it. First, the Euler characteristic of the
base surface $\Omega$ must be equal to two, therefore we have a graph $\Gamma$
on a Riemann sphere. Then only two sets of free parameters should remain.
Let’s look at all the possibilities, there are three of them:

(1) two sets among $\mathbf{p}^{a}$ are free, all others have a form

$\mathbf{p}^{i}=\mathbf{p}(a_{i}):=(a_{i},a_{i},\dots),\quad i=1,\dots,k$
(4.4)

(in the considered case $k\leq\textsc{v}$) or coinside with

$\mathbf{p}^{i}=\mathbf{p}_{\infty}:=(1,0,0,\dots),\quad
i=1,\dots,\textsc{v}-k$ (4.5)

where $a$ is a complex number. In this case the spectrum of matrices
$W^{*}_{b}$, $b=1,\dots,\textsc{f}$ should be

$\texttt{Spectr}\left[W_{b}(\textsc{I})\right]=(\underbrace{1,\dots,1}_{N_{b}},\underbrace{0,\dots,0}_{N-N_{b}})$
(4.6)

with some $N_{b}$. In this case the right hand side of (4.2) is as follows

$\sum_{\lambda}N^{-n|\lambda|}s_{\lambda}(\mathbf{p}^{1})s_{\lambda}(\mathbf{p}^{2})\prod_{(i,j)\in\lambda}\prod_{c=1}^{k}(a_{c}+j-i)\prod_{b=1}^{\textsc{v}}(N_{b}+j-i)$
(4.7)

About specializations (4.4),(4.5) see [27]), [17], and see [15],[16] on (4.6).

Among these tau functions, one can be found, which was studied in [29]. The
authors use rectangular matrices $Z$, in our approach this is the choice of
source matrices which may degenerated. The graph $\Gamma^{*}$ in this case is
just an open chain.

Tau functions of type (4.7) for the first time were used in [30], then in [31]
in the context of combinatorial problems and also in [32], [33],[34] as the
generating functions of Hurwitz numbers of certain type.

(2) We have one free set, say, $\mathbf{p}^{1}$, all others have forms (4.5)
and (4.4). We have a single free monodromy, say, $W^{*}_{1}$, all others have
their spectrum as in (4.6). In this case the right hand side of (4.2) is

$\sum_{\lambda}N^{-n|\lambda|}s_{\lambda}(\mathbf{p}^{1})s_{\lambda}(\mathbf{p}^{2}(W^{*}_{1}(\textsc{I})))\prod_{(i,j)\in\lambda}\prod_{c=1}^{k+1}(a_{c}+j-i)\prod_{b=2}^{\textsc{v}}(N_{b}+j-i)$
(4.8)

It is written as a function of the set of the higher times $\mathbf{p}^{1}$
and as a function in the so-called Miwa variables given by the spectrum of
$W^{*}_{1}(\textsc{I})$

(3) We two free monodromies, say, these are $W^{*}_{1}(\textsc{I}))$ and
$W^{*}_{2}(\textsc{I}))$ while the spectrum of all others are given by (4.6).
Each $\mathbf{p}^{a}$ is either (4.4) or (4.5). In this case we have the right
hand side of (4.2) written in Miwa variables as

$\sum_{\lambda}N^{-n|\lambda|}s_{\lambda}(\mathbf{p}^{1}(W^{*}_{1}(\textsc{I})))s_{\lambda}(\mathbf{p}^{2}(W^{*}_{2}(\textsc{I})))\prod_{(i,j)\in\lambda}\prod_{c=1}^{k+2}(a_{c}+j-i)\prod_{b=3}^{\textsc{v}}(N_{b}+j-i)$
(4.9)

Acknowledgements. I would like to thank Andrei D. Mironov for the inspiring
discussion. The work was supported by the Russian Science Foundation (Grant
No.20-12-00195).

## References

  * [1] A. K. Pogrebkov and V. N. Sushko, Quantization of the $(sin\psi)_{2}$ interaction in terms of fermion variables, Translated from Teoreticheskaya i Mathematicheskaya Fizika, 24 (1975) 425-429 (September, 1975). The original paper was submitted on May 15, 1975
  * [2] I.G. Macdonald, Symmetric Functions and Hall Polynomials, Clarendon Press, Oxford, (1995).
  * [3] A.V. Mikhailov, “On the Integrability of two-dimensional Generalization of the Toda Lattice”, Letters in Journal of Experimental and Theoretical Physic s, v.30, p. 443-448, 1979; A.V.Mikhailov, M. A. Olshanetski, A.M.Perelomov, Two-dimentional generalized Toda lattice, Comm.Math.Phys 79 (1981) no. 4 473-488
  * [4] M. Jimbo and T. Miwa, “Solitons and infinite dimensional Lie algebras”, Publ. RIMS Kyoto Univ. 19, pp. 943–1001 (1983)
  * [5] K.Ueno and K.Takasaki, Toda lattice hierarchy, Adv. Stud. Pure Math. 4, 1-95 (1984).
  * [6] K. Takasaki, Initial value problem for the Toda lattice hierarchy, Adv. Stud. Pure Math. 4 (1984) 139–163
  * [7] T. Takebe, Representation Theoretical Meaning of Initial Value Problem for the Toda Lattice Hierarchy I, LMP 21 (1991) 77–84
  * [8] A. Mironov, A. Morozov, Hook variables: cut-and-join operators and tau functions, Phys.Lett. B804 (2020) 135362; arxiv:1912.00635
  * [9] M. L. Mehta “Random Matrices”, 3nd edition (Elsevier, Academic, San Diego CA, 2004)
  * [10] E. Strahov, Differential equations for singular values of products of Ginibre random matrices, Journal of Physics A Mathematical and Theoretical 47(32) (2014)
  * [11] Ipsen, J.R. Products of Independent Gaussian Random Matrices. Doctoral dissertation, arXiv:1510.06128
  * [12] S. K. Lando, A. K. Zvonkin, “Graphs on Surfaces and their Applications”, Encyclopaedia of Mathematical Sciences, Volume 141, with appendix by D. Zagier, Springer, N.Y. (2004)
  * [13] S. M. Natanzon and A. Yu. Orlov, Hurwitz numbers from matrix integrals over Gaussian measure, Proc.Symposium in Pure Math. 103 1 (2021) Integrability, Quantization and Geometry, pp 337-377; arXiv:2002.00466.
  * [14] S. M. Natanzon and A. Yu. Orlov,“Hurwitz numbers from Feynman diagrams”, Theoretical and Mathematical Physics 204 (3), (2020) 1166-1194; arXiv:2006.07396
  * [15] A.Yu.Orlov, 2D Yang-Mills Theory and Tau Functions in: Geometric Methods in Physics XXXVIII, Editors: P. Kielanowski, Anatol Odzijewicz, E. Previato 221-250 (2020)
  * [16] N Amburg, A Orlov, D Vasiliev, On products of random matrices, Entropy 22 (9), 972 (2020)
  * [17] A. Yu. Orlov, New solvable matrix integrals , Intern. J. Mod. Phys. A 19  (suppl 02), pp. 276-93 (2004).
  * [18] Alexandrov, A.: Matrix models for random partitions. Nucl. Phys. B 851, 620-650 (2011)
  * [19] A. D. Mironov, A. Yu. Morozov and S. M. Natanzon, Complect set of cut-and-join operators in the Hurwitz-Kontsevich theory, Theor. and Math.Phys. 166:1,(2011), pp.1-22; arXiv:0904.4227
  * [20] A. D. Mironov, A. Yu. Morozov and S. M. Natanzon, Integrability properties of Hurwitz partition functions. II. Multiplication of cut-and-join operators and WDVV equations, JHEP 11 (2011) 097
  * [21] A. D. Mironov, A. Yu. Morozov and S. M. Natanzon, Integrability of Hurwitz Partition Functions. I. Summary, J. Phys. A: Math. Theor. 45 (2012) 045209
  * [22] A. D. Mironov, A. Yu. Morozov and S. M. Natanzon, Algebra of differential operators associated with Young diagramms, J.Geom.and Phys. n.62(2012), pp. 148-155
  * [23] A. Alexandrov, A. Mironov, A. Morozov and S. Natanzon, Integrability of Hurwitz Partition Functions. I. Summary, J.Phys.A: Math.Theor.45(2012) 045209, arXiv: 1103.4100
  * [24] A. Alexandrov, A. Mironov, A. Morozov and S. Natanzon, On KP-integrable Hurwitz functions,JHEP 11(2014) 080, arXiv: 1405.1395
  * [25] S. Kharchev, A. Marshakov, A. Mironov and A. Morozov, Generalized Kazakov-Migdal-Kontsevich Model: group theory aspects, International Journal of Mod Phys A10 (1995) p.2015
  * [26] S. Kharchev, “Kadomtsev-Petviashvili Hierarchy and Generalized Kontsevich Model”, arXiv:hep-th/9810091
  * [27] A. Yu. Orlov and D. Scherbin, Fermionic representation for basic hypergeometric functions related to Schur polynomials, arXiv preprint nlin/0001001
  * [28] A. Yu. Orlov and D. Scherbin, Hypergeometric solutions of soliton equations, Theoretical and Mathematical Physics 128 (1), pp. 906-926 (2001)
  * [29] J. Ambjorn and L. Chekhov The matrix model for hypergeometric Hurwitz number, Theoret. and Math. Phys., 1 81:3 (2014), 1486-1498; arXiv:1409.3553
  * [30] I. P. Goulden and D. M. Jackson, The KP hierarchy, branched covers, and triangulations, Advances in Mathematics, 219 pp. 932-951, 2008
  * [31] M. Kazarian and S. Lando, Combinatorial solutions to integrable hierarchies, Uspekhi Mat. Nauk 70 (2015), no. 3(423), pp. 77-106. English translation: 2015 Russ. Math. Surv. 70, pp. 453-482; arXiv:1512.07172
  * [32] M. Guay-Paquet, J. Harnad, 2D Toda $\tau$-functions as combinatorial generating functions, Lett. Math. Phys. 105, 827-852 (2015)
  * [33] J. Harnad and A. Yu. Orlov, “Hypergeometric $\tau$-functions, Hurwitz numbers and enumeration of paths”, Commun. Math. Phys. 338 (2015) pp. 267–284, arxiv: math.ph/1407.7800
  * [34] S.M. Natanzon, A.Yu. Orlov, BKP and projective Hurwitz numbers, Letters in Mathematical Physics 107 (6), 1065-1109 (2017)

##  Appendix A Partitions and Schur functions

The partition is a set of integers
$\lambda=(\lambda_{1},\dots,\lambda_{\ell})$, where
$\lambda_{1}\geq\cdots\lambda_{\ell}>0$ and $\lambda_{i}$ are called parts of
$\lambda$. The sum of the parts is called the weight of $\lambda$ and is
usually denoted as $|\lambda|$. There is a well-known notion of Young diagrams
$\lambda$ whose row lengths are equal to parts of $\lambda$. We accept the
agreement as in [2]. Let us denote the horizontal lengths of the part of the
strings that starts at the $ii$ node (and does not include it) and goes to the
right as $\alpha_{i}$, and the vertical part of the column that starts at (not
including) the same node and goes down we denote $\beta_{i}$. In Frobenius
coordinates, $\lambda$ is written as
$(\alpha|\beta)=(\alpha_{1},\dots,\alpha_{k}|\beta_{1},\dots,\beta_{k})$,
where $k$ is the number of nodes on the main diagonal of the Young diagram,
see [2] for details.

The following relations define Schur polynomials (Schur functions) as a
function of power sum variables $\mathbf{p}=(p_{1},p_{2},\dots)$:

$e^{\sum_{k>0}\frac{1}{k}p_{k}x^{k}}=\sum_{k\geq 0}s_{(k)}(\mathbf{p})$
$s_{\lambda}(\mathbf{p})=\mathrm{det}\left[s_{(\lambda_{i}-i+j)}(\mathbf{p})\right]_{i,j}$
${\rm In\,case}\quad p_{m}=\mathrm{tr}(X)^{m}\quad{\rm one\,can\,write}\quad
s_{\lambda}(\mathbf{p}(X))=:s_{\lambda}(X)$

We have

$s_{\lambda}(\mathbf{p})=\frac{\dim\lambda}{|\lambda|!}\sum_{\mu}\varphi_{\lambda}(\mu)\mathbf{p}_{\mu},\quad\dim\lambda=\frac{\prod_{i<j}(\lambda_{i}-\lambda_{j}-i+j)}{\prod_{i}(\lambda_{i}-i+\ell)!}$

The known orthogonality relations for the factors $\varphi_{\lambda}$ (which
are the normalized characters of the symmetric group in representation
$\lambda$) [2] are as follows:

$\zeta_{\Delta}\sum_{\lambda}\left(\frac{{\rm\dim}\lambda}{d!}\right)^{2}\varphi_{\lambda}(\mu)\varphi_{\lambda}(\Delta)=\delta_{\Delta,\mu}$
(A.1)

and

$\left(\frac{{\rm\dim}\lambda}{d!}\right)^{2}\sum_{\Delta}\zeta_{\Delta}\varphi_{\lambda}(\Delta)\varphi_{\mu}(\Delta)=\delta_{\lambda,\mu}$
(A.2)

where $d=|\Delta|=|\lambda|$ and if each part $i$ enters $m_{i}$ times the
partition $\Delta$ we define

$\zeta_{\Delta}=\prod_{i}m_{i}!i^{m_{i}}$ (A.3)

(The same number was denoted $|\texttt{Aut}(\Delta)|$ in the paragraph
”Integrals of products of the Schur functions”). All details can be found in
[2].

The Frobenius formula for Hurwitz numbers mentions in Section 3 is

${\cal{H}}_{\Omega}(\Delta^{1},\dots,\Delta^{m})=\sum_{\lambda}\left(\frac{\dim\lambda}{|\lambda|!}\right)^{\textsc{e}}\varphi_{\lambda}(\Delta^{1})\cdots\varphi_{\lambda}(\Delta^{m})$

where e is the Euler characteristic of the base surface $\Omega$, the
partitions $\Delta^{1},\dots,\Delta^{m}$ is the set of the ramification
profiles in branch points, see the textbook [12] for the details.

##  Appendix B Fermions. Multicomponent tau function

I will only briefly give the notation and some formulas. Fermion operators
$\psi^{(a)}_{n},\psi^{{\dagger}(a)}_{n}$ where $a=1,\dots,\textsc{N}$,
$n\in\mathbb{Z}$, satisfy the relations

$[\psi^{(a)}_{n},\psi^{(b)}_{m}]_{+}=0=[\psi^{(a){\dagger}}_{n},\psi^{(b{\dagger})}_{m}]_{+},\quad[\psi^{(a){\dagger}}_{n},\psi^{(b)}_{m}]_{+}=\delta_{a,b}\delta_{n,m}$

and

$\langle 0|\psi^{(a)}_{n}=\langle
0|\psi^{{\dagger}(a)}_{-n-1}=0=\psi^{(a)}_{-n-1}|0\rangle=\psi^{{\dagger}(a)}_{n}|0\rangle,\quad
n>0$

The wonderful observation of Kyoto school (see for instance [4]) is the
following relation

$s_{(\alpha|\beta)}(\mathbf{p}^{a})=(-1)^{\sum_{i=1}^{k}\beta_{i}}\langle
0|e^{\sum_{k>0}\frac{1}{k}p^{(a)}_{k}J_{k}^{(a)}}\sum_{i=1}^{k}\psi_{\alpha_{i}}^{(a)}\psi^{{\dagger}(a)}_{-\beta_{i}-1}|0\rangle,\quad$
(B.1)

which is also known as bosonization formula (a version of this formula). Here

$J_{n}^{(a)}=\sum_{i\in\mathbb{Z}}\psi^{(a)}_{i}\psi^{{\dagger}(a)}_{n+i}$

. We also need the property
$s_{(\alpha|\beta)}(\mathbf{p})=(-1)^{|(\alpha|\beta))|}s_{(\beta|\alpha)}(-\mathbf{p})$.

Following [4] we introduce $2\textsc{N}$-component KP (N-component Toda) tau
function in form of a fermionic expectation value which we choose as follows:

$\langle
0|e^{\sum_{a}\sum_{k>0}(-1)^{a+1}\frac{1}{k}p_{k}J^{(a)}_{k}}e^{\sum_{1\leq
a\leq
2\textsc{N}}^{\circlearrowleft}\sum_{i}\psi^{{\dagger}(a)}_{-i-1}\psi^{(a+1)}_{i}}|0\rangle=:\tau_{o}(\mathbf{p}^{1},\dots,\mathbf{p}^{\textsc{N}})$
(B.2)

where $\circlearrowleft$ means that
$\psi^{(2\textsc{N}+1)}_{-i-1}:=\psi^{(1)}_{-i-1}$. We call it round dance tau
function. Expanding the exponential in the last expression in Taylor series
and using the relations above, we get (2.2). In case $\textsc{N}=1$ we get the
right hand side of (2.1).

##  Appendix C Complex Ginibre ensemble

On this subject there is an extensive literature, for instance see lists of
references in [10, 11, 18, 29, 13, 14, 16].

Let us consider integrals over $N\times N$ complex matrices
$Z_{1},\dots,Z_{n}$ where the measure is defined as

$d\Omega(Z_{1},\dots,Z_{n})=c_{N}^{n}\prod_{i=1}^{n}\prod_{a,b=1}^{N}d\Re(Z_{i})_{ab}d\Im(Z_{i})_{ab}\text{e}^{-N|(Z_{i})_{ab}|^{2}}$
(C.1)

where the integration domain is
$\mathbb{C}^{N^{2}}\times\cdots\times\mathbb{C}^{N^{2}}$ and where $c_{N}^{n}$
is the normalization constant defined via $\int d\Omega(Z_{1},\dots,Z_{n})=1$.
We get (3.1) for the correlation functions of the entries of the matrices
$Z_{1},\dots,Z_{n}$.

The set of $n$ $N\times N$ complex matrices with measure (C.1) is called the
set of $n$ independent complex Ginibre ensembles. Such ensembles have wide
applications in physics and in information transfer theory.